1 logical function bas_version() 2c $Id$ 3c 4c: Routine that calclulates the size of the common block structures 5c used in the basis set object and the mapped representation object. 6c:input none 7c:output always true. 8c 9 implicit none 10#include "nwc_const.fh" 11#include "basP.fh" 12#include "geobasmapP.fh" 13#include "stdio.fh" 14c 15 integer cdata,idata,rdata 16 integer mapidata, total4, total8 17c 18c.. character data 19 cdata = 256*nbasis_bsmx ! bs_name 20 cdata = cdata + 16*ntags_bsmx*nbasis_bsmx ! bs_tags 21 cdata = cdata + 80*ntags_bsmx*nbasis_bsmx ! bs_stdname 22 cdata = cdata + 256*nbasis_bsmx ! bs_trans 23 cdata = cdata + 256*nbasis_rtdb_mx ! bs_names_rtdb 24 cdata = cdata + 256*nbasis_bsmx*nbasis_assoc_max ! name_assoc 25c 26c.. real data 27 rdata = 1 ! bsversion 28 rdata = 8*rdata 29c 30c.. integer data in basis set object common 31 idata = 3*nbasis_bsmx ! exndcf 32 idata = idata + ndbs_head*nbasis_bsmx ! infbs_head 33 idata = idata + 34 & ndbs_tags*ntags_bsmx*nbasis_bsmx ! infbs_tags 35 idata = idata + 36 & ndbs_ucont*(nucont_bsmx+1)*nbasis_bsmx ! infbs_cont 37 idata = idata + nbasis_bsmx ! len_bs_name 38 idata = idata + nbasis_bsmx ! len_bs_trans 39 idata = idata + nbasis_rtdb_mx ! len_bs_rtdb 40 idata = idata + nbasis_bsmx ! bsactive (int==logical) 41 idata = idata + nbasis_bsmx ! bas_spherical (int==logical) 42 idata = idata + nbasis_bsmx ! bas_any_gc (int==logical) 43 idata = idata + nbasis_bsmx ! bas_any_sp_shell (int==logical) 44 idata = idata + nbasis_bsmx ! bas_norm_id 45 idata = idata + nbasis_bsmx ! angular_bs 46 idata = idata + nbasis_bsmx ! nbfmax_bs 47 idata = idata + nbasis_assoc_max*nbasis_bsmx ! handle_assoc 48 idata = idata + nbasis_assoc_max*nbasis_bsmx ! parent_assoc 49 idata = idata + 1 ! nbasis_rtdb 50 idata = 4*idata 51c 52c.. integer data in the mapped object. 53 mapidata = 3*nbasis_bsmx ! ibs_cn2ucn 54 mapidata = mapidata + 3*nbasis_bsmx ! ibs_cn2ce 55 mapidata = mapidata + 3*nbasis_bsmx ! ibs_ce2uce 56 mapidata = mapidata + 3*nbasis_bsmx ! ibs_cn2bfr 57 mapidata = mapidata + 3*nbasis_bsmx ! ibs_ce2cnr 58 mapidata = mapidata + nbasis_bsmx ! ncont_tot_gb 59 mapidata = mapidata + nbasis_bsmx ! nprim_tot_gb 60 mapidata = mapidata + nbasis_bsmx ! nbf_tot_gb 61 mapidata = mapidata + nbasis_bsmx ! ibs_geom 62 mapidata = 4*mapidata 63c 64c.. total space 65 total4 = idata + mapidata 66 total8 = 2*total4 + rdata + cdata 67 total4 = total4 + rdata + cdata 68c 69 write(LuOut,'(////1x,a,f5.2,a)') 70 & ' **** basis set version ',bsversion,' ****' 71 write(LuOut,'(1x,a,i20,a)') 72 & ' character data in-core ',cdata,' bytes' 73 write(LuOut,'(1x,a,i20,a)') 74 & ' real data in-core ',rdata,' bytes' 75 write(LuOut,'(1x,a,i20,a)') 76 & ' integer*4 data in-core ',idata,' bytes' 77 write(LuOut,'(1x,a,i20,a)') 78 & 'or integer*8 data in-core ',(2*idata),' bytes' 79 write(LuOut,'(1x,a,i20,a)') 80 & ' integer*4 mapping data in-core ', 81 & mapidata,' bytes' 82 write(LuOut,'(1x,a,i20,a/)') 83 & 'or integer*8 mapping data in-core ', 84 & (2*mapidata),' bytes' 85 write(LuOut,*)' total(4) = ',total4,' bytes' 86 write(LuOut,*)' total(8) = ',total8,' bytes' 87c 88c.. convert to kilobytes 89c 90 cdata = (cdata + 999) / 1000 91 rdata = (rdata + 999) / 1000 92 idata = (idata + 999) / 1000 93 mapidata = (mapidata + 999) / 1000 94 total4 = (total4 + 999) / 1000 95 total8 = (total8 + 999) / 1000 96 write(LuOut,'(///1x,a,f5.2,a)') 97 & ' **** basis set version ',bsversion,' ****' 98 write(LuOut,'(1x,a,i20,a)') 99 & ' character data in-core ',cdata,' Kbytes' 100 write(LuOut,'(1x,a,i20,a)') 101 & ' real data in-core ',rdata,' Kbytes' 102 write(LuOut,'(1x,a,i20,a)') 103 & ' integer*4 data in-core ',idata,' Kbytes' 104 write(LuOut,'(1x,a,i20,a)') 105 & 'or integer*8 data in-core ',(2*idata),' Kbytes' 106 write(LuOut,'(1x,a,i20,a)') 107 & ' integer*4 mapping data in-core ', 108 & mapidata,' Kbytes' 109 write(LuOut,'(1x,a,i20,a/)') 110 & 'or integer*8 mapping data in-core ', 111 & (2*mapidata),' Kbytes' 112 write(LuOut,*)' total(4) = ',total4,' Kbytes' 113 write(LuOut,*)' total(8) = ',total8,' Kbytes' 114 write(LuOut,'(////)') 115c 116 bas_version = .true. 117 end 118*..................................................................... 119C> \ingroup bas 120C> @{ 121c 122C> \brief Creates a new basis instance 123c 124C> \return Returns .true. if successfull, and .false. otherwise 125c 126 logical function bas_create(basis,name) 127c 128c creates a handle and marks it active in the in-core data structure 129c 130 implicit none 131#include "mafdecls.fh" 132#include "nwc_const.fh" 133#include "basP.fh" 134#include "geobasmapP.fh" 135#include "inp.fh" 136#include "basdeclsP.fh" 137#include "bas_exndcf_dec.fh" 138#include "stdio.fh" 139c::functions 140c ifill from util 141c dfill from util 142c::passed 143c 144 integer basis !< [Output] the basis handle 145 character*(*) name !< [Input] the name of basis set 146c 147*:debug: integer iiii, jjjj 148 integer ii ! dummy loop counter 149c 150 external basis_data ! This for the T3D linker 151c 152#include "bas_exndcf_sfn.fh" 153c 154 do 00100 basis=1,nbasis_bsmx 155 if (.not.bsactive(basis)) goto 01000 15600100 continue 157c 158 write(LuOut,*)' bas_create: no free basis handles for ',name 159 bas_create = .false. 160 return 161c 16201000 continue 163c 164c store some information in basis data structure 165c (NOTE: name discarded in LOAD operation) 166c 167 bs_name(basis) = name 168 len_bs_name(basis) = inp_strlen(name) 169c 170c Initialize basis info to be empty 171c 172 bs_trans(basis) = ' ' 173 call ifill(ndbs_head, 0, infbs_head(1,basis), 1) 174 exndcf(H_exndcf ,basis) = -1 ! handle 175 exndcf(K_exndcf ,basis) = 0 ! index 176 exndcf(SZ_exndcf,basis) = 0 ! allocated size 177 call ifill(ndbs_tags*ntags_bsmx, 0, 178 $ infbs_tags(1,1,basis), 1) 179 call ifill(ndbs_ucont*nucont_bsmx, 0, 180 $ infbs_cont(1,1,basis), 1) 181 do ii = 1,ntags_bsmx 182 bs_stdname(ii,basis) = ' ' 183 bs_tags(ii,basis) = ' ' 184 enddo 185 do ii = 1,nbasis_assoc_max 186 name_assoc(ii,basis) = ' ' 187 enddo 188 call ifill(nbasis_assoc_max,0,handle_assoc(1,basis),1) 189 call ifill(nbasis_assoc_max,0,parent_assoc(1,basis),1) 190 bas_nassoc(basis) = 0 191c 192c Initialize geo-basis info to empty 193c 194 call ifill(3, 0, ibs_cn2ucn(1,basis), 1) 195 call ifill(3, 0, ibs_cn2ce (1,basis), 1) 196 call ifill(3, 0, ibs_ce2uce(1,basis), 1) 197 call ifill(3, 0, ibs_cn2bfr(1,basis), 1) 198 call ifill(3, 0, ibs_ce2cnr(1,basis), 1) 199 nbfmax_bs(basis) = -565 200 ncont_tot_gb(basis) = 0 201 nprim_tot_gb(basis) = 0 202 ncoef_tot_gb(basis) = 0 203 nbf_tot_gb(basis) = 0 204 ibs_geom(basis) = 0 205 bas_norm_id(basis) = BasNorm_UN 206c 207c Mark basis as active and return info 208c 209 bsactive(basis) = .true. 210 basis = basis - Basis_Handle_Offset 211 bas_create = .true. 212c 213c debug print all basis sets that are active 214c 215*:debug: write(LuOut,*)' bas_create: active basis sets' 216*:debug: do iiii = 1,nbasis_bsmx 217*:debug: if (bsactive(iiii)) then 218*:debug: jjjj = inp_strlen(bs_name(iiii)) 219*:debug: write(LuOut,*)'bas_create:',bs_name(iiii)(1:jjjj),' ',iiii 220*:debug: endif 221*:debug: enddo 222c 223 end 224*..................................................................... 225c 226C> \brief Destroys a basis instance 227c 228C> \return Returns .true. if successfull, and .false. otherwise 229c 230 logical function bas_destroy(basisin) 231 implicit none 232c::functions 233 logical bas_get_ecp_handle, ecp_check_handle, bas_do_destroy 234 external bas_get_ecp_handle, ecp_check_handle, bas_do_destroy 235 logical bas_get_so_handle, so_check_handle 236 external bas_get_so_handle, so_check_handle 237c::passed 238 integer basisin !< [Input] the basis handle 239c::local 240 logical ignore 241 integer ecpid, soid 242c 243 ignore = bas_get_ecp_handle(basisin,ecpid) 244 ignore = bas_get_so_handle(basisin,soid) 245 bas_destroy = .true. 246 if (ecp_check_handle(ecpid,'bas_destroy')) then 247 bas_destroy = bas_destroy.and.bas_do_destroy(ecpid) 248 endif 249 if (so_check_handle(soid,'bas_destroy')) then 250 bas_destroy = bas_destroy.and.bas_do_destroy(soid) 251 endif 252 bas_destroy = bas_destroy .and. bas_do_destroy(basisin) 253 end 254*..................................................................... 255C> @} 256 logical function bas_do_destroy(basisin) 257c 258c destroys information about an active incore basis 259c and the associated mapping arrays. 260c 261 implicit none 262#include "errquit.fh" 263#include "mafdecls.fh" 264#include "basdeclsP.fh" 265#include "nwc_const.fh" 266#include "basP.fh" 267#include "bas_exndcf_dec.fh" 268#include "ecpso_decP.fh" 269#include "geomP.fh" 270#include "geobasmapP.fh" 271#include "bas_ibs_dec.fh" 272#include "stdio.fh" 273c::functions used 274c::bas 275 logical bas_check_handle 276 logical gbs_map_clear 277 logical geom_ncent 278 logical ecp_get_num_elec 279 external gbs_map_clear 280 external bas_check_handle 281 external geom_ncent 282 external ecp_get_num_elec 283c::passed 284 integer basisin ![input] basis set handle to be destroyed 285c::local 286 integer idbstag 287 integer i, nat, num_elec 288 integer geom 289 integer basis 290 integer h_tmp 291c 292#include "bas_ibs_sfn.fh" 293#include "bas_exndcf_sfn.fh" 294#include "ecpso_sfnP.fh" 295c 296 bas_do_destroy = bas_check_handle(basisin,'bas_do_destroy') 297 if (.not. bas_do_destroy) return 298 299 basis = basisin + Basis_Handle_Offset 300 301c 302c must restore active geometry data appropriately 303 geom = ibs_geom(basis) 304 if (active(geom)) then 305 if (Is_ECP(basis)) then 306 if (.not.geom_ncent(geom,nat)) 307 & call errquit('bas_do_destroy:geom_ncent failed',911, 308 & BASIS_ERR) 309 do i = 1,nat 310 idbstag = sf_ibs_ce2uce(i,basis) 311 if (ecp_get_num_elec(basisin, 312 & bs_tags(idbstag,basis),num_elec)) then 313 charge(i,geom) = charge(i,geom) + dble(num_elec) 314 endif 315 enddo 316 endif 317 endif 318c 319 if(.not.gbs_map_clear(basisin)) then 320 write(LuOut,*)' error clearing map ' 321 bas_do_destroy = .false. 322 return 323 endif 324 325c 326 angular_bs(basis) = -565 327 bas_norm_id(basis) = -565 328 nbfmax_bs(basis) = -565 329 bsactive(basis) = .false. 330 bas_spherical(basis) = .false. 331 bas_any_gc(basis) = .false. 332 bas_any_sp_shell(basis) = .false. 333c 334 h_tmp = exndcf(H_exndcf,basis) 335 if (h_tmp .ne. -1) then 336 if (.not.ma_free_heap(h_tmp)) call errquit 337 & ('bas_do_destroy: error freeing heap',911, MEM_ERR) 338 endif 339 exndcf(H_exndcf ,basis) = -1 340 exndcf(K_exndcf ,basis) = 0 341 exndcf(SZ_exndcf,basis) = 0 342c 343 bas_do_destroy = .true. 344 end 345*..................................................................... 346C> \ingroup bas 347C> @{ 348c 349C> \brief Checks whether a given handle refers to a valid basis instance 350c 351C> \return Return .true. if basisin is a valid basis instance, and 352C> .false. otherwise 353c 354 logical function bas_check_handle(basisin,msg) 355c 356c Checks to see if a basis set handle is valid 357c 358 implicit none 359#include "nwc_const.fh" 360#include "basP.fh" 361#include "stdio.fh" 362c:passed 363 integer basisin !< [Input] the basis handle 364 character*(*) msg !< [Input] an error message 365c::local 366 integer basis 367c 368 basis = basisin + Basis_Handle_Offset 369 bas_check_handle = basis.gt.0 .and. basis.le.nbasis_bsmx 370 if (bas_check_handle) 371 & bas_check_handle = bas_check_handle .and. bsactive(basis) 372c 373* user's responsibility to deal with status 374* if (.not. bas_check_handle) then 375* write(LuOut,*)msg,': basis handle is invalid ' 376* write(LuOut,*)'basis_check_handle: lexical handle ',basis 377* write(LuOut,*)'basis_check_handle: handle ',basisin 378* endif 379 return 380 end 381*..................................................................... 382c 383C> \brief Checks whether a given handle refers to a valid ECP basis 384C> instance 385c 386C> \return Return .true. if ecpidin is a valid ECP basis instance, and 387C> .false. otherwise 388c 389 logical function ecp_check_handle(ecpidin,msg) 390c 391c Checks to see if an ECP basis set handle is valid 392c 393 implicit none 394#include "basdeclsP.fh" 395#include "nwc_const.fh" 396#include "basP.fh" 397#include "ecpso_decP.fh" 398#include "inp.fh" 399c:functions 400 logical bas_check_handle 401 external bas_check_handle 402c:passed 403 integer ecpidin !< [Input] the ECP basis handle 404 character*(*) msg !< [Input] an error message 405c::local 406 character*255 newmsg 407c 408#include "ecpso_sfnP.fh" 409c 410 newmsg(1:17) = 'ecp_check_handle:' 411 newmsg(18:) = msg 412 ecp_check_handle = bas_check_handle(ecpidin,newmsg) 413 if(.not.ecp_check_handle) return 414 ecp_check_handle = ecp_check_handle .and. 415 & Is_ECP_in(ecpidin) 416 return 417 end 418*..................................................................... 419c 420C> \brief Checks whether a given handle refers to a valid spin-orbit 421C> potential instance 422c 423C> \return Return .true. if soidin is a valid spin-orbit potential 424C> instance, and .false. otherwise 425c 426 logical function so_check_handle(soidin,msg) 427c 428c Checks to see if an ECP so basis set handle is valid 429c 430 implicit none 431#include "basdeclsP.fh" 432#include "nwc_const.fh" 433#include "basP.fh" 434#include "ecpso_decP.fh" 435#include "inp.fh" 436c:functions 437 logical bas_check_handle 438 external bas_check_handle 439c:passed 440 integer soidin !< [Input] the spin-orbit ECP basis handle 441 character*(*) msg !< [Input] an error message 442c::local 443 character*255 newmsg 444c 445#include "ecpso_sfnP.fh" 446c 447 newmsg(1:17) = 'so_check_handle:' 448 newmsg(18:) = msg 449 so_check_handle = bas_check_handle(soidin,newmsg) 450 if(.not.so_check_handle) return 451 so_check_handle = so_check_handle .and. 452 & Is_SO_in(soidin) 453 return 454 end 455*..................................................................... 456c 457C> \brief Print a summary of a basis instance 458c 459C> \return Return .true. if successfull, and .false. otherwise 460c 461 logical function bas_summary_print(basisin) 462 implicit none 463#include "errquit.fh" 464#include "stdio.fh" 465#include "mafdecls.fh" 466#include "nwc_const.fh" 467#include "basP.fh" 468#include "basdeclsP.fh" 469#include "inp.fh" 470#include "ecpso_decP.fh" 471#include "bas_starP.fh" 472c::-functions 473 logical bas_check_handle 474 external bas_check_handle 475 integer nbf_from_ucont 476 external nbf_from_ucont 477c::-passed 478 integer basisin !< [Input] the basis set handle 479c::-local 480 character*16 dum_tag 481 character*255 dum_stdtag 482 character*12 polynomial ! string for spherical/cartesian 483 character*80 dum_type, dum_string 484 integer basis ! lexical index 485 integer i_tag ! loop counter 486 integer i_cont ! loop counter 487 integer mytags ! dummy 488 integer myfcont ! loop range value (first) 489 integer mylcont ! loop range value (last) 490 integer mycont ! number of contractions 491 integer mynbf ! number of functions 492 integer tmp1, tmp2, tmp3 493 integer myngen 494 integer mytype 495 integer jtype 496 character*1 ctype(0:6) 497 integer cnt_type(0:6) 498* 499#include "ecpso_sfnP.fh" 500* 501 bas_summary_print = 502 & bas_check_handle(basisin,'bas_summary_print') 503* 504 if (Is_ECP_in(basisin).or.Is_SO_in(basisin)) then 505 bas_summary_print = .true. 506 return 507 endif 508* 509 ctype(0)='s' 510 ctype(1)='p' 511 ctype(2)='d' 512 ctype(3)='f' 513 ctype(4)='g' 514 ctype(5)='h' 515 ctype(6)='i' 516* 517 basis = basisin + Basis_Handle_Offset 518* 519 if (bas_spherical(basis)) then 520 polynomial = ' (spherical)' 521 else 522 polynomial = ' (cartesian)' 523 endif 524 write(luout,'(/)') 525 write(luout,10000) 526 & bs_name(basis)(1:inp_strlen(bs_name(basis))), 527 & bs_trans(basis)(1:inp_strlen(bs_trans(basis))), 528 & polynomial 529 mytags = infbs_head(Head_Ntags,basis) 530 do i_tag = 1,mytags 531 call ifill(7,0,cnt_type,1) 532 myfcont = infbs_tags(Tag_Fcont,i_tag,basis) 533 mylcont = infbs_tags(Tag_Lcont,i_tag,basis) 534 mycont = mylcont - myfcont + 1 535 mynbf = 0 536 do i_cont = myfcont, mylcont 537 mynbf = mynbf + nbf_from_ucont(i_cont,basisin) 538 mytype = infbs_cont(Cont_Type,i_cont,basis) 539 myngen = infbs_cont(Cont_Ngen,i_cont,basis) 540 if (myngen.lt.1) call errquit ( 541 & 'bas_summary_print: fatal error myngen:',myngen, 542 & BASIS_ERR) 543 if (mytype.ge.0) then 544 cnt_type(mytype) = cnt_type(mytype) + myngen 545 else 546 do jtype = 0,abs(mytype) 547 cnt_type(jtype) = cnt_type(jtype) + 1 548 enddo 549 endif 550 enddo 551 dum_tag = bs_tags(i_tag,basis) 552 tmp1 = inp_strlen(bs_stdname(i_tag,basis)) 553 if (tmp1 .lt. (30-1)) then 554 tmp2 = (30-tmp1)/2 555 else 556 tmp2 = 1 557 endif 558 dum_stdtag = ' ' 559 dum_stdtag(tmp2:) = bs_stdname(i_tag,basis) 560 dum_type = ' ' 561 dum_string = ' ' 562 do jtype = 0,6 563 dum_type = dum_string 564 tmp1 = inp_strlen(dum_type) 565 if (tmp1.lt.1) tmp1 = 1 566 if (cnt_type(jtype).gt.0) then 567 write(dum_string,'(a,i5,a)')dum_type(1:tmp1), 568 & cnt_type(jtype),ctype(jtype) 569 endif 570 enddo 571 tmp1 = 0 572 dum_type = dum_string 573 do jtype = 1,inp_strlen(dum_type) 574 if (dum_type(jtype:jtype).ne.char(0).and. 575 & dum_type(jtype:jtype).ne.' ') then 576 tmp1 = tmp1 + 1 577 dum_string(tmp1:tmp1) = dum_type(jtype:jtype) 578 endif 579 enddo 580 tmp1 = tmp1 + 1 581 do jtype = tmp1,len(dum_string) 582 dum_string(jtype:jtype) = ' ' 583 enddo 584 tmp1 = inp_strlen(dum_string) 585 write(luout,10001) 586 & dum_tag, 587 & dum_stdtag(1:30), 588 & mycont, 589 & mynbf,dum_string(1:tmp1) 590 enddo 591 do i_tag=1,star_nr_tags 592 tmp1 = inp_strlen(star_bas_typ(i_tag)) 593 if (tmp1 .lt. (30-1)) then 594 tmp2 = (30-tmp1)/2 595 else 596 tmp2 = 1 597 endif 598 dum_stdtag = ' ' 599 dum_stdtag(tmp2:) = star_bas_typ(i_tag) 600 tmp1 = 1 601 if (i_tag .gt. 1) tmp1 = star_nr_excpt(i_tag-1) + 1 602 if ((star_nr_excpt(i_tag) - tmp1) .gt. -1) then 603 write(LuOut,10002) star_tag(i_tag), 604 & dum_stdtag(1:30), 'all atoms except', 605 & (star_excpt(tmp3)(1:inp_strlen(star_excpt(tmp3))), 606 & tmp3=tmp1,star_nr_excpt(i_tag)) 607 else 608 write(LuOut,10002) star_tag(i_tag), 609 & dum_stdtag(1:30), 'on all atoms ' 610 endif 611 enddo 612 613 write(luout,'(/)') 614 call util_flush(luout) 615 bas_summary_print = .true. 616* 61710000 format(' Summary of "',a,'" -> "',a,'"',a/ 618 $ ' ---------------------------------------------', 619 & '---------------------------------'/ 620 $ ' Tag Description ', 621 & 'Shells Functions and Types'/ 622 $ ' ---------------- ------------------------------ ', 623 & '------ ---------------------') 62410001 format(1x,a16,1x,a30,2x,i4,5x,i4,3x,a) 62510002 format(1x,a16,1x,a30,2x,a17,1x,10(a)) 626* 627 end 628*..................................................................... 629c 630C> \brief Print the contents of a basis instance 631c 632C> \return Return .true. if successfull, and .false. otherwise 633c 634 logical function bas_print(basisin) 635c 636c routine to print unique basis information that is in core 637c 638 implicit none 639#include "mafdecls.fh" 640#include "nwc_const.fh" 641#include "basP.fh" 642#include "basdeclsP.fh" 643#include "inp.fh" 644#include "geom.fh" 645#include "ecpso_decP.fh" 646#include "stdio.fh" 647#include "bas_starP.fh" 648c 649c function declarations 650c 651 logical bas_check_handle, ecp_print, so_print 652 external bas_check_handle, ecp_print, so_print 653c:: passed 654 integer basisin !< [Input] the basis set handle 655c:: local 656 integer mytags, myucont, myprim, mycoef, basis 657 integer i,j,k,l, ifcont, mygen, mytype, iexptr, icfptr 658 integer atn, len_tag, len_ele 659 character*2 symbol 660 character*16 element 661 character*3 ctype(0:6),cltype(2) 662 character*3 shell_type 663*. . . . . . . . . . . ! Room for tag+space+(+element+) = 16+1+1+16+1 664 character*35 buffer 665 character*12 polynomial 666c 667#include "bas_exndcf.fh" 668#include "ecpso_sfnP.fh" 669c 670 if (Is_ECP_in(basisin)) then 671 bas_print = ecp_print(basisin) 672 return 673 endif 674 if (Is_SO_in(basisin)) then 675 bas_print = so_print(basisin) 676 return 677 endif 678c 679 ctype(0)='S' 680 ctype(1)='P' 681 ctype(2)='D' 682 ctype(3)='F' 683 ctype(4)='G' 684 ctype(5)='H' 685 ctype(6)='I' 686 cltype(1)='SP' 687 cltype(2)='SPD' 688 bas_print = .true. 689 basis = basisin + Basis_Handle_Offset 690c 691 bas_print = bas_check_handle(basisin,'bas_print') 692 if (.not. bas_print) return 693c 694c print basis set information 695c 696 if (bas_spherical(basis)) then 697 polynomial = ' (spherical)' 698 else 699 polynomial = ' (cartesian)' 700 endif 701 write(LuOut,1)bs_name(basis)(1:inp_strlen(bs_name(basis))), 702 $ bs_trans(basis)(1:inp_strlen(bs_trans(basis))), polynomial 703 1 format(' Basis "',a,'" -> "',a,'"',a/ 704 $ ' -----') 705 mytags = infbs_head(HEAD_NTAGS,basis) 706 if (mytags.le.0) then 707 write(LuOut,*)'No explicit basis set is defined !' 708 write(LuOut,*) 709c 710c there could be star tags defined, so check that before returning 711c 712 goto 00010 713 endif 714c 715 myucont = infbs_head(HEAD_NCONT,basis) 716 myprim = infbs_head(HEAD_NPRIM,basis) 717 mycoef = infbs_head(HEAD_NCOEF,basis) 718c 719 do 00100 i=1,mytags 720 721 if (geom_tag_to_element(bs_tags(i,basis), symbol, element, 722 $ atn)) then 723 len_tag = inp_strlen(bs_tags(i,basis)) 724 len_ele = inp_strlen(element) 725 write(buffer,'(a,'' ('',a,'')'')') 726 $ bs_tags(i,basis)(1:len_tag), element(1:len_ele) 727 else 728 buffer = bs_tags(i,basis) 729 endif 730 len_tag = inp_strlen(buffer) 731 call util_print_centered(LuOut, buffer, len_tag/2 + 1, .true.) 732 733 myucont = infbs_tags(TAG_NCONT,i,basis) 734c 735 ifcont = infbs_tags(TAG_FCONT,i,basis) 736c 737 write(LuOut,6) 738 6 format( 739 $ ' Exponent Coefficients '/ 740 $ ' -------------- ',57('-')) 741 do 00200 j=1,myucont 742 myprim = infbs_cont(CONT_NPRIM,ifcont,basis) 743 mygen = infbs_cont(CONT_NGEN,ifcont,basis) 744 745 mytype = infbs_cont(CONT_TYPE, ifcont, basis) 746 if (mytype.lt.0) then 747 shell_type = cltype(abs(mytype)) 748 else 749 shell_type = ctype(mytype) 750 endif 751 iexptr = infbs_cont(CONT_IEXP,ifcont,basis) - 1 752 icfptr = infbs_cont(CONT_ICFP,ifcont,basis) - 1 753 do 00300 k=1,myprim 754 write(LuOut,7) j, shell_type(1:2), 755 & sf_exndcf((iexptr+k),basis), 756 & (sf_exndcf((icfptr+k+(l-1)*myprim),basis),l=1,mygen) 757 7 format(1x,i2,1x,a2,1x,1pE14.8,0p20f10.6) 75800300 continue 759 write(LuOut,*) 760 ifcont = ifcont + 1 76100200 continue 76200100 continue 763c 764c Check if we have star tag definitions in the basis set 765c 76600010 if (star_nr_tags .gt. 0) then 767 write(LuOut,*) 768 write(LuOut,8) 769 8 format(' In addition, one or more string tags have been', 770 & ' defined containing a * .'/' These tags, and ', 771 & 'their exceptions list are printed below.'// 772 & ' Tag ',12(' '),' Description ',18(' '), 773 & ' Excluding '/' ',16('-'),' ',30('-'),' ',30('-')) 774 do i=1,star_nr_tags 775 k = 1 776 if (i .gt. 1) k = star_nr_excpt(i-1) + 1 777 write(LuOut,9) star_tag(i), star_bas_typ(i)(1:30), 778 & (star_excpt(j)(1:inp_strlen(star_excpt(j))), 779 & j=k,star_nr_excpt(i)) 780 9 format(1x,a16,1x,a30,1x,10(a)) 781 enddo 782 write(LuOut,*) 783 write(LuOut,*) 784 endif 785c 786c If geom is set print out the info about total basis info 787c associated with the geometry also 788c 789c ... not done yet 790c 791 return 792 end 793*..................................................................... 794c 795C> \brief Load a basis set from the RTDB 796c 797C> Retrieve the basis set from the RTDB using the name specified. 798C> The data is stored in the basis set instance specified. Furthermore 799C> some tables are constructed linking basis functions to atoms, etc. 800C> For this purpose a geometry needs to be provided as well. 801c 802C> \return Return .true. if the basis was loaded successfully, and 803C> .false. otherwise 804c 805 logical function bas_rtdb_load(rtdb, geom, basisin, name) 806 implicit none 807#include "errquit.fh" 808#include "nwc_const.fh" 809#include "global.fh" 810#include "rtdb.fh" 811#include "mafdecls.fh" 812#include "context.fh" 813#include "stdio.fh" 814#include "inp.fh" 815#include "basP.fh" 816#include "basdeclsP.fh" 817c::functions 818 logical bas_rtdb_do_load, bas_get_ecp_name, bas_set_ecp_name 819 logical bas_create, bas_do_destroy, bas_set_ecp_handle 820 logical ecp_set_parent_handle, bas_summary_print 821 logical bas_name_exist_rtdb, bas_print 822 external bas_rtdb_do_load, bas_get_ecp_name, bas_set_ecp_name 823 external bas_create, bas_do_destroy, bas_set_ecp_handle 824 external ecp_set_parent_handle, bas_summary_print 825 external bas_name_exist_rtdb, bas_print 826 logical bas_get_so_name, bas_set_so_name, bas_set_so_handle 827 logical so_set_parent_handle 828 external bas_get_so_name, bas_set_so_name, bas_set_so_handle 829 external so_set_parent_handle 830c::passed 831 integer rtdb !< [Input] the RTDB handle 832 integer geom !< [Input] the geometry handle 833 integer basisin !< [Input] the basis set handle 834 character*(*) name !< [Input] the name of the basis set on the 835c !< RTDB 836c 837c::local 838 character*256 ecp_name, so_name 839 integer ecpid, soid 840 logical status 841 logical status_ecp_cr, status_so_cr 842 logical status_ecp_load, status_so_load 843 logical status_ecp_ph, status_so_ph 844 logical status_ecp, status_so 845c 846*:debug: write(LuOut,*)' bas_rtdb_load:rtdb,geom,basis ', 847*:debug: & rtdb,geom,basisin 848*:debug: write(LuOut,*)' bas_rtdb_load:name ',name 849 bas_rtdb_load = 850 & bas_rtdb_do_load(rtdb, geom, basisin, name) 851 if (.not.bas_rtdb_load) return 852 if (.not.inp_compare(.false.,'ao basis',name(1:8))) 853 & goto 112 854* 855 if (.not.bas_get_ecp_name(basisin,ecp_name)) call errquit 856 & ('bas_rtdb_load: bas_get_ecp_name failed',911, RTDB_ERR) 857 if (ecp_name .eq. ' ') then 858 ecp_name = 'ecp basis' 859 endif 860 861 status = bas_name_exist_rtdb(rtdb,ecp_name) 862 863 if (status) then 864* 865* ecp_name is on rtdb so load it 866* 867 status_ecp_cr = bas_create(ecpid,ecp_name) 868 if (.not.status_ecp_cr) then 869 call errquit 870 & ('bas_rtdb_load: bas_create failed for ecpid',911, 871 & RTDB_ERR) 872 endif 873 874 status_ecp_load = bas_rtdb_do_load(rtdb,geom,ecpid,ecp_name) 875 if (status_ecp_load) then 876*.... everything is okay 877 if (.not.bas_set_ecp_name(basisin,ecp_name)) call errquit 878 & ('bas_rtdb_load: bas_set_ecp_name failed',911, 879 & RTDB_ERR) 880 881 status_ecp = bas_set_ecp_handle(basisin,ecpid) 882 if (.not.status_ecp) then 883 call errquit 884 & ('bas_rtdb_load: bas_set_ecp_handle failed',911, 885 & RTDB_ERR) 886 endif 887 888 status_ecp_ph = ecp_set_parent_handle(ecpid,basisin) 889 if (.not.status_ecp_ph) then 890 call errquit 891 & ('bas_rtdb_load: ecp_set_parent_handle failed',911, 892 & RTDB_ERR) 893 endif 894 else 895*... ecp exists but is not used by current geometry so destroy it ! 896 if (.not.bas_do_destroy(ecpid)) then 897 write(luout,*)' unused ecp basis failed to be destroyed' 898 call errquit('bas_rtdb_load: bas_do_destroy failed',911, 899 & RTDB_ERR) 900 endif 901 endif 902 endif 903 904 if (.not.bas_get_so_name(basisin,so_name)) call errquit 905 & ('bas_rtdb_load: bas_get_so_name failed',911, RTDB_ERR) 906 if (so_name .eq. ' ') then 907 so_name = 'so potential' 908 endif 909 910 status = bas_name_exist_rtdb(rtdb,so_name) 911 912 if (status) then 913* 914* so_name is on rtdb so load it 915* 916 status_so_cr = bas_create(soid,so_name) 917 if (.not.status_so_cr) then 918 call errquit 919 & ('bas_rtdb_load: bas_create failed for soid',911, 920 & RTDB_ERR) 921 endif 922 923 status_so_load = bas_rtdb_do_load(rtdb,geom,soid,so_name) 924 if (status_so_load) then 925*.... everything is okay 926 if (.not.bas_set_so_name(basisin,so_name)) call errquit 927 & ('bas_rtdb_load: bas_set_so_name failed',911, 928 & RTDB_ERR) 929 930 status_so = bas_set_so_handle(basisin,soid) 931 if (.not.status_so) then 932 call errquit 933 & ('bas_rtdb_load: bas_set_so_handle failed',911, 934 & RTDB_ERR) 935 endif 936 937 status_so_ph = so_set_parent_handle(soid,basisin) 938 if (.not.status_so_ph) then 939 call errquit 940 & ('bas_rtdb_load: so_set_parent_handle failed',911, 941 & RTDB_ERR) 942 endif 943 else 944*... so exists but is not used by current geometry so destroy it ! 945 if (.not.bas_do_destroy(soid)) then 946 write(luout,*)' unused so basis failed to be destroyed' 947 call errquit('bas_rtdb_load: bas_do_destroy failed',911, 948 & RTDB_ERR) 949 endif 950 endif 951 endif 952c 953c Check if we need to print the basis set. The basis set might not 954c have been printed, this happens when the user asks for the basis 955c to be printed but the basis set contains star tags, which only can 956c can be resolved at the task level. 957c 958c Reusing ecp_name and so_name and status 959c 960 112 if (.not.context_rtdb_match(rtdb,name,ecp_name)) 961 $ ecp_name = name 962 so_name = 'basisprint:'//ecp_name 963 $ (1:inp_strlen(ecp_name)) 964 if (rtdb_get(rtdb,so_name,mt_log,1,status)) then 965 if (status) then 966 if (ga_nodeid() .eq. 0) then 967 if (.not. bas_print(basisin)) 968 $ call errquit('bas_rtdb_load: print failed', 0, 969 & RTDB_ERR) 970 if (.not.bas_summary_print(basisin)) call errquit 971 $ ('bas_rtdb_load: basis summary print failed',911, 972 & RTDB_ERR) 973 endif 974 status = .false. 975 if (.not. rtdb_put(rtdb,so_name,mt_log,1,status)) 976 $ call errquit('bas_rtdb_load: rtdb_put failed',911, 977 & RTDB_ERR) 978 endif 979 endif 980 call bas_ecce_print_basis(basisin,'bas_input') 981c 982 end 983*..................................................................... 984C> @} 985 logical function bas_rtdb_do_load(rtdb, geom, basisin, name) 986 implicit none 987#include "errquit.fh" 988c 989c routine that loads a basis set from the rtdb and using the 990c geometry information builds the mapping arrays to contractions/ 991c shells, basis functions, and centers. 992c 993#include "rtdb.fh" 994#include "mafdecls.fh" 995#include "nwc_const.fh" 996#include "geomP.fh" 997#include "geom.fh" 998#include "context.fh" 999#include "basP.fh" 1000#include "basdeclsP.fh" 1001#include "geobasmapP.fh" 1002#include "inp.fh" 1003#include "global.fh" 1004#include "stdio.fh" 1005#include "bas_starP.fh" 1006#include "bas.fh" 1007c 1008c function declarations 1009c 1010c::function 1011c:bas 1012 logical bas_geobas_build 1013 logical bas_add_ucnt_init 1014 logical bas_add_ucnt_tidy 1015 external bas_geobas_build 1016 external bas_add_ucnt_init 1017 external bas_add_ucnt_tidy 1018c:: passed 1019 integer rtdb ! [input] rtdb handle 1020 integer geom ! [input] geometry handle with info loaded 1021 integer basisin ! [input] basis handle 1022 character*(*) name ! [input] basis set name that must be on 1023*. . . . . . . . . . . . . the rtdb 1024c:: local 1025 integer ecXtra ! amount of extra space required for zero shell info 1026 parameter (ecXtra = 2) 1027 character*25 nameex 1028 integer lentmp, basis, nexcf 1029 character*256 tmp 1030 logical rtdb_status 1031 integer h_tmp, k_tmp, h_new, k_new 1032 integer iunique, uniquecent, istar , ntag_read 1033 integer istart, iexcpt, ilen, junique 1034 integer i_array,l_array 1035 logical oIsexcept 1036 character*16 tag_to_add, tag_in_lib 1037c:: statement functions 1038#include "bas_exndcf.fh" 1039c:: initalize local 1040c 1041 rtdb_status = .true. 1042c 1043c 1044c check geom and basis handles returns false if either is invalid 1045c 1046 bas_rtdb_do_load = geom_check_handle(geom,'bas_rtdb_do_load') 1047 if (.not.bas_rtdb_do_load) return 1048 bas_rtdb_do_load = bas_check_handle(basisin,'bas_rtdb_do_load') 1049 if (.not.bas_rtdb_do_load) return 1050c 1051 basis = basisin + Basis_Handle_Offset 1052c 1053c store geom tag with basis map info 1054c 1055 ibs_geom(basis) = geom 1056c 1057c translate "name" to current "context" 1058c 1059 bs_name(basis) = name 1060 len_bs_name(basis) = inp_strlen(name) 1061 if (.not.context_rtdb_match(rtdb,name,bs_trans(basis))) 1062 & bs_trans(basis) = name 1063 len_bs_trans(basis) = inp_strlen(bs_trans(basis)) 1064c 1065c generate rtdb names and load information 1066c 1067 tmp = 'basis:'//bs_trans(basis)(1:len_bs_trans(basis)) 1068 lentmp = inp_strlen(tmp) + 1 1069c 1070 ntag_read=0 1071 tmp(lentmp:) = ' ' 1072 tmp(lentmp:) = ':bs_nr_tags' 1073 rtdb_status = rtdb_status .and. 1074 & rtdb_get(rtdb,tmp,mt_int,1,ntag_read) 1075 if (ntag_read .gt. 0) then 1076 tmp(lentmp:) = ' ' 1077 tmp(lentmp:) = ':bs_tags' 1078 rtdb_status = rtdb_status .and. 1079 & rtdb_cget(rtdb, tmp, ntags_bsmx, bs_tags(1,basis)) 1080c 1081 tmp(lentmp:) = ' ' 1082 tmp(lentmp:) = ':bs_stdname' 1083 rtdb_status = rtdb_status .and. 1084 & rtdb_cget(rtdb, tmp, ntags_bsmx, bs_stdname(1,basis)) 1085 endif 1086c 1087 tmp(lentmp:) = ' ' 1088 tmp(lentmp:) = ':assoc ecp name' 1089 rtdb_status = rtdb_status .and. 1090 & rtdb_cget(rtdb, tmp, 2, name_assoc(1,basis)) 1091c 1092 nexcf = 0 1093 tmp(lentmp:) = ' ' 1094 tmp(lentmp:) = ':number of exps and coeffs' 1095 rtdb_status = rtdb_status .and. 1096 & rtdb_get(rtdb,tmp,mt_int,1,nexcf) 1097c 1098 write(nameex,'(a23,i2)')' basis exps and coeffs ',basis 1099c 1100 if (exndcf(H_exndcf,basis) .ne. -1) then 1101 if (.not. ma_free_heap(exndcf(H_exndcf,basis))) 1102 $ call errquit('bas_rtdb_do_load: ma is corrupted', 1103 $ exndcf(H_exndcf,basis), RTDB_ERR) 1104 endif 1105 1106 if (.not.ma_alloc_get(mt_dbl,(nexcf+ecXtra),nameex, 1107 & h_tmp, k_tmp)) then 1108 write(LuOut,*)' not enough memory' 1109 call errquit 1110 & (' bas_rtdb_do_load: error allocating space'// 1111 & ' for exndcf',911, RTDB_ERR) 1112 else 1113 call dfill((nexcf+ecXtra),0.0d00,dbl_mb(k_tmp),1) 1114 exndcf(H_exndcf,basis) = h_tmp 1115 exndcf(K_exndcf,basis) = k_tmp 1116 exndcf(SZ_exndcf,basis) = (nexcf+ecXtra) 1117 endif 1118c 1119 if (nexcf .gt. 0) then 1120 tmp(lentmp:) = ' ' 1121 tmp(lentmp:) = ':exps and coeffs' 1122 rtdb_status = rtdb_status .and. 1123 & rtdb_get( 1124 & rtdb, tmp, mt_dbl, nexcf, dbl_mb(mb_exndcf(1,basis))) 1125 endif 1126c 1127 tmp(lentmp:) = ' ' 1128 tmp(lentmp:) = ':header' 1129 rtdb_status = rtdb_status .and. 1130 & rtdb_get( 1131 & rtdb, tmp, mt_int, ndbs_head, infbs_head(1,basis)) 1132c 1133 tmp(lentmp:) = ' ' 1134 tmp(lentmp:) = ':tags info' 1135 rtdb_status = rtdb_status .and. 1136 & rtdb_get( 1137 & rtdb, tmp, mt_int, 1138 & ndbs_tags*ntags_bsmx, infbs_tags(1,1,basis)) 1139c 1140 tmp(lentmp:) = ' ' 1141 tmp(lentmp:) = ':contraction info' 1142 rtdb_status = rtdb_status .and. 1143 & rtdb_get( 1144 & rtdb, tmp, mt_int, 1145 & ndbs_ucont*nucont_bsmx, infbs_cont(1,1,basis)) 1146c 1147c now deal with the star tags from the input. We have the geometry, so 1148c we extract the tags from the geometry, load the star tag info, and 1149c add the basis sets 1150c 1151c first check if we have any star tags to deal with at all 1152c account for old rtdbs without star tag info 1153c 1154 tmp(lentmp:) = ' ' 1155 tmp(lentmp:) = ':star nr tags' 1156 star_nr_tags = 0 1157 if (rtdb_status) rtdb_status = rtdb_status .and. rtdb_get(rtdb, 1158 & tmp,mt_int,1,star_nr_tags) 1159c 1160c read remaining star tag info and analyze only if we have actually a star 1161c tag input (hence, if star_nr_tags > 0) 1162c 1163 if (star_nr_tags .gt. 0) then 1164 tmp(lentmp:) = ' ' 1165 tmp(lentmp:) = ':star tag names' 1166 rtdb_status = rtdb_status .and. rtdb_cget( 1167 & rtdb,tmp,star_nr_tags,star_tag) 1168 tmp(lentmp:) = ' ' 1169 tmp(lentmp:) = ':star tag_in_lib' 1170 rtdb_status = rtdb_status .and. rtdb_cget( 1171 & rtdb,tmp,star_nr_tags,star_in_lib) 1172 tmp(lentmp:) = ' ' 1173 tmp(lentmp:) = ':star bas type' 1174 rtdb_status = rtdb_status .and. rtdb_cget( 1175 & rtdb,tmp,star_nr_tags,star_bas_typ) 1176 tmp(lentmp:) = ' ' 1177 tmp(lentmp:) = ':star filename' 1178 rtdb_status = rtdb_status .and. rtdb_cget( 1179 & rtdb,tmp,star_nr_tags,star_file) 1180 tmp(lentmp:) = ' ' 1181 tmp(lentmp:) = ':star tot except' 1182 rtdb_status = rtdb_status .and. rtdb_get( 1183 & rtdb,tmp,mt_int,1,star_tot_excpt) 1184 if (star_tot_excpt .gt. 0) then 1185 tmp(lentmp:) = ' ' 1186 tmp(lentmp:) = ':star except' 1187 rtdb_status = rtdb_status .and. rtdb_cget( 1188 & rtdb,tmp,star_tot_excpt,star_excpt) 1189 endif 1190 tmp(lentmp:) = ' ' 1191 tmp(lentmp:) = ':star nr except' 1192 rtdb_status = rtdb_status .and. rtdb_get( 1193 & rtdb,tmp,mt_int,star_nr_tags,star_nr_excpt) 1194 tmp(lentmp:) = ' ' 1195 tmp(lentmp:) = ':star rel' 1196 rtdb_status = rtdb_status .and. rtdb_get( 1197 & rtdb,tmp,mt_log,star_nr_tags,star_rel) 1198 tmp(lentmp:) = ' ' 1199 tmp(lentmp:) = ':star segment' 1200 rtdb_status = rtdb_status .and. rtdb_get( 1201 & rtdb,tmp,mt_log,1,star_segment) 1202c 1203c get a list of tags from the current geometry 1204c 1205 if (.not. geom_ncent_unique(geom,uniquecent)) 1206 & call errquit('bas_rtdb_do_load: geom_ncent_unique',211, 1207 & RTDB_ERR) 1208 if (.not.MA_alloc_Get(mt_int,uniquecent, 'iarray', 1209 & l_array, i_array)) 1210 & call errquit('basrtdbdoload: cannot allocate ',0, MA_ERR) 1211 if (.not. geom_uniquecent_get(geom,uniquecent, 1212 , int_mb(i_array))) 1213 & call errquit('bas_rtdb_do_load: geom_uniquecent_get',211, 1214 & RTDB_ERR) 1215c 1216c We got the tag_to_add info from the geometry object 1217c 1218 do iunique = 1, uniquecent 1219 if (.not. geom_cent_tag(geom, 1220 , int_mb(i_array+iunique-1),tag_to_add)) 1221 & call errquit('bas_rtdb_do_load: geom_cen_tag',211, 1222 & RTDB_ERR) 1223c 1224c First check if we did this geometry tag already 1225c 1226 do junique = 1, iunique - 1 1227 if (.not. geom_cent_tag(geom, 1228 , int_mb(i_array+junique-1),tag_in_lib)) 1229 & call errquit('bas_rtdb_do_load: geom_cen_tag',211, 1230 & RTDB_ERR) 1231 if (inp_compare(.true.,tag_to_add,tag_in_lib) 1232 & .and. (inp_strlen(tag_to_add) .eq. 1233 & inp_strlen(tag_in_lib))) goto 00012 1234 enddo 1235c 1236c 1237c Now for each star tag, check if tag_to_add matches and add basis set 1238c 1239 do istar = 1, star_nr_tags 1240c 1241c For * input, directly check exceptions list 1242c For aa*, first check if tag_to_add matches star tag itself 1243c For Bq, if B* is used and geometry tag is Bq, skip, only when Bq* 1244c is used we need to go on 1245c Skip and do not load a basis for dummy center X 1246c 1247 if (inp_contains(.true.,'*',star_tag(istar),ilen)) then 1248 if (ilen .gt. 1) then 1249 if (tag_to_add(1:ilen-1) .ne. 1250 & star_tag(istar)(1:ilen-1)) goto 00011 1251 endif 1252 if (inp_compare(.false.,'Bq',tag_to_add(1:2)).and..not. 1253 & inp_compare(.false.,'Bq',star_tag(istar)(1:2))) 1254 & goto 00011 1255 if (inp_compare(.false.,'X',tag_to_add(1:1)).and..not. 1256 & inp_compare(.false.,'Xe',tag_to_add(1:2))) 1257 & goto 00011 1258 endif 1259c 1260c There is a match between the tag_to_add and the star_tag, check for 1261c matches in exceptions list 1262c 1263 oIsexcept = .false. 1264 istart = 1 1265 if (istar .gt. 1) istart = star_nr_excpt(istar-1) + 1 1266 do iexcpt = istart, star_nr_excpt(istar) 1267 ilen = inp_strlen(tag_to_add) 1268 if (inp_compare(.true.,tag_to_add,star_excpt(iexcpt)) 1269 & .and. (inp_strlen(tag_to_add) .eq. 1270 & inp_strlen(star_excpt(iexcpt)))) 1271 & oIsexcept = .true. 1272 enddo 1273c 1274c If not in exceptions list, add basis set to the basis object 1275c 1276 if (.not. oIsexcept) then 1277 tag_in_lib = star_in_lib(istar) 1278 if (inp_contains(.true.,'*',star_in_lib(istar),ilen)) 1279 & tag_in_lib = tag_to_add 1280 call bas_tag_lib(basisin,star_segment,tag_to_add, 1281 & tag_in_lib, star_bas_typ(istar), star_file(istar), 1282 & star_rel(istar)) 1283 endif 128400011 continue 1285 enddo 128600012 continue 1287 enddo 1288 if (.not.MA_Free_Heap(l_array)) call 1289 E errquit(' basrtdbdoload: failed freeheap',0, MA_ERR) 1290c 1291c We messed around with the memory. The basis set sould have room for 1292c two (2) extra variables, which are added below. Add them now. 1293c 1294 h_tmp = exndcf(H_exndcf,basis) 1295 k_tmp = exndcf(K_exndcf,basis) 1296 ilen = exndcf(SZ_exndcf,basis) 1297 if (.not. ma_alloc_get(mt_dbl,ilen+2,nameex,h_new,k_new)) call 1298 & errquit('bas_rtdb_do_load: ma_alloc_get failed',ilen+2, 1299 & MA_ERR) 1300 exndcf(H_exndcf,basis) = h_new 1301 exndcf(K_exndcf,basis) = k_new 1302 exndcf(SZ_exndcf,basis) = ilen + 2 1303 call dcopy(ilen,dbl_mb(k_tmp),1,dbl_mb(k_new),1) 1304 if (.not.ma_free_heap(h_tmp)) call errquit 1305 & ('bas_rtdb_do_load: error freeing old exponents',211, 1306 & MEM_ERR) 1307c 1308 endif 1309c 1310 star_nr_tags = 0 1311c 1312c read the basis now get check status of read operations 1313c 1314 if (.not.rtdb_status) then 1315 if (exndcf(H_exndcf,basis) .ne. -1) then 1316 if (.not. ma_free_heap(exndcf(H_exndcf,basis))) 1317 $ call errquit('bas_rtdb_load: ma free failed?',0, 1318 & MA_ERR) 1319 exndcf(H_exndcf,basis) = -1 1320 exndcf(K_exndcf,basis) = 0 1321 exndcf(SZ_exndcf,basis)= 0 1322 endif 1323c 1324c rjh ... can be quiet now since the application should 1325c whine if it really did need the basis set 1326c 1327* if(ga_nodeid().eq.0) then 1328* write(LuOut,*) 1329* write(LuOut,*) ' bas_rtdb_do_load: basis not present "', 1330* $ bs_name(basis)(1:inp_strlen(bs_name(basis))), 1331* & '" -> "', 1332* $ bs_trans(basis)(1:inp_strlen(bs_trans(basis))), 1333* & '"' 1334* write(LuOut,*) 1335* endif 1336 bas_rtdb_do_load = .false. 1337c.....add diagnostics later 1338 return 1339 endif 1340c 1341c compute internal information and geobas maps 1342c 1343 bas_rtdb_do_load = bas_geobas_build(basisin) 1344c 1345c Add zero exponent S function for sp code/texas interface, 1346* incore information only 1347c 1348 if (bas_rtdb_do_load) then 1349 k_tmp = infbs_head(HEAD_EXCFPTR,basis) 1350 infbs_cont(CONT_TYPE, 0,basis) = 0 1351 infbs_cont(CONT_NPRIM,0,basis) = 1 1352 infbs_cont(CONT_NGEN, 0,basis) = 1 1353 infbs_cont(CONT_TAG, 0,basis) = -1 1354 k_tmp = k_tmp + 1 1355 infbs_cont(CONT_IEXP, 0,basis) = k_tmp 1356 dbl_mb(mb_exndcf(k_tmp,basis)) = 0.0d00 1357 k_tmp = k_tmp + 1 1358 infbs_cont(CONT_ICFP, 0,basis) = k_tmp 1359 dbl_mb(mb_exndcf(k_tmp,basis)) = 1.0d00 1360 infbs_head(HEAD_EXCFPTR,basis) = k_tmp 1361* 1362 call bas_ecce_print_basis(basisin,'load_basis') 1363 endif 1364 end 1365*..................................................................... 1366 logical function bas_geobas_build(basisin) 1367 implicit none 1368#include "errquit.fh" 1369c 1370#include "mafdecls.fh" 1371#include "rtdb.fh" 1372#include "nwc_const.fh" 1373#include "geomP.fh" 1374#include "geom.fh" 1375#include "context.fh" 1376#include "basP.fh" 1377#include "basdeclsP.fh" 1378#include "geobasmapP.fh" 1379#include "inp.fh" 1380#include "global.fh" 1381#include "bas_ibs_dec.fh" 1382#include "ecpso_decP.fh" 1383#include "stdio.fh" 1384c::passed 1385 integer basisin 1386c::local 1387 integer basis 1388 integer geom 1389 integer nat 1390 integer i, idum_cont, idum_at 1391 integer j, jstart, jend, jsize 1392 integer kstart, kend, ksize, lsize, icount 1393 integer nbf, iu_cont, myang 1394 integer my_gen, my_type 1395 integer atn 1396 character*16 element 1397 character*2 symbol 1398 logical status 1399 logical foundit 1400 logical found_any 1401 logical is_bq 1402 logical ecpORso, relbas 1403 integer int_dummy, num_elec 1404 integer h_tmp, k_tmp 1405 character*2 tag12 1406 character*16 name_tmp 1407 double precision erep_save 1408 integer uce, idbstag 1409*debug:mem integer inode 1410c::functions 1411 logical bas_high_angular 1412 external bas_high_angular 1413 integer nbf_from_ucont 1414 external nbf_from_ucont 1415 logical ecp_get_num_elec 1416 external ecp_get_num_elec 1417 logical bas_match_tags 1418 external bas_match_tags 1419 logical basis_is_rel 1420 external basis_is_rel 1421c 1422#include "bas_ibs_sfn.fh" 1423#include "ecpso_sfnP.fh" 1424c 1425 basis = basisin + Basis_Handle_Offset 1426 geom = ibs_geom(basis) 1427 ecpORso = Is_ECP(basis).or.Is_SO(basis) 1428 relbas = basis_is_rel(basisin) 1429c 1430*debug:mem do inode = 0,(ga_nnodes() - 1) 1431*debug:mem if (inode.eq.ga_nodeid()) call MA_summarize_allocated_blocks() 1432*debug:mem call ga_sync() 1433*debug:mem enddo 1434c 1435 1436 status = geom_ncent(geom, nat) 1437 if (nat.eq.0.or..not.status) then 1438 write(LuOut,*)' bas_geobas_build: ERROR ' 1439 write(LuOut,*)' number of centers is zero or weird' 1440 write(LuOut,*)' nat = ',nat 1441 bas_geobas_build = .false. 1442c..... add diagnostics later 1443 return 1444 endif 1445c.... set spherical flag 1446 if (infbs_head(HEAD_SPH,basis).eq.1) then 1447 bas_spherical(basis) = .true. 1448 else 1449 bas_spherical(basis) = .false. 1450 endif 1451c.... set flag if any general contractions are present 1452 bas_any_gc(basis) = .false. 1453 if(.not. ecpORso) then 1454 do i = 1,(infbs_head(Head_Ncont,basis)) 1455 if (.not.bas_any_gc(basis)) then 1456 my_gen = infbs_cont(Cont_Ngen,i,basis) 1457 my_type = infbs_cont(Cont_Type,i,basis) 1458 if (my_gen.gt.1.and.my_type.ge.0) 1459 & bas_any_gc(basis) = .true. 1460 endif 1461 enddo 1462 endif 1463c.... set flag if any sp (spd,spdf) shells are present 1464 bas_any_sp_shell(basis) = .false. 1465 if(.not. ecpORso) then 1466 do i = 1,(infbs_head(Head_Ncont,basis)) 1467 if (.not.bas_any_sp_shell(basis)) then 1468 my_gen = infbs_cont(Cont_Ngen,i,basis) 1469 my_type = infbs_cont(Cont_Type,i,basis) 1470 if (my_type.lt.0) then 1471 if (my_gen.ne.2) then 1472 write(luout,*)' sp shell with n_gen = ',my_gen 1473 call errquit ('bas_geobas_build: fatal error',911, 1474 & BASIS_ERR) 1475 endif 1476 bas_any_sp_shell(basis) = .true. 1477 endif 1478 endif 1479 enddo 1480 endif 1481c 1482c... clear old ibs_ce2uce if it exists 1483 if (ibs_ce2uce(SZ_ibs,basis).gt.0) then 1484*debug:mem write(LuOut,*)' clearing old ce2uce data ',ga_nodeid() 1485*debug:mem call util_flush(LuOut) 1486 h_tmp = ibs_ce2uce(H_ibs,basis) 1487 if (.not.ma_free_heap(h_tmp)) call errquit 1488 & ('bas_geobas_build: error freeing ibs_ce2uce',911, 1489 & BASIS_ERR) 1490 ibs_ce2uce(H_ibs,basis) = 0 1491 ibs_ce2uce(K_ibs,basis) = 0 1492 ibs_ce2uce(SZ_ibs,basis) = 0 1493 endif 1494c 123456789012 1495 write(name_tmp,'(a12,i4)') ' ibs_ce2uce ',basis 1496 if (.not.ma_alloc_get(mt_int,nat,name_tmp,h_tmp,k_tmp)) then 1497 call errquit 1498 & ('bas_geobas_build: error ma_alloc ibs_ce2uce',911, 1499 & MA_ERR) 1500 else 1501*debug:mem write(LuOut,*)' generating ce2uce data ',ga_nodeid() 1502*debug:mem call util_flush(LuOut) 1503 call ifill(nat,0,int_mb(k_tmp),1) 1504 ibs_ce2uce(H_ibs,basis) = h_tmp 1505 ibs_ce2uce(K_ibs,basis) = k_tmp 1506 ibs_ce2uce(SZ_ibs,basis) = nat 1507*debug:mem do inode = 0,(ga_nnodes() - 1) 1508*debug:mem if (inode.eq.ga_nodeid()) call MA_summarize_allocated_blocks() 1509*debug:mem call ga_sync() 1510*debug:mem call util_flush(LuOut) 1511*debug:mem enddo 1512 endif 1513c 1514c build center to unique center map 1515c 1516 if (nat.gt.nat_mx) then 1517 write(LuOut,*)' nat = ',nat 1518 write(LuOut,*)' nat max = ',nat_mx 1519 call errquit ('bas_geobas_build: nat.gt.nat_mx',911, BASIS_ERR) 1520 endif 1521 found_any = .false. 1522 do 00100 i=1,nat 1523 foundit = .false. 1524*before match_tags: do 00101 j = 1,infbs_head(HEAD_NTAGS,basis) 1525*before match_tags: if(bas_match_tags(tags(i,geom),bs_tags(j,basis))) then 1526*before match_tags: int_mb(mb_ibs_ce2uce(i,basis)) = j 1527*before match_tags: foundit = .true. 1528*before match_tags: goto 00102 1529*before match_tags: endif 1530*before match_tags:00101 continue 1531 if (bas_match_tags(tags(i,geom),basisin,j)) then 1532 int_mb(mb_ibs_ce2uce(i,basis)) = j 1533 foundit = .true. 1534 found_any = .true. 1535 goto 00102 1536 endif 1537 if (.not. foundit .and. .not. (ecpORso .or. relbas)) then 1538 if (geom_tag_to_element(tags(i,geom), symbol, element, 1539 $ atn)) then 1540 if (ga_nodeid().eq.0) 1541 & write(LuOut,10) 1542 & tags(i,geom)(1:inp_strlen(tags(i,geom))), 1543 $ element(1:inp_strlen(element)), 1544 $ bs_name(basis)(1:len_bs_name(basis)) 1545 10 format(/' ERROR: geometry tag ',a,' (',a, 1546 & ') is an atom ', 1547 $ 'but has no functions in basis "',a,'"'/ 1548 $ ' ERROR: only bq* centers can have no functions') 1549 if (ga_nodeid().eq.0) call util_flush(LuOut) 1550 call errquit 1551 & ('bas_geobas_build: basis/geometry mismatch', 0, 1552 & BASIS_ERR) 1553 else 1554 tag12 = tags(i,geom)(1:2) 1555 is_bq = inp_compare(.false.,'bq',tag12) .or. 1556 $ (inp_compare(.false.,'X',tags(i,geom)(1:1)) .and. 1557 $ (.not. inp_compare(.false.,'e',tags(i,geom)(2:2)))) 1558 if (ga_nodeid().eq.0 .and.(.not.is_bq)) 1559 & write(LuOut,11) i, 1560 & tags(i,geom)(1:inp_strlen(tags(i,geom))), 1561 $ bs_name(basis)(1:len_bs_name(basis)) 1562 11 format(/'WARNING: geometry tag ',i4, ' ', a, 1563 $ ' not found in basis "',a,'"'/) 1564 int_mb(mb_ibs_ce2uce(i,basis)) = 0 1565 if (ga_nodeid().eq.0) call util_flush(LuOut) 1566 endif 1567 endif 156800102 continue 156900100 continue 1570* 1571 if (.not.found_any) then 1572 if (ecpORso) then 1573 bas_geobas_build = .false. 1574 return 1575 else 1576 if (ga_nodeid().eq.0) then 1577 write(luout,*)' none of the geometry tags matched any ', 1578 & 'basis set tag in the basis "', 1579 & bs_name(basis)(1:len_bs_name(basis)), 1580 & '"' 1581 endif 1582 call errquit('bas_geobas_build: fatal error',911, 1583 & BASIS_ERR) 1584 endif 1585 endif 1586c 1587c build total # of contractions 1588c 1589*debug:mem do inode = 0,(ga_nnodes() - 1) 1590*debug:mem if (inode.eq.ga_nodeid()) then 1591*debug:mem call MA_summarize_allocated_blocks() 1592*debug:mem if (MA_verify_allocator_stuff()) then 1593*debug:mem write(LuOut,*)' no errors' 1594*debug:mem else 1595*debug:mem write(LuOut,*)' errors' 1596*debug:mem endif 1597*debug:mem endif 1598*debug:mem call ga_sync() 1599*debug:mem call util_flush(LuOut) 1600*debug:mem enddo 1601*debug:mem call bas_print_allocated_info('bas_geobas_build 1') 1602 ncont_tot_gb(basis) = 0 1603 do 10200 i=1,nat 1604 uce = sf_ibs_ce2uce(i,basis) 1605*debug:mem write(LuOut,*)' myuce = ',uce,i,ga_nodeid() 1606*debug:mem call util_flush(LuOut) 1607 if (uce.gt.0) then 1608 idum_cont = infbs_tags(TAG_NCONT,uce,basis) 1609 ncont_tot_gb(basis) = idum_cont + ncont_tot_gb(basis) 1610 endif 161110200 continue 1612c 1613c allocate space for center -> contraction range map 1614c 1615c... clear old ibs_ce2cnr if it exists 1616 if (ibs_ce2cnr(SZ_ibs,basis).gt.0) then 1617*debug:mem write(LuOut,*)' clearing old ce2cnr data ',ga_nodeid() 1618*debug:mem call util_flush(LuOut) 1619 h_tmp = ibs_ce2cnr(H_ibs,basis) 1620 if (.not.ma_free_heap(h_tmp)) call errquit 1621 & ('bas_geobas_build: error freeing ibs_ce2cnr',911, 1622 & MEM_ERR) 1623 ibs_ce2cnr(H_ibs,basis) = 0 1624 ibs_ce2cnr(K_ibs,basis) = 0 1625 ibs_ce2cnr(SZ_ibs,basis) = 0 1626 endif 1627c 123456789012 1628 write(name_tmp,'(a12,i4)') ' ibs_ce2cnr ',basis 1629 if (.not.ma_alloc_get(mt_int,(2*nat), 1630 & name_tmp,h_tmp,k_tmp)) then 1631 call errquit 1632 & ('bas_geobas_build: error ma_alloc ibs_ce2cnr',911, MA_ERR) 1633 else 1634*debug:mem write(LuOut,*)' generating ce2cnr data ',ga_nodeid() 1635*debug:mem call util_flush(LuOut) 1636 call ifill((2*nat),0,int_mb(k_tmp),1) 1637 ibs_ce2cnr(H_ibs,basis) = h_tmp 1638 ibs_ce2cnr(K_ibs,basis) = k_tmp 1639 ibs_ce2cnr(SZ_ibs,basis) = 2*nat 1640*debug:mem do inode = 0,(ga_nnodes() - 1) 1641*debug:mem if (inode.eq.ga_nodeid()) call MA_summarize_allocated_blocks() 1642*debug:mem call ga_sync() 1643*debug:mem call util_flush(LuOut) 1644*debug:mem enddo 1645 endif 1646c 1647c build center -> contraction range map 1648c 1649 int_dummy = 0 1650 do 00200 i=1,nat 1651 if (sf_ibs_ce2uce(i,basis).gt.0) then 1652 idum_cont = 1653 & infbs_tags(TAG_NCONT,sf_ibs_ce2uce(i,basis),basis) 1654 int_mb(mb_ibs_ce2cnr(1,i,basis)) = int_dummy + 1 1655 int_mb(mb_ibs_ce2cnr(2,i,basis)) = int_dummy + idum_cont 1656 int_dummy = idum_cont + int_dummy 1657 else 1658*. . . . . . . . . . . . . . . . . . . . ! No functions on this center 1659 int_mb(mb_ibs_ce2cnr(1,i,basis)) = 0 1660 int_mb(mb_ibs_ce2cnr(2,i,basis)) = -1 1661 endif 166200200 continue 1663c 1664 if (ncont_tot_gb(basis) .eq. 0) call errquit 1665 $ ('bas_geobas_build: no functions in basis set', 0, 1666 & BASIS_ERR) 1667 if (ncont_tot_gb(basis) .gt. ncont_mx) then 1668 write(LuOut,*)' number of contractions = ', 1669 & ncont_tot_gb(basis) 1670 write(LuOut,*)' number of contractions max = ',ncont_mx 1671 call errquit ('bas_geobas_build: ncont.gt.ncont_mx ',911, 1672 & BASIS_ERR) 1673 endif 1674c 1675c allocate space for contraction -> center map 1676c 1677c... clear old ibs_cn2ce if it exists 1678 if (ibs_cn2ce(SZ_ibs,basis).gt.0) then 1679*debug:mem write(LuOut,*)' clearing old cn2ce data ',ga_nodeid() 1680*debug:mem call util_flush(LuOut) 1681 h_tmp = ibs_cn2ce(H_ibs,basis) 1682 if (.not.ma_free_heap(h_tmp)) call errquit 1683 & ('bas_geobas_build: error freeing ibs_cn2ce',911, 1684 & MEM_ERR) 1685 ibs_cn2ce(H_ibs,basis) = 0 1686 ibs_cn2ce(K_ibs,basis) = 0 1687 ibs_cn2ce(SZ_ibs,basis) = 0 1688 endif 1689c 123456789012 1690 write(name_tmp,'(a12,i4)') ' ibs_cn2ce ',basis 1691 if (.not.ma_alloc_get(mt_int,(1+ncont_tot_gb(basis)), 1692 & name_tmp,h_tmp,k_tmp)) then 1693 call errquit('bas_geobas_build: error ma_alloc ibs_cn2ce',911, 1694 & MA_ERR) 1695 else 1696*debug:mem write(LuOut,*)' generating cn2ce data ',ga_nodeid() 1697*debug:mem call util_flush(LuOut) 1698 call ifill((1+ncont_tot_gb(basis)),0,int_mb(k_tmp),1) 1699 ibs_cn2ce(H_ibs,basis) = h_tmp 1700 ibs_cn2ce(K_ibs,basis) = k_tmp 1701 ibs_cn2ce(SZ_ibs,basis) = 1+ncont_tot_gb(basis) 1702 endif 1703c 1704c build contraction -> center map 1705c 1706 do 00300 i=1,nat 1707 if (sf_ibs_ce2uce(i,basis).gt.0) then 1708 jstart = sf_ibs_ce2cnr(1,i,basis) 1709 jend = sf_ibs_ce2cnr(2,i,basis) 1710 do 00400 j=jstart,jend 1711 int_mb(mb_ibs_cn2ce(j,basis)) = i 171200400 continue 1713 endif 171400300 continue 1715c 1716c set zero element of cn2ce to something useless 1717c 1718 int_mb(mb_ibs_cn2ce(0,basis)) = -1 1719c 1720c allocate space for ibs_cn2ucn map 1721c 1722c... clear old ibs_cn2ucn if it exists 1723 if (ibs_cn2ucn(SZ_ibs,basis).gt.0) then 1724*debug:mem write(LuOut,*)' clearing old cn2ucn data ',ga_nodeid() 1725*debug:mem call util_flush(LuOut) 1726 h_tmp = ibs_cn2ucn(H_ibs,basis) 1727 if (.not.ma_free_heap(h_tmp)) call errquit 1728 & ('bas_geobas_build: error freeing ibs_cn2ucn',911, 1729 & BASIS_ERR) 1730 ibs_cn2ucn(H_ibs,basis) = 0 1731 ibs_cn2ucn(K_ibs,basis) = 0 1732 ibs_cn2ucn(SZ_ibs,basis) = 0 1733 endif 1734c 123456789012 1735 write(name_tmp,'(a12,i4)') ' ibs_cn2ucn ',basis 1736 if (.not.ma_alloc_get(mt_int,(1+ncont_tot_gb(basis)), 1737 & name_tmp,h_tmp,k_tmp)) then 1738 call errquit 1739 & ('bas_geobas_build: error ma_alloc ibs_cn2ucn',911, 1740 & MA_ERR) 1741 else 1742*debug:mem write(LuOut,*)' generating cn2ucn data ',ga_nodeid() 1743*debug:mem call util_flush(LuOut) 1744 call ifill((1+ncont_tot_gb(basis)),0,int_mb(k_tmp),1) 1745 ibs_cn2ucn(H_ibs,basis) = h_tmp 1746 ibs_cn2ucn(K_ibs,basis) = k_tmp 1747 ibs_cn2ucn(SZ_ibs,basis) = 1+ncont_tot_gb(basis) 1748 endif 1749c 1750c build contraction -> unique contraction map 1751c 1752 do 00500 i=1,nat 1753 jstart = sf_ibs_ce2cnr(1,i,basis) 1754 jend = sf_ibs_ce2cnr(2,i,basis) 1755 jsize = jend - jstart + 1 1756 if (jsize .gt. 0) then 1757 idum_at = sf_ibs_ce2uce(i,basis) 1758 kstart = infbs_tags(TAG_FCONT,idum_at,basis) 1759 kend = infbs_tags(TAG_LCONT,idum_at,basis) 1760 ksize = kend - kstart + 1 1761 lsize = infbs_tags(TAG_NCONT,idum_at,basis) 1762 if (jsize.eq.ksize.and.ksize.eq.lsize) then 1763 icount = 0 1764 do 00600 j=jstart,jend 1765 int_mb(mb_ibs_cn2ucn(j,basis)) = kstart + icount 1766 icount = icount + 1 176700600 continue 1768 else 1769 write(LuOut,*)' bas_geobas_build: ERROR ' 1770 write(LuOut,*)' contraction range size mismatch' 1771 write(LuOut,*)' cont. range (',jstart,':',jend,')' 1772 write(LuOut,*)' unique cont. range (',kstart,':',kend,')' 1773 write(LuOut,*)' cont. size: ',jsize 1774 write(LuOut,*)' calculated unique cont. size: ',ksize 1775 write(LuOut,*)' lookup unique cont. size: ',lsize 1776 bas_geobas_build = .false. 1777 return 1778 endif 1779 endif 178000500 continue 1781c 1782c set zero element 1783c 1784 int_mb(mb_ibs_cn2ucn(0,basis)) = 0 1785c 1786c allocate space for ibs_cn2bfr map 1787c 1788c... clear old ibs_cn2bfr if it exists 1789 if (ibs_cn2bfr(SZ_ibs,basis).gt.0) then 1790*debug:mem write(LuOut,*)' clearing old cn2bfr data ',ga_nodeid() 1791*debug:mem call util_flush(LuOut) 1792 h_tmp = ibs_cn2bfr(H_ibs,basis) 1793 if (.not.ma_free_heap(h_tmp)) call errquit 1794 & ('bas_geobas_build: error freeing ibs_cn2bfr',911, 1795 & BASIS_ERR) 1796 ibs_cn2bfr(H_ibs,basis) = 0 1797 ibs_cn2bfr(K_ibs,basis) = 0 1798 ibs_cn2bfr(SZ_ibs,basis) = 0 1799 endif 1800c 123456789012 1801 write(name_tmp,'(a12,i4)') ' ibs_cn2bfr ',basis 1802 if (.not.ma_alloc_get(mt_int,(2*(1+ncont_tot_gb(basis))), 1803 & name_tmp,h_tmp,k_tmp)) then 1804 call errquit 1805 & ('bas_geobas_build: error ma_alloc ibs_cn2bfr',911, 1806 & BASIS_ERR) 1807 else 1808*debug:mem write(LuOut,*)' generating cn2bfr data ',ga_nodeid() 1809*debug:mem call util_flush(LuOut) 1810 call ifill((2*(1+ncont_tot_gb(basis))),0,int_mb(k_tmp),1) 1811 ibs_cn2bfr(H_ibs,basis) = h_tmp 1812 ibs_cn2bfr(K_ibs,basis) = k_tmp 1813 ibs_cn2bfr(SZ_ibs,basis) = 2*(1+ncont_tot_gb(basis)) 1814 endif 1815c 1816c build nprim_tot_gb, ncoef_tot_gb, nbf_tot_gb, and 1817c contraction -> basis function range map 1818c find nbfmax for basis (initialized in block data statement) 1819c 1820 nbf_tot_gb(basis) = 0 1821 nprim_tot_gb(basis) = 0 1822 do 00700 i = 1,ncont_tot_gb(basis) 1823 iu_cont = sf_ibs_cn2ucn(i,basis) 1824c 1825 nbf = nbf_from_ucont(iu_cont,basisin) 1826 nbfmax_bs(basis) = max(nbfmax_bs(basis),nbf) 1827c 1828 int_mb(mb_ibs_cn2bfr(1,i,basis)) = nbf_tot_gb(basis) + 1 1829 int_mb(mb_ibs_cn2bfr(2,i,basis)) = nbf_tot_gb(basis) + nbf 1830 1831 nbf_tot_gb(basis) = nbf_tot_gb(basis) + nbf 1832 nprim_tot_gb(basis) = nprim_tot_gb(basis) + 1833 & infbs_cont(CONT_NPRIM,iu_cont,basis) 1834 ncoef_tot_gb(basis) = ncoef_tot_gb(basis) + 1835 & infbs_cont(CONT_NPRIM,iu_cont,basis)* 1836 & infbs_cont(CONT_NGEN,iu_cont,basis) 183700700 continue 1838c 1839c set zero elements of cn2bfr 1840c 1841 int_mb(mb_ibs_cn2bfr(1,0,basis)) = 0 1842 int_mb(mb_ibs_cn2bfr(2,0,basis)) = 0 1843c 1844c build high angular momentum of this loaded <basis|geom> pair 1845c note angular_bs(*) initialized in block data function 1846c 1847 if (.not.bas_high_angular(basisin,myang))call errquit 1848 & ('bas_geobas_build: error bas_high_angular',911, 1849 & BASIS_ERR) 1850* 1851 if (Is_ECP(basis)) then 1852* if ECP basis must modify geometry data appropriately 1853 do i = 1,nat 1854 oecpcent(i,geom)=.false. 1855 idbstag = sf_ibs_ce2uce(i,basis) 1856 if (ecp_get_num_elec(basisin, 1857 & bs_tags(idbstag,basis),num_elec)) then 1858 oecpcent(i,geom) = .true. 1859 charge(i,geom) = charge(i,geom) - dble(num_elec) 1860 endif 1861 enddo 1862 erep_save = erep(geom) 1863 call geom_compute_values(geom) 1864* 1865* cannot print this here 1866* 1867* write(luout,*) 1868* & ' nuclear repulsion energy without ECP: ',erep_save 1869* write(luout,*) 1870* & ' nuclear repulsion energy with ECP: ',erep(geom) 1871 endif 1872* 1873 bas_geobas_build = .true. 1874* 1875 end 1876*..................................................................... 1877C> \ingroup bas 1878C> @{ 1879c 1880C> \brief Store a basis set on the RTDB under the given name 1881c 1882C> \return Return .true. if successfull, and .false. otherwise 1883c 1884 logical function bas_rtdb_store(rtdb, name, basisin) 1885 implicit none 1886c 1887c routine that does an incore basis set store to the rtdb 1888c 1889#include "mafdecls.fh" 1890#include "basdeclsP.fh" 1891#include "nwc_const.fh" 1892#include "basP.fh" 1893#include "ecpso_decP.fh" 1894#include "bas_starP.fh" 1895c::passed 1896 integer rtdb !< [Input] the RTDB handle 1897 character*(*) name !< [Input] the name to use when storing basis 1898 integer basisin !< [Input] the basis set handle 1899c::functions 1900 logical bas_check_handle, bas_rtdb_do_store 1901 external bas_check_handle, bas_rtdb_do_store 1902c 1903c Store basis set (not geometry) related info about specified 1904c basis in into the rtdb with the given name 1905c 1906c::local 1907 integer basis ! Actual index into basis set arrays 1908 integer size_ex 1909c 1910cc AJL/Begin/SPIN-POLARISED ECPs 1911 integer channels 1912 logical ecp_get_high_chan 1913 external ecp_get_high_chan 1914cc AJL/End 1915c 1916c:: statement functions 1917#include "errquit.fh" 1918#include "bas_exndcf.fh" 1919#include "ecpso_sfnP.fh" 1920c 1921cc AJL/Begin/SPIN-POLARISED ECPs 1922 channels = 1 1923cc AJL/End 1924c 1925 bas_rtdb_store = bas_check_handle(basisin,'bas_rtdb_store') 1926 if (.not. bas_rtdb_store) return 1927 basis = basisin + Basis_Handle_Offset 1928c 1929 if (Is_ECP(basis).or.Is_SO(basis)) then 1930 size_ex = (2*infbs_head(HEAD_NPRIM,basis))+ 1931 & infbs_head(HEAD_NCOEF,basis) 1932c 1933cc AJL/Begin/SPIN-POLARISED ECPs 1934cc Use this to put a note of spin_polarised_ecps in rtdb 1935cc For use in nwdft/dft_fockbld.F 1936cc 1937cc By using this function we deal with anomalies for all 1938cc applications of the spin-polarised ECPs being for both channels, 1939cc as ecp_get_high_chan will treat that as spin-independent 1940cc calculation of the fock matrices 1941 if (Is_ECP(basis)) then 1942 if (.not.ecp_get_high_chan(basisin,channels)) 1943 & call errquit('bas_rtdb_store error',911, BASIS_ERR) 1944 endif 1945cc AJL/End 1946c 1947 else 1948 size_ex = infbs_head(HEAD_NPRIM,basis)+ 1949 & infbs_head(HEAD_NCOEF,basis) 1950 endif 1951 bas_rtdb_store = bas_rtdb_do_store(rtdb, name, 1952 $ bs_tags(1,basis), infbs_head(1,basis), 1953 $ infbs_tags(1,1,basis), 1954 $ infbs_cont(1,1,basis), 1955 & dbl_mb(mb_exndcf(1,basis)), 1956 $ infbs_head(HEAD_NTAGS,basis), 1957 $ infbs_head(HEAD_NCONT,basis), 1958 $ size_ex, 1959 & bs_stdname(1,basis), 1960 & name_assoc(1,basis), 1961 & channels, ! AJL/SPIN-POLARISED ECPs 1962 & star_nr_tags, star_tag, star_in_lib, star_bas_typ, 1963 & star_file, star_nr_excpt, star_excpt, star_tot_excpt, 1964 & star_rel, star_segment) 1965c 1966 end 1967C> @} 1968*..................................................................... 1969 logical function bas_rtdb_do_store(rtdb, name, tagsin, 1970 & head_array, tags_array, ucont_array, excfin, ntagsin, 1971 & nucontin, nexcf, stdnamein, ecp_name, ecp_channels, 1972 & star_nr_tags, star_tag, star_in_lib, star_bas_typ, 1973 & star_file, star_nr_excpt, star_excpt, star_tot_excpt, 1974 & star_rel, star_segment) 1975 implicit none 1976c 1977c stores from argument data structures the basis set information to 1978c the rtdb. 1979c 1980#include "mafdecls.fh" 1981#include "basdeclsP.fh" 1982#include "nwc_const.fh" 1983#include "basP.fh" 1984#include "inp.fh" 1985#include "rtdb.fh" 1986#include "stdio.fh" 1987c 1988c This routine stores the basis set information in the appropriate 1989c data structure on the run-time-data-base (rtdb). 1990c 1991c This is a private routine called by the user level routine 1992c bas_rtdb_store(rtdb, name, basis) 1993c 1994c::: functions 1995 logical bas_rtdb_in 1996 logical bas_rtdb_add 1997 external bas_rtdb_in 1998 external bas_rtdb_add 1999c::: passed 2000 integer rtdb ! [input] rtdb handle 2001 character*(*) name ! [input] name of basis set 2002 integer ntagsin ! [input] number of tags 2003 integer nucontin ! [input] number of unique contractions 2004 integer nexcf ! [input] number of exponents and 2005 character*16 tagsin(ntagsin) ! [input] name of tags 2006 character*80 stdnamein(ntagsin) ! [input] names of basis on a tag 2007 integer head_array(ndbs_head) ! [input] head data 2008 integer tags_array(ndbs_tags,ntagsin) ! [input] tag data 2009 integer ucont_array(ndbs_ucont,nucontin) ! [input] unique 2010*. . . . . . . . . . . . . . . . . . . . . . . contraction data 2011 character*(*) ecp_name(2) ! [input] ecp 2012*. . . . . . . . . . . . . . . . . . . . . . . associated with 2013*. . . . . . . . . . . . . . . . . . . . . . . normal basis 2014 double precision excfin(nexcf) ! [input] exponents 2015c. . . . . . . . . . . . . . . . . . . contractions coeffs. 2016 integer star_nr_tags ! [input] number of star containing tags 2017 character*16 star_tag(star_nr_tags) ! [input] star tag data 2018 character*16 star_in_lib(star_nr_tags) ! [input] star in lib data 2019 character*255 star_bas_typ(star_nr_tags) ! [input] star basis data 2020 character*255 star_file(star_nr_tags) ! [input] star filename data 2021 logical star_rel(star_nr_tags) ! [input] star relativistic data 2022 logical star_segment ! [input] star segment or not ? 2023 integer star_tot_excpt ! [input] star total number of except tags 2024 integer star_nr_excpt(star_nr_tags) ! [input] number of except per tag 2025 character*16 star_excpt(star_nr_tags) ! [input] except tag list 2026c 2027cc AJL/Begin/SPIN-POLARISED ECPs 2028 integer ecp_channels ! [input] number of channels for ECP (1 or 2). 2029c............................ *Should* be 1 in all other scenarios 2030cc AJL/End 2031c 2032c::: local 2033 character*256 tmp 2034 integer len_name, lentmp 2035 logical status 2036c 2037 bas_rtdb_do_store = .true. 2038c 2039 status = bas_rtdb_in(rtdb) 2040c 2041c generate rtdb names and store information 2042c 2043 len_name = inp_strlen(name) 2044 tmp = 'basis:'//name(1:len_name) 2045 lentmp = inp_strlen(tmp) + 1 2046 2047 status = .true. 2048 status = status.and.bas_rtdb_add(rtdb,name) 2049 tmp(lentmp:) = ' ' 2050 tmp(lentmp:) = ':bs_nr_tags' 2051 status = status .and. rtdb_put(rtdb,tmp,mt_int,1,ntagsin) 2052 if (ntagsin .gt. 0) then 2053 tmp(lentmp:) = ' ' 2054 tmp(lentmp:) = ':bs_tags' 2055 status = status .and. rtdb_cput(rtdb, 2056 & tmp,ntagsin,tagsin) 2057 2058 tmp(lentmp:) = ' ' 2059 tmp(lentmp:) = ':bs_stdname' 2060 status = status .and. rtdb_cput(rtdb, 2061 & tmp,ntagsin,stdnamein) 2062 endif 2063 2064 tmp(lentmp:) = ' ' 2065 tmp(lentmp:) = ':assoc ecp name' 2066 status = status .and. rtdb_cput(rtdb,tmp,2,ecp_name) 2067c 2068 tmp(lentmp:) = ' ' 2069 tmp(lentmp:) = ':number of exps and coeffs' 2070 status = status .and. rtdb_put(rtdb,tmp,mt_int,1,nexcf) 2071 2072 if (nexcf .gt. 0) then 2073 tmp(lentmp:) = ' ' 2074 tmp(lentmp:) = ':exps and coeffs' 2075 status = status .and. rtdb_put(rtdb,tmp,mt_dbl,nexcf,excfin) 2076 endif 2077 2078 tmp(lentmp:) = ' ' 2079 tmp(lentmp:) = ':header' 2080 status = status .and. 2081 & rtdb_put(rtdb,tmp,mt_int,ndbs_head,head_array) 2082 2083 tmp(lentmp:) = ' ' 2084 tmp(lentmp:) = ':tags info' 2085 status = status .and. rtdb_put( 2086 & rtdb,tmp,mt_int,(ndbs_tags*ntags_bsmx), tags_array) 2087 2088 tmp(lentmp:) = ' ' 2089 tmp(lentmp:) = ':contraction info' 2090 status = status .and. rtdb_put( 2091 & rtdb,tmp,mt_int,(ndbs_ucont*nucont_bsmx),ucont_array) 2092 tmp(lentmp:) = ' ' 2093 tmp(lentmp:) = ':star nr tags' 2094 status = status .and. rtdb_put( 2095 & rtdb,tmp,mt_int,1,star_nr_tags) 2096c 2097c store remaining star tag info only if we have actually a star 2098c tag input (hence, if star_nr_tags > 0) 2099c 2100 if (star_nr_tags .gt. 0) then 2101 tmp(lentmp:) = ' ' 2102 tmp(lentmp:) = ':star tag names' 2103 status = status .and. rtdb_cput(rtdb,tmp,star_nr_tags, 2104 & star_tag) 2105 tmp(lentmp:) = ' ' 2106 tmp(lentmp:) = ':star tag_in_lib' 2107 status = status .and. rtdb_cput(rtdb,tmp,star_nr_tags, 2108 & star_in_lib) 2109 tmp(lentmp:) = ' ' 2110 tmp(lentmp:) = ':star bas type' 2111 status = status .and. rtdb_cput(rtdb,tmp,star_nr_tags, 2112 & star_bas_typ) 2113 tmp(lentmp:) = ' ' 2114 tmp(lentmp:) = ':star filename' 2115 status = status .and. rtdb_cput(rtdb,tmp,star_nr_tags, 2116 & star_file) 2117 tmp(lentmp:) = ' ' 2118 tmp(lentmp:) = ':star tot except' 2119 status = status .and. rtdb_put( 2120 & rtdb,tmp,mt_int,1,star_tot_excpt) 2121 if ( star_tot_excpt .gt. 0) then 2122 tmp(lentmp:) = ' ' 2123 tmp(lentmp:) = ':star except' 2124 status = status .and. rtdb_cput(rtdb,tmp,star_tot_excpt, 2125 & star_excpt) 2126 endif 2127 tmp(lentmp:) = ' ' 2128 tmp(lentmp:) = ':star nr except' 2129 status = status .and. rtdb_put( 2130 & rtdb,tmp,mt_int,star_nr_tags,star_nr_excpt) 2131 tmp(lentmp:) = ' ' 2132 tmp(lentmp:) = ':star rel' 2133 status = status .and. rtdb_put( 2134 & rtdb,tmp,mt_log,star_nr_tags,star_rel) 2135 tmp(lentmp:) = ' ' 2136 tmp(lentmp:) = ':star segment' 2137 status = status .and. rtdb_put( 2138 & rtdb,tmp,mt_log,1,star_segment) 2139 endif 2140c 2141c reset the stag tag data array for reuse 2142c 2143 star_nr_tags = 0 2144c 2145cc AJL/Begin/SPIN-POLARISED ECPs 2146cc Add if we are using spin_polarised_ecps 2147cc There may be a better place for this information 2148 if (ecp_channels.gt.1) then 2149 status = status .and. rtdb_put( 2150 & rtdb,'dft:spin_polarised_ecps',mt_int,1,ecp_channels) 2151 endif 2152cc AJL/End 2153c 2154c read the basis now get check status of read operations 2155c 2156 if (.not.status) then 2157 write(LuOut,*)' bas_rtdb_store: ERROR ' 2158 write(LuOut,*)' one or more put operations failed ' 2159 bas_rtdb_do_store = .false. 2160c..... add diagnostics later 2161 return 2162 endif 2163 return 2164 end 2165*..................................................................... 2166C> \ingroup bas 2167C> @{ 2168c 2169C> \brief Retrieve the highest angular momentum in the specified basis 2170c 2171C> \return Return .true. if successfull, and .false. otherwise 2172c 2173 logical function bas_high_angular(basisin,high_angular) 2174 implicit none 2175c 2176c calculate, return and store high angular momentem function 2177c for given basis. 2178c 2179#include "nwc_const.fh" 2180#include "basP.fh" 2181#include "basdeclsP.fh" 2182#include "stdio.fh" 2183c::functions 2184 logical bas_check_handle 2185 external bas_check_handle 2186c:: passed 2187 integer basisin !< [Input] the basis set handle 2188 integer high_angular !< [Output] the maximum angular momentum 2189c:local 2190 integer basis, i, myang 2191 integer myutag 2192c 2193 bas_high_angular = bas_check_handle(basisin,'bas_high_angular') 2194 if (.not. bas_high_angular ) then 2195 write(luout,*) 'bas_high_angular: basis handle not valid ' 2196 return 2197 endif 2198c 2199 basis = basisin + Basis_Handle_Offset 2200 if (angular_bs(basis) .gt. -565) then 2201 high_angular = angular_bs(basis) 2202 bas_high_angular = .true. 2203 return 2204 endif 2205c 2206 myutag = infbs_head(head_ntags,basis) 2207 high_angular = -565 2208 do i = 1,myutag 2209 myang = abs(infbs_tags(tag_high_ang,i,basis)) 2210 high_angular = max(high_angular,myang) 2211 enddo 2212*is this needed here? 2213 angular_bs(basis) = high_angular 2214c 2215 bas_high_angular = .true. 2216 return 2217 end 2218C> @} 2219*..................................................................... 2220cc AJL/Begin 2221C> \ingroup ecp 2222C> @{ 2223c 2224C> \brief Retrieve the number of channels in the ecp basis 2225c 2226C> \return Return .true. if successfull, and .false. otherwise 2227c 2228 logical function ecp_get_high_chan(ecpidin,channels) 2229 implicit none 2230c 2231c calculate and return high channel for given ecp. 2232c 2233#include "nwc_const.fh" 2234#include "basP.fh" 2235#include "basdeclsP.fh" 2236#include "stdio.fh" 2237c::functions 2238 logical ecp_check_handle 2239 external ecp_check_handle 2240c:: passed 2241 integer ecpidin !< [Input] the ECP set handle 2242 integer channels !< [Output] the number of channels 2243c:local 2244 integer ecpid, i, j, mychannel 2245 integer myucont, mytags 2246c 2247 ecpid = ecpidin + Basis_Handle_Offset 2248c 2249 ecp_get_high_chan = ecp_check_handle(ecpidin,'ecp_get_high_chan') 2250 if (.not. ecp_get_high_chan) then 2251 write (6,*) 'This is not an ECP Basis' 2252 return 2253 endif 2254c 2255 mytags = infbs_head(HEAD_NTAGS,ecpid) 2256 if (mytags.le.0) return 2257c 2258 channels = 0 2259 do i=1,mytags 2260 myucont = infbs_tags(TAG_NCONT,i,ecpid) 2261 do j = 1,myucont 2262 mychannel = infbs_cont(CONT_CHANNEL,j,ecpid) 2263 channels = max(channels,mychannel) 2264 enddo 2265 enddo 2266c If all values are 0 we have a spin independent calculation 2267 if (channels.eq.0) then 2268 channels = 1 2269c Somewhere we have initialised spin in the alpha of beta channels 2270 else if (channels.gt.0) then 2271 channels = 2 2272 end if 2273 2274c Save value process to be added 2275c 2276 return 2277 end 2278C> @} 2279cc AJL/End 2280*..................................................................... 2281 logical function gbs_map_clear(basisin) 2282 implicit none 2283#include "errquit.fh" 2284#include "mafdecls.fh" 2285#include "nwc_const.fh" 2286#include "basP.fh" 2287#include "geobasmapP.fh" 2288#include "bas_ibs_dec.fh" 2289#include "stdio.fh" 2290c 2291c routine to clear online map information and basis information 2292c 2293c::functions 2294 logical bas_check_handle 2295 external bas_check_handle 2296c:util 2297c ifill 2298c::passed 2299 integer basisin ! [input] basis set handle 2300c::local 2301 integer basis 2302 integer h_tmp 2303c 2304#include "bas_ibs_sfn.fh" 2305c 2306 gbs_map_clear = bas_check_handle(basisin,'gbs_map_clear') 2307 if (.not. gbs_map_clear ) then 2308 write(LuOut,*) ' basis handle not valid ' 2309 return 2310 endif 2311c 2312 basis = basisin + Basis_Handle_Offset 2313c... clear old ibs_ce2uce if it exists 2314 if (ibs_ce2uce(SZ_ibs,basis).gt.0) then 2315 h_tmp = ibs_ce2uce(H_ibs,basis) 2316 if (.not.ma_free_heap(h_tmp)) call errquit 2317 & ('gbs_map_clear: error freeing ibs_ce2uce',911, MEM_ERR) 2318 ibs_ce2uce(H_ibs,basis) = 0 2319 ibs_ce2uce(K_ibs,basis) = 0 2320 ibs_ce2uce(SZ_ibs,basis) = 0 2321 endif 2322c... clear old ibs_ce2cnr if it exists 2323 if (ibs_ce2cnr(SZ_ibs,basis).gt.0) then 2324 h_tmp = ibs_ce2cnr(H_ibs,basis) 2325 if (.not.ma_free_heap(h_tmp)) call errquit 2326 & ('gbs_map_clear: error freeing ibs_ce2cnr',911, MEM_ERR) 2327 ibs_ce2cnr(H_ibs,basis) = 0 2328 ibs_ce2cnr(K_ibs,basis) = 0 2329 ibs_ce2cnr(SZ_ibs,basis) = 0 2330 endif 2331c... clear old ibs_cn2ce if it exists 2332 if (ibs_cn2ce(SZ_ibs,basis).gt.0) then 2333 h_tmp = ibs_cn2ce(H_ibs,basis) 2334 if (.not.ma_free_heap(h_tmp)) call errquit 2335 & ('gbs_map_clear: error freeing ibs_cn2ce',911, MEM_ERR) 2336 ibs_cn2ce(H_ibs,basis) = 0 2337 ibs_cn2ce(K_ibs,basis) = 0 2338 ibs_cn2ce(SZ_ibs,basis) = 0 2339 endif 2340c... clear old ibs_cn2ucn if it exists 2341 if (ibs_cn2ucn(SZ_ibs,basis).gt.0) then 2342 h_tmp = ibs_cn2ucn(H_ibs,basis) 2343 if (.not.ma_free_heap(h_tmp)) call errquit 2344 & ('gbs_map_clear: error freeing ibs_cn2ucn',911, MEM_ERR) 2345 ibs_cn2ucn(H_ibs,basis) = 0 2346 ibs_cn2ucn(K_ibs,basis) = 0 2347 ibs_cn2ucn(SZ_ibs,basis) = 0 2348 endif 2349c... clear old ibs_cn2bfr if it exists 2350 if (ibs_cn2bfr(SZ_ibs,basis).gt.0) then 2351 h_tmp = ibs_cn2bfr(H_ibs,basis) 2352 if (.not.ma_free_heap(h_tmp)) call errquit 2353 & ('gbs_map_clear: error freeing ibs_cn2bfr',911, MEM_ERR) 2354 ibs_cn2bfr(H_ibs,basis) = 0 2355 ibs_cn2bfr(K_ibs,basis) = 0 2356 ibs_cn2bfr(SZ_ibs,basis) = 0 2357 endif 2358c 2359 call ifill(3,0,ibs_cn2ucn(1,basis),1) 2360 call ifill(3,0,ibs_cn2ce (1,basis),1) 2361 call ifill(3,0,ibs_cn2bfr(1,basis),1) 2362 call ifill(3,0,ibs_ce2uce(1,basis),1) 2363 call ifill(3,0,ibs_ce2cnr(1,basis),1) 2364 ncont_tot_gb(basis) = 0 2365 nprim_tot_gb(basis) = 0 2366 ncoef_tot_gb(basis) = 0 2367 nbf_tot_gb(basis) = 0 2368c 2369 gbs_map_clear = .true. 2370 return 2371 end 2372*..................................................................... 2373 logical function gbs_map_print(basisin) 2374 implicit none 2375c 2376c prints the basis set <-> geometry mapping information 2377c 2378#include "mafdecls.fh" 2379#include "nwc_const.fh" 2380#include "basP.fh" 2381#include "geobasmapP.fh" 2382#include "basdeclsP.fh" 2383#include "geom.fh" 2384#include "bas_ibs_dec.fh" 2385#include "stdio.fh" 2386c::function 2387 logical bas_check_handle 2388 external bas_check_handle 2389c::passed 2390 integer basisin ! [input] basis set handle 2391c::local 2392 integer nat, basis, i 2393 integer myfirst, mylast, mysize, mycenter, myucont 2394 integer mygeom 2395 logical status 2396c 2397#include "bas_ibs_sfn.fh" 2398c 2399 basis = basisin + Basis_Handle_Offset 2400c 2401c check geom and basis handles returns false if either is invalid 2402c 2403 mygeom = ibs_geom(basis) 2404 gbs_map_print=geom_check_handle(mygeom,'gbs_map_print') 2405 if (.not.gbs_map_print) return 2406 gbs_map_print=bas_check_handle(basisin,'gbs_map_print') 2407 if (.not.gbs_map_print) return 2408c 2409c find number of atoms 2410 status = geom_ncent(mygeom, nat) 2411c 2412 if (nat.eq.0.or..not.status) then 2413 write(LuOut,*)' gbs_map_print: ERROR ' 2414 write(LuOut,*)' number of centers is zero or weird' 2415 write(LuOut,*)' nat = ',nat 2416 gbs_map_print = .false. 2417c..... add diagnostics later 2418 return 2419 endif 2420c 2421c print global information 2422c 2423 write(LuOut,*)'<<< GBS_MAP_PRINT >>>' 2424 write(LuOut,*)' total number of atoms :',nat 2425 write(LuOut,*)' total number of contractions :', 2426 & ncont_tot_gb(basis) 2427 write(LuOut,*)' total number of primitives :', 2428 & nprim_tot_gb(basis) 2429 write(LuOut,*)' total number of coefficients :', 2430 & ncoef_tot_gb(basis) 2431 write(LuOut,*)' total number of basis functions :', 2432 & nbf_tot_gb(basis) 2433c 2434c print center based mapping information. 2435c 2436 write(LuOut,*)' ' 2437 write(LuOut,*)'=================================================', 2438 & '===============================' 2439 write(LuOut,*)' center -> unique center map <ibs_ce2uce>' 2440 write(LuOut,*)' -> contraction range map <ibs_ce2cnr>' 2441 write(LuOut,*)'=================================================', 2442 & '===============================' 2443 do 00100 i=1,nat 2444 write(LuOut,'(1x,a,i4,2x,a,i3)') 2445 & 'center:',i,'maps to unique center:', 2446 & sf_ibs_ce2uce(i,basis) 2447 myfirst = sf_ibs_ce2cnr(1,i,basis) 2448 mylast = sf_ibs_ce2cnr(2,i,basis) 2449 mysize = mylast - myfirst + 1 2450 write(LuOut,'(14x,a,i4,2x,a,i4,a,i4,a,/)') 2451 & 'has',mysize,'contractions <first:', 2452 & myfirst,'> <last:',mylast,'>' 245300100 continue 2454c 2455c print contraction based mapping information 2456c 2457 write(LuOut,*)' ' 2458 write(LuOut,*)'=================================================', 2459 & '===============================' 2460 write(LuOut,*)' contraction -> center map ', 2461 & ' <ibs_cn2ce>' 2462 write(LuOut,*)' -> unique contraction in basis set', 2463 & ' <ibs_cn2ucn>' 2464 write(LuOut,*)' -> basis function range ', 2465 & ' <ibs_cn2bfr>' 2466 write(LuOut,*)'=================================================', 2467 & '===============================' 2468c 2469 do 00200 i=1,ncont_tot_gb(basis) 2470 mycenter = sf_ibs_cn2ce(i,basis) 2471 myucont = sf_ibs_cn2ucn(i,basis) 2472 myfirst = sf_ibs_cn2bfr(1,i,basis) 2473 mylast = sf_ibs_cn2bfr(2,i,basis) 2474 mysize = mylast - myfirst + 1 2475 write(LuOut,'(1x,a,i4,2x,a,i3)') 2476 & 'contraction',i,'is on center:',mycenter 2477 write(LuOut,'(18x,a,i3)') 2478 & 'is represented by unique contraction:',myucont 2479 write(LuOut,'(18x,a,i5,a,i5,a,i5,a,/)') 2480 & 'has',mysize,' basis functions <first:',myfirst, 2481 & '> <last:',mylast,'>' 248200200 continue 2483c 2484 gbs_map_print = .true. 2485 return 2486 end 2487*..................................................................... 2488C> \ingroup bas 2489C> @{ 2490c 2491C> \brief Store the names of all current basis set instances on the RTDB 2492c 2493C> Write the number of current basis set instances and the names 2494C> of all basis set instances to the RTDB. No actual contents of the 2495C> basis sets is stored by this function. The bas_rtdb_store routine 2496C> can be used to store the actual basis set data instead. 2497c 2498C> \return Return .true. if successfull, and .false. otherwise 2499c 2500 logical function bas_rtdb_out(rtdb) 2501c 2502c output to rtdb info about known basis sets 2503c 2504 implicit none 2505#include "mafdecls.fh" 2506#include "nwc_const.fh" 2507#include "basP.fh" 2508#include "rtdb.fh" 2509#include "stdio.fh" 2510c 2511c::passed 2512 integer rtdb !< [Input] the RTDB handle 2513c 2514 bas_rtdb_out = 2515 & rtdb_put(rtdb, 'basis:nbasis', 2516 & MT_INT, 1, nbasis_rtdb) 2517 & .and. 2518 & rtdb_cput(rtdb, 'basis:names', nbasis_rtdb, 2519 & bs_names_rtdb) 2520 if (.not. bas_rtdb_out) 2521 & write(LuOut,*) ' bas_rtdb_out: rtdb is corrupt ' 2522c 2523 end 2524*..................................................................... 2525c 2526C> \brief Add another basis set name to the list of known basis sets on 2527C> the RTDB 2528c 2529C> Just add another name to the list of known basis set instances. No 2530C> actual basis set data is stored. 2531c 2532C> \return Return .true. if successfull, and .false. otherwise 2533c 2534 logical function bas_rtdb_add(rtdb, name) 2535 implicit none 2536c 2537c add basis set name to known basis set list on rtdb 2538c 2539#include "mafdecls.fh" 2540#include "nwc_const.fh" 2541#include "basP.fh" 2542#include "rtdb.fh" 2543#include "inp.fh" 2544#include "stdio.fh" 2545c 2546 integer rtdb !< [Input] the RTDB handle 2547 character*(*) name !< [Input] the name of basis set to add 2548 integer basis 2549 logical status 2550 integer ln 2551 logical bas_rtdb_in, bas_rtdb_out 2552 external bas_rtdb_in, bas_rtdb_out 2553c 2554c See if name is on the rtdb already 2555c 2556 ln = inp_strlen(name) 2557 status = bas_rtdb_in(rtdb) 2558 bas_rtdb_add = .true. 2559 do 00100 basis = 1, nbasis_rtdb 2560 if (name(1:ln) .eq. 2561 & bs_names_rtdb(basis)(1:len_bs_rtdb(basis))) return 256200100 continue 2563c 2564c Name is not present ... add and rewrite info 2565c 2566 if (nbasis_rtdb .eq. nbasis_rtdb_mx) then 2567 write(LuOut,*) ' bas_rtdb_add: too many basis tries on rtdb ', 2568 & name 2569 bas_rtdb_add = .false. 2570 return 2571 endif 2572 nbasis_rtdb = nbasis_rtdb + 1 2573 bs_names_rtdb(nbasis_rtdb) = name 2574 len_bs_rtdb(nbasis_rtdb) = ln 2575c 2576 bas_rtdb_add = bas_rtdb_out(rtdb) 2577 if (.not. bas_rtdb_add) then 2578 write(LuOut,*) ' bas_rtdb_add: rtdb error adding ', name(1:ln) 2579 return 2580 endif 2581c 2582 bas_rtdb_add = .true. 2583c 2584 end 2585*..................................................................... 2586c 2587C> \brief Print contents of all current basis set instances 2588c 2589C> \return Return .true. if successful, and .false. otherwise 2590c 2591 logical function bas_print_all() 2592c 2593c routine to print active all basis set(s) information 2594c 2595 implicit none 2596#include "nwc_const.fh" 2597#include "basP.fh" 2598c::function 2599 logical bas_print 2600 external bas_print 2601c::local 2602 integer basis,basin 2603c 2604 bas_print_all = .true. 2605 do 00100 basis=1,nbasis_bsmx 2606 if(bsactive(basis)) then 2607 basin = basis - Basis_Handle_Offset 2608 bas_print_all = bas_print_all .and. bas_print(basin) 2609 endif 261000100 continue 2611c 2612 return 2613 end 2614c 2615C> \brief Print the names of all basis set instances on the RTDB 2616c 2617C> \return Return .true. if successful, and .false. otherwise 2618c 2619 subroutine bas_print_known_bases(rtdb) 2620 implicit none 2621#include "mafdecls.fh" 2622#include "nwc_const.fh" 2623#include "basP.fh" 2624#include "basdeclsP.fh" 2625#include "rtdb.fh" 2626#include "inp.fh" 2627#include "global.fh" 2628#include "stdio.fh" 2629 integer rtdb !< [Input] the RTDB handle 2630c 2631c Print list of known basis sets 2632c 2633 logical bas_rtdb_in 2634 external bas_rtdb_in 2635 logical ignore 2636 integer basis, ma_type, natom, nelem 2637 character*26 date 2638 character*32 name32 2639 character*128 key 2640 integer head_array(ndbs_head) 2641c 2642 ignore = bas_rtdb_in(rtdb) 2643 if (ga_nodeid() .eq. 0) then 2644 write(LuOut,*) 2645 call util_print_centered(LuOut,'Basis sets in the database', 2646 $ 23,.true.) 2647 write(LuOut,*) 2648 if (nbasis_rtdb .le. 0) then 2649 write(LuOut,*) ' There are no basis sets in the database' 2650 write(LuOut,*) 2651 else 2652 if (nbasis_rtdb .gt. 0) write(LuOut,3) 2653 3 format( 2654 $ 1x,4x,2x,'Name',28x,2x,'Natoms',2x, 2655 $ 'Last Modified',/, 2656 $ 1x,5x,2x,32('-'),2x,6('-'),2x,24('-')) 2657 do basis = 1, nbasis_rtdb 2658 key = ' ' 2659 write(key,'(''basis:'',a,'':header'')') 2660 $ bs_names_rtdb(basis)(1:len_bs_rtdb(basis)) 2661 if (.not. rtdb_get(rtdb, key, mt_int, 2662 $ ndbs_head, head_array)) then 2663 write(LuOut,*) ' Warning: basis set ', basis, 2664 $ ' may be corrupt' 2665 natom = -1 2666 else 2667 natom = head_array(HEAD_NTAGS) 2668 endif 2669 if (.not. rtdb_get_info(rtdb, key, ma_type, 2670 $ nelem, date)) then 2671 write(LuOut,*) ' Warning: basis set ', basis, 2672 $ ' may be corrupt' 2673 date = 'unknown' 2674 endif 2675 name32 = bs_names_rtdb(basis)(1:len_bs_rtdb(basis)) 2676 write(LuOut,4) basis, name32, natom, date 2677 4 format(1x,i4,2x,a32,2x,i6,2x,a26) 2678 end do 2679 if (nbasis_rtdb .gt. 0) then 2680 if (.not. rtdb_cget(rtdb,'ao basis',1,key)) 2681 $ key = 'ao basis' 2682 write(LuOut,*) 2683 write(LuOut,5) key(1:inp_strlen(key)) 2684 5 format(2x,'The basis set named "',a, 2685 $ '" is the default AO basis for restart') 2686 endif 2687 write(LuOut,*) 2688 write(LuOut,*) 2689 endif 2690 call util_flush(LuOut) 2691 endif 2692c 2693 end 2694C> @} 2695*..................................................................... 2696 logical function bas_rtdb_in(rtdb) 2697c 2698c load in info about known basis sets ... this is more 2699c for diagnostic and debugging purposes 2700c 2701 implicit none 2702#include "mafdecls.fh" 2703#include "nwc_const.fh" 2704#include "basP.fh" 2705#include "rtdb.fh" 2706#include "inp.fh" 2707#include "stdio.fh" 2708c 2709 integer rtdb ! [input] run time data base handle 2710 integer bas 2711 bas_rtdb_in = .false. 2712 nbasis_rtdb = 0 2713 if (rtdb_get(rtdb, 'basis:nbasis', MT_INT, 1, nbasis_rtdb)) 2714 $ then 2715 if (.not. rtdb_cget(rtdb,'basis:names', nbasis_rtdb_mx, 2716 $ bs_names_rtdb)) then 2717 write(LuOut,*) 'bas_rtdb_in: rtdb corrupt' 2718 else 2719 do 00100 bas = 1, nbasis_rtdb 2720 len_bs_rtdb(bas) = inp_strlen(bs_names_rtdb(bas)) 272100100 continue 2722 bas_rtdb_in = .true. 2723 endif 2724 endif 2725c 2726 return 2727 end 2728*..................................................................... 2729C> \ingroup bas 2730C> @{ 2731c 2732C> \brief Retrieve the center rank of a given shell in a basis set 2733c 2734C> After the basis set has been mapped to a geometry a complete list 2735C> of shells (or contractions) is available. This routine retrieves 2736C> the center rank of a shell from the complete list of shells. 2737c 2738C> \return Return .true. if successful, and .false. otherwise. 2739c 2740 logical function bas_cn2ce(basisin,cont,center) 2741 implicit none 2742c 2743c returns the center for a given mapped contraction 2744c 2745#include "mafdecls.fh" 2746#include "nwc_const.fh" 2747#include "basP.fh" 2748#include "geobasmapP.fh" 2749#include "bas_ibs_dec.fh" 2750#include "stdio.fh" 2751c::function 2752 logical bas_check_handle 2753 external bas_check_handle 2754c::passed 2755 integer basisin !< [Input] the basis set handle 2756 integer cont !< [Input] the mapped contraction index 2757 integer center !< [Output] the center rank 2758c::local 2759 integer basis 2760#include "bas_ibs_sfn.fh" 2761c 2762 bas_cn2ce = bas_check_handle(basisin,'bas_cn2ce') 2763 if(.not.bas_cn2ce) return 2764c 2765 basis = basisin + Basis_Handle_Offset 2766 bas_cn2ce = cont.ge.0 .and. cont.le.ncont_tot_gb(basis) 2767 if (.not.bas_cn2ce) then 2768 write(LuOut,*)' bas_cn2ce: invalid contraction information ' 2769 write(LuOut,*)' contraction range is 0:',ncont_tot_gb(basis) 2770 write(LuOut,*)' input contraction was: ',cont 2771 return 2772 endif 2773 center = sf_ibs_cn2ce(cont,basis) 2774c 2775 return 2776 end 2777*..................................................................... 2778c 2779C> \brief Retrieve the rank of the symmetry unique center of a shell in 2780C> a basis set 2781c 2782C> After the basis set has been mapped to a geometry a complete list 2783C> of shells (or contractions) is available. This routine retrieves 2784C> the center rank of the symmetry unique center of a shell from the 2785C> complete list of shells. 2786c 2787C> \return Return .true. if successful, and .false. otherwise. 2788c 2789 logical function bas_cn2uce(basisin,cont,ucenter) 2790 implicit none 2791c 2792c returns the UNIQUE center for a given mapped contraction 2793c 2794#include "mafdecls.fh" 2795#include "nwc_const.fh" 2796#include "basP.fh" 2797#include "geobasmapP.fh" 2798#include "bas_ibs_dec.fh" 2799#include "stdio.fh" 2800c::function 2801 logical bas_check_handle 2802 external bas_check_handle 2803c::passed 2804 integer basisin !< [Input] the basis set handle 2805 integer cont !< [Input] the mapped contraction index 2806 integer ucenter !< [Output] the symmetry unique center rank 2807c::local 2808 integer basis 2809#include "bas_ibs_sfn.fh" 2810c 2811 bas_cn2uce = bas_check_handle(basisin,'bas_cn2uce') 2812 if(.not.bas_cn2uce) return 2813c 2814 basis = basisin + Basis_Handle_Offset 2815 bas_cn2uce = cont.ge.0 .and. cont.le.ncont_tot_gb(basis) 2816 if (.not.bas_cn2uce) then 2817 write(LuOut,*)' bas_cn2uce: invalid contraction information ' 2818 write(LuOut,*)' contraction range is 0:',ncont_tot_gb(basis) 2819 write(LuOut,*)' input contraction was: ',cont 2820 return 2821 endif 2822c 2823 ucenter = sf_ibs_cn2ce(cont,basis) 2824 ucenter = sf_ibs_ce2uce(ucenter,basis) 2825c 2826 end 2827*..................................................................... 2828c 2829C> \brief Retrieve the range of basis functions corresponding to a shell 2830c 2831C> A shell (or equivalently contraction) may contain one or more basis 2832C> functions. Hence in the total molecular basis a given shell 2833C> corresponds to a range of basis functions. This routine retrieve 2834C> such a range of basis functions in that it 2835C> 2836C> - stores the first basis function of a shell in ifirst 2837C> 2838C> - stores the last basis function of a shell in ilast 2839C> 2840C> \return Return .true. if successfull, and .false. otherwise. 2841c 2842 logical function bas_cn2bfr(basisin,cont,ifirst,ilast) 2843c 2844c returns the first basis function index of a mapped contraction 2845c in ifirst and the last basis function index in ilast 2846c 2847 implicit none 2848#include "mafdecls.fh" 2849#include "nwc_const.fh" 2850#include "basP.fh" 2851#include "geobasmapP.fh" 2852#include "bas_ibs_dec.fh" 2853#include "stdio.fh" 2854c::function 2855 logical bas_check_handle 2856 external bas_check_handle 2857c::passed 2858 integer basisin !< [Input] the basis set handle 2859 integer cont !< [Input] the mapped contraction index 2860 integer ifirst !< [Output] the first basis function 2861 integer ilast !< [Output] the last basis function 2862c::local 2863 integer basis 2864c 2865#include "bas_ibs_sfn.fh" 2866c 2867 bas_cn2bfr = .true. 2868#ifdef BASIS_DEBUG 2869 bas_cn2bfr = bas_check_handle(basisin,'bas_cn2bfr') 2870 if(.not.bas_cn2bfr) return 2871#endif 2872c 2873 basis = basisin + Basis_Handle_Offset 2874 bas_cn2bfr = cont.ge.0 .and. cont.le.ncont_tot_gb(basis) 2875 if (.not.bas_cn2bfr) then 2876 write(LuOut,*)' bas_cn2bfr: invalid contraction information ' 2877 write(LuOut,*)' contraction range is 0:',ncont_tot_gb(basis) 2878 write(LuOut,*)' input contraction was: ',cont 2879 return 2880 endif 2881c 2882 ifirst = sf_ibs_cn2bfr(1,cont,basis) 2883 ilast = sf_ibs_cn2bfr(2,cont,basis) 2884c 2885 return 2886 end 2887*..................................................................... 2888c 2889C> \brief Retrieve the range of basis functions corresponding to a 2890C> center 2891c 2892C> A center may contain zero or more basis functions. Hence in the total 2893C> molecular basis a given center corresponds to a range of basis 2894C> functions. This routine retrieves such a range of basis functions in 2895C> that it 2896C> 2897C> - stores the first basis function of a center in ibflo 2898C> 2899C> - stores the last basis function of a center in ibfhi 2900C> 2901C> In case that a center does not have any basis functions at all 2902C> 2903C> - ibflo is set to 0 2904C> 2905C> - ibfhi is set to -1 2906C> 2907C> \return Return .true. if successfull, and .false. otherwise. 2908c 2909 logical function bas_ce2bfr(basis, icent, ibflo, ibfhi) 2910c 2911c returns the basis function range for a given center 2912c 2913 implicit none 2914#include "stdio.fh" 2915 integer basis !< [Input] the basis set handle 2916 integer icent !< [Input] the center rank 2917 integer ibflo !< [Output] the first basis function on the center 2918 integer ibfhi !< [Output] the last basis function on the center 2919c 2920 integer cnlo, cnhi, tmp 2921 logical status 2922 logical bas_ce2cnr, bas_cn2bfr 2923 external bas_ce2cnr, bas_cn2bfr 2924c 2925 status = .true. 2926 status = status .and. bas_ce2cnr(basis, icent, cnlo, cnhi) 2927 if (cnhi .gt. 0) then 2928 status = status .and. bas_cn2bfr(basis, cnlo, ibflo, tmp) 2929 status = status .and. bas_cn2bfr(basis, cnhi, tmp, ibfhi) 2930 else 2931 ibflo = 0 2932 ibfhi = -1 2933 endif 2934c 2935 bas_ce2bfr = status 2936c 2937 end 2938*..................................................................... 2939c 2940C> \brief Retrieves the range of shells for a given center 2941C> 2942C> A center may have zero or more shells (or equivalently contractions) 2943C> associated with it. This routine retrieves the range of shells 2944C> associated with a center in that it 2945C> 2946C> - stores the rank of the first shell in ifirst 2947C> 2948C> - stores the rank of the last shell in ilast 2949C> 2950C> In case center does not have any shells (e.g. a point charge) the 2951C> results are undefined. 2952C> 2953C> \return Return .true. if successfull, and .false. otherwise. 2954c 2955 logical function bas_ce2cnr(basisin,center,ifirst,ilast) 2956c 2957c returns the mapped contraction range on a given center 2958c 2959 implicit none 2960#include "mafdecls.fh" 2961#include "nwc_const.fh" 2962#include "basP.fh" 2963#include "geobasmapP.fh" 2964#include "geom.fh" 2965#include "bas_ibs_dec.fh" 2966#include "stdio.fh" 2967c::function 2968 logical bas_check_handle 2969 external bas_check_handle 2970c::passed 2971 integer basisin !< [Input] the basis set handle 2972 integer center !< [Input] the center rank 2973 integer ifirst !< [Output] the first mapped contraction 2974 integer ilast !< [Output] the last mapped contraction 2975c::local 2976 integer basis, nat 2977#include "bas_ibs_sfn.fh" 2978c 2979 bas_ce2cnr = .true. 2980#ifdef BASIS_DEBUG 2981 bas_ce2cnr = bas_check_handle(basisin,'bas_ce2cnr') 2982 if(.not.bas_ce2cnr) return 2983#endif 2984c 2985 basis = basisin + Basis_Handle_Offset 2986 bas_ce2cnr = geom_ncent(ibs_geom(basis),nat) 2987 if (nat.eq.0.or..not.bas_ce2cnr) then 2988 bas_ce2cnr = .false. 2989 write(LuOut,*)' bas_ce2cnr: ERROR ' 2990 write(LuOut,*)' number of centers is zero or weird' 2991 write(LuOut,*)' nat = ',nat 2992c..... add diagnostics later 2993 return 2994 endif 2995 2996 bas_ce2cnr = center.gt.0 .and. center.le.nat 2997 if (.not.bas_ce2cnr) then 2998 write(LuOut,*)' bas_ce2cnr: invalid center information ' 2999 write(LuOut,*)' center range is 1:',nat 3000 write(LuOut,*)' input center was : ',center 3001 return 3002 endif 3003c 3004 ifirst = sf_ibs_ce2cnr(1, center, basis) 3005 ilast = sf_ibs_ce2cnr(2, center, basis) 3006c 3007 return 3008 end 3009*..................................................................... 3010c 3011C> \brief Retrieves the center a given basis function is sited on 3012C> 3013C> Every Gaussian basis function is sited at one particular expansion 3014C> center. This routine routine retrieves the particular center on 3015C> which a given basis function is sited. 3016C> 3017C> \return Return .true. if successfull, and .false. otherwise. 3018c 3019 logical function bas_bf2ce(basisin,testbf,center) 3020c 3021c routine to return the center of a given basis function 3022c 3023 implicit none 3024#include "mafdecls.fh" 3025#include "nwc_const.fh" 3026#include "basP.fh" 3027#include "geobasmapP.fh" 3028#include "geom.fh" 3029#include "bas_ibs_dec.fh" 3030#include "stdio.fh" 3031c::function 3032 logical bas_check_handle 3033 external bas_check_handle 3034c::passed 3035 integer basisin !< [Input] the basis set handle 3036 integer testbf !< [Input] the basis function index 3037 integer center !< [Output] the center rank 3038c::local 3039 integer basis, nat, iat, ibflo, ibfhi, istart, iend 3040c 3041#include "bas_ibs_sfn.fh" 3042c 3043 bas_bf2ce = bas_check_handle(basisin,'bas_bf2ce') 3044 if (.not. bas_bf2ce) return 3045c 3046 basis = basisin + Basis_Handle_Offset 3047 bas_bf2ce = geom_ncent(ibs_geom(basis),nat) 3048 if (.not.bas_bf2ce .or. nat.le.0) then 3049 bas_bf2ce = .false. 3050 write(LuOut,*)' bas_bf2ce: ERROR ' 3051 write(LuOut,*)' number of centers is zero or weird' 3052 write(LuOut,*)' nat = ',nat 3053 return 3054 endif 3055c 3056c... linear search through atoms 3057c 3058 center = -1 3059 do 00100 iat = 1,nat 3060 istart = sf_ibs_ce2cnr(1,iat,basis) 3061 iend = sf_ibs_ce2cnr(2,iat,basis) 3062 if ((iend - istart + 1 ) .le. 0) goto 00100 3063 ibflo = sf_ibs_cn2bfr(1,istart,basis) 3064 ibfhi = sf_ibs_cn2bfr(2,iend,basis) 3065 if (testbf.ge.ibflo.and.testbf.le.ibfhi) then 3066 center = iat 3067 return 3068 endif 306900100 continue 3070c 3071 end 3072*..................................................................... 3073c 3074C> \brief Retrieves the shell a given basis function is part off 3075C> 3076C> Every Gaussian basis function is part of one particular shell 3077C> (or equivalently contraction). This routine routine retrieves the 3078C> particular shell a given basis function is part off. 3079C> 3080C> \return Return .true. if successfull, and .false. otherwise. 3081c 3082 logical function bas_bf2cn(basisin,testbf,cont) 3083c 3084c returns the mapped contraction index that contains the given 3085c basis function index 3086c 3087 implicit none 3088#include "mafdecls.fh" 3089#include "nwc_const.fh" 3090#include "basP.fh" 3091#include "geobasmapP.fh" 3092#include "bas_ibs_dec.fh" 3093#include "stdio.fh" 3094c::function 3095 logical bas_numcont 3096 logical bas_check_handle 3097 external bas_numcont 3098 external bas_check_handle 3099c::passed 3100 integer basisin !< [Input] the basis set handle 3101 integer testbf !< [Input] the basis function index 3102 integer cont !< [Output] the mapped shell index 3103c::local 3104 integer basis, icont, ibflo, ibfhi 3105 integer numcont 3106c 3107#include "bas_ibs_sfn.fh" 3108c 3109 bas_bf2cn = bas_check_handle(basisin,'bas_bf2cn') 3110 if (.not. bas_bf2cn) return 3111c 3112 bas_bf2cn = bas_numcont(basisin,numcont) 3113 if (.not.bas_bf2cn) then 3114 write(LuOut,*)'bas_bf2cn: could not get number of contractions' 3115 return 3116 endif 3117c 3118 basis = basisin + Basis_Handle_Offset 3119c 3120c... linear search through contractions 3121 cont = -1 3122 do 00100 icont = 1,numcont 3123 ibflo = sf_ibs_cn2bfr(1,icont,basis) 3124 ibfhi = sf_ibs_cn2bfr(2,icont,basis) 3125 if(testbf.ge.ibflo.and.testbf.le.ibfhi) then 3126 cont = icont 3127 return 3128 endif 312900100 continue 3130c 3131 end 3132*..................................................................... 3133c 3134C> \brief Retrieves the number of basis functions in the molecular 3135C> basis set 3136C> 3137C> The basis set itself contains only the unique definitions of the 3138C> combinations of primitive basis functions that generate the actual 3139C> basis functions. I.e. it contains only the contraction coefficients, 3140C> the exponents and the kind of atom it should be used for. 3141C> 3142C> The molecular basis set is constructed from this information by 3143C> mapping this specification onto a particular structure. This process 3144C> Generates the actual specification of basis functions for the atoms. 3145C> In molecular calculations various properties of the molecular basis 3146C> set are essential. 3147C> 3148C> This routine extracts to total number of basis functions in the 3149C> molecular basis set. 3150C> 3151C> \return Return .true. if successfull, and .false. otherwise. 3152c 3153 logical function bas_numbf(basisin,nbf) 3154c 3155c returns the total number of basis functions of the mapped basis set. 3156c 3157 implicit none 3158#include "nwc_const.fh" 3159#include "basP.fh" 3160#include "geobasmapP.fh" 3161c::function 3162 logical bas_check_handle 3163 external bas_check_handle 3164c::passed 3165 integer basisin !< [Input] the basis set handle 3166 integer nbf !< [Output] the number of basis functions 3167c::local 3168 integer basis 3169c 3170 nbf = -6589 3171 bas_numbf = bas_check_handle(basisin,'bas_numbf') 3172 if (.not. bas_numbf) return 3173 3174 basis = basisin + Basis_Handle_Offset 3175 nbf = nbf_tot_gb(basis) 3176 bas_numbf = .true. 3177 return 3178 end 3179*..................................................................... 3180c 3181C> \brief Retrieve the total number of primitive basis functions in the 3182C> molecular basis set 3183C> 3184C> A shell of basis functions consists of a number of contracted 3185C> Gaussians of a particular angular momentum. The total number of basis 3186C> functions corresponds to the number of contracted functions. Each 3187C> individual Gaussian function in a contraction is referred to as 3188C> a primitive function. Clearly then there is a number of primitive 3189C> basis functions associated with the molecular basis set as well. 3190C> This routine retrieves the latter number of primitive basis functions 3191C> from the basis set. 3192C> 3193C> \return Returns .true. if successfull, and .false. otherwise. 3194c 3195 logical function bas_num_prim(basisin,nprim) 3196c 3197c returns the total number of primitives of the mapped basis set. 3198c 3199 implicit none 3200#include "nwc_const.fh" 3201#include "basP.fh" 3202#include "geobasmapP.fh" 3203c::function 3204 logical bas_check_handle 3205 external bas_check_handle 3206c::passed 3207 integer basisin !< [Input] the basis set handle 3208 integer nprim !< [Output] the number of primitives in mapped basis 3209c::local 3210 integer basis 3211c 3212 nprim = -6589 3213 bas_num_prim = bas_check_handle(basisin,'bas_num_prim') 3214 if (.not. bas_num_prim) return 3215 3216 basis = basisin + Basis_Handle_Offset 3217 nprim = nprim_tot_gb(basis) 3218 bas_num_prim = .true. 3219 return 3220 end 3221*..................................................................... 3222c 3223C> \brief Retrieve the total number of contraction coefficients of the 3224C> molecular basis set 3225C> 3226C> Shells consist of contracted Gaussian basis functions. With each 3227C> primitive Gaussian there is associated an exponent and at least 3228C> one contraction coefficient. 3229C> 3230C> This routine retrieves the total number of contraction coefficients 3231C> in the molecular basis set. 3232C> 3233C> \return Returns .true. if successfull, and .false. otherwise. 3234c 3235 logical function bas_num_coef(basisin,ncoef) 3236c 3237c returns the total number of coeffsof the mapped basis set. 3238c 3239 implicit none 3240#include "nwc_const.fh" 3241#include "basP.fh" 3242#include "geobasmapP.fh" 3243c::function 3244 logical bas_check_handle 3245 external bas_check_handle 3246c::passed 3247 integer basisin !< [Input] the basis set handle 3248 integer ncoef !< [Output] the number of contraction coefficients 3249 !< in the mapped basis 3250c::local 3251 integer basis 3252c 3253 ncoef = -6589 3254 bas_num_coef = bas_check_handle(basisin,'bas_num_coef') 3255 if (.not. bas_num_coef) return 3256 3257 basis = basisin + Basis_Handle_Offset 3258 ncoef = ncoef_tot_gb(basis) 3259 bas_num_coef = .true. 3260 return 3261 end 3262*..................................................................... 3263c 3264C> \brief Retrieve the names of a basis set instance 3265C> 3266C> Each basis set instance has two names. By one name the basis set 3267C> in memory is known. This name is defined at the time the basis set 3268C> instance is created. The other name is the one used to store the 3269C> basis set on the RTDB. 3270C> 3271C> This routine retrieves both basis set names. 3272C> 3273C> \return Returns .true. if successfull, and .false. otherwise. 3274c 3275 logical function bas_name(basisin,basis_name,trans_name) 3276c 3277c returns the name and translated name of the basis set 3278c 3279 implicit none 3280#include "nwc_const.fh" 3281#include "basP.fh" 3282#include "inp.fh" 3283c::functions 3284c inp_strlen() from inp 3285 logical bas_check_handle 3286 external bas_check_handle 3287c::passed 3288 integer basisin !< [Input] the basis set handle 3289 character*(*) basis_name !< [Output] the basis set name when loaded 3290 character*(*) trans_name !< [Output] the basis set name in context 3291c::local 3292 integer basis ! actual offset into basis arrays 3293* integer lenofit ! length of name 3294c 3295 bas_name = bas_check_handle(basisin,'bas_name') 3296 if (.not. bas_name) return 3297c 3298 basis = basisin + Basis_Handle_Offset 3299c 3300* lenofit = inp_strlen(bs_name(basis)) 3301* basis_name(1:lenofit) = bs_name(basis)(1:lenofit) 3302* lenofit = inp_strlen(bs_trans(basis)) 3303* basis_name(1:lenofit) = bs_trans(basis)(1:lenofit) 3304 basis_name = bs_name(basis) 3305 trans_name = bs_trans(basis) 3306c 3307 end 3308*..................................................................... 3309c 3310C> \brief Retrieves the tag of a given shell 3311C> 3312C> \return Returns .true. if successfull, and .false. otherwise. 3313c 3314 logical function bas_cont_tag(basisin,icont,tagout) 3315 implicit none 3316#include "nwc_const.fh" 3317#include "basP.fh" 3318#include "stdio.fh" 3319#include "mafdecls.fh" 3320#include "geobasmapP.fh" 3321#include "bas_ibs_dec.fh" 3322c::-function 3323 logical bas_check_handle 3324 external bas_check_handle 3325c::-passed 3326 integer basisin !< [Input] the basis set handle 3327 integer icont !< [Input] the shell index 3328 character*(*) tagout !< [Output] the shell tag 3329c::-local 3330 integer center, ucenter, basis 3331 integer len_tagout 3332#include "bas_ibs_sfn.fh" 3333c 3334 bas_cont_tag = bas_check_handle(basisin,'bas_cont_tag') 3335 if (.not.bas_cont_tag) return 3336 3337 basis = basisin + Basis_Handle_Offset 3338 3339 center = sf_ibs_cn2ce(icont,basis) 3340 ucenter = sf_ibs_ce2uce(center,basis) 3341 len_tagout = len(tagout) 3342 tagout = bs_tags(ucenter,basis)(1:len_tagout) 3343 end 3344*..................................................................... 3345c 3346C> \brief Retrieves generic information of a given shell 3347C> 3348C> \return Returns .true. if successfull, and .false. otherwise. 3349c 3350 logical function bas_continfo(basisin,icont, 3351 & type,nprimo,ngeno,sphcart) 3352c 3353c returns the generic information about the given mapped contraction 3354c 3355 implicit none 3356#include "mafdecls.fh" 3357#include "nwc_const.fh" 3358#include "basP.fh" 3359#include "geobasmapP.fh" 3360#include "basdeclsP.fh" 3361#include "bas_ibs_dec.fh" 3362#include "stdio.fh" 3363c::function 3364 logical bas_check_handle 3365 external bas_check_handle 3366c::passed 3367 integer basisin !< [Input] the basis handle 3368 integer icont !< [Input] the shell index 3369 integer type !< [Output] the shell type (sp/s/p/d/..) 3370 integer nprimo !< [Output] the number of primitives 3371 integer ngeno !< [Output] the number of contractions 3372 integer sphcart !< [Output] 0/1 for cartesian/spherical harmonic 3373 !< basis functions 3374c::local 3375 integer basis,myucont,icontmax 3376c 3377#include "bas_ibs_sfn.fh" 3378c 3379 nprimo = -123 3380 ngeno = -456 3381 sphcart = -789 3382c 3383 bas_continfo = bas_check_handle(basisin,'bas_continfo') 3384 if (.not.bas_continfo) return 3385 3386 basis = basisin + Basis_Handle_Offset 3387c 3388 icontmax = ncont_tot_gb(basis) 3389c 3390 if (.not.(icont.ge.0.and.icont.le.icontmax)) then 3391 write(LuOut,*)' bas_continfo: ERROR ' 3392 write(LuOut,*)' contraction range for basis is 0:', 3393 & icontmax 3394 write(LuOut,*)' information requested for contraction:',icont 3395 bas_continfo = .false. 3396 return 3397 endif 3398c 3399 myucont = sf_ibs_cn2ucn(icont,basis) 3400c... 3401 if (bas_spherical(basis)) then 3402 sphcart = 1 3403 else 3404 sphcart = 0 3405 endif 3406 type = infbs_cont(CONT_TYPE, myucont,basis) 3407 nprimo = infbs_cont(CONT_NPRIM,myucont,basis) 3408 ngeno = infbs_cont(CONT_NGEN, myucont,basis) 3409 bas_continfo=.true. 3410 return 3411 end 3412*..................................................................... 3413c 3414C> \brief Retrieves the total number of shells in the molecular 3415C> basis set 3416C> 3417C> \return Returns .true. if successfull, and .false. otherwise. 3418c 3419 logical function bas_numcont(basisin,numcont) 3420c 3421c returns the total number of mapped contractions/shells for the 3422c given basis set 3423c 3424 implicit none 3425#include "nwc_const.fh" 3426#include "basP.fh" 3427#include "geobasmapP.fh" 3428#include "basdeclsP.fh" 3429c::function 3430 logical bas_check_handle 3431 external bas_check_handle 3432c::passed 3433 integer basisin !< [Input] the basis set handle 3434 integer numcont !< [Output] the number of mapped shells 3435c::local 3436 integer basis 3437c 3438 numcont = -6589 3439 bas_numcont = bas_check_handle(basisin,'bas_numcont') 3440 if (.not.bas_numcont) return 3441 3442 basis = basisin + Basis_Handle_Offset 3443 3444 numcont = ncont_tot_gb(basis) 3445 3446 bas_numcont = .true. 3447 return 3448 end 3449*..................................................................... 3450c 3451C> \brief Retrieves the maximum number of basis functions in any shell 3452C> 3453C> This function scans through all the shells in a basis set and finds 3454C> the maximum number of basis functions. It deals with Cartesian and 3455C> spherical harmonic basis functions as well as segmented and 3456C> generally contracted basis sets. 3457C> 3458C> \return Returns .true. if successful, and .false. otherwise 3459c 3460 logical function bas_nbf_cn_max(basisin,nbf_max) 3461 implicit none 3462c 3463c calculate, return and store maximum basis function block size 3464c for all contractions in a given basis. 3465c 3466#include "nwc_const.fh" 3467#include "basP.fh" 3468#include "basdeclsP.fh" 3469#include "stdio.fh" 3470c::functions 3471 logical bas_check_handle 3472 external bas_check_handle 3473 integer nbf_from_ucont 3474 external nbf_from_ucont 3475c:: passed 3476 integer basisin !< [Input] the basis set handle 3477 integer nbf_max !< [Output] the maximum number of basis functions 3478 !< in any shell 3479c:local 3480 integer basis, myucont, i, mynbf 3481c 3482 bas_nbf_cn_max = bas_check_handle(basisin,'bas_nbf_cn_max') 3483 if (.not. bas_nbf_cn_max ) then 3484 write(LuOut,*) 'bas_nbf_cn_max: basis handle not valid ' 3485 return 3486 endif 3487c 3488 basis = basisin + Basis_Handle_Offset 3489 if (nbfmax_bs(basis) .gt. -565) then 3490 nbf_max = nbfmax_bs(basis) 3491 bas_nbf_cn_max = .true. 3492 return 3493 endif 3494c 3495 myucont = infbs_head(HEAD_NCONT,basis) 3496 nbf_max = -565 3497c 3498 do 00100 i = 1,myucont 3499 mynbf = nbf_from_ucont(i,basisin) 3500 nbf_max = max(nbf_max, mynbf) 350100100 continue 3502c 3503 nbfmax_bs(basis) = nbf_max 3504 bas_nbf_cn_max = .true. 3505 return 3506 end 3507*..................................................................... 3508C> 3509C> \brief Retrieves the maximum number of basis functions on any center 3510C> 3511C> \return Return .true. if successfull, and .false. otherwise 3512c 3513 logical function bas_nbf_ce_max(basisin,nbf_max) 3514 implicit none 3515c 3516c calculate, return and store maximum basis function block size 3517c for all contractions in a given basis. 3518c 3519#include "mafdecls.fh" 3520#include "nwc_const.fh" 3521#include "basP.fh" 3522#include "geom.fh" 3523#include "basdeclsP.fh" 3524#include "geobasmapP.fh" 3525#include "bas_ibs_dec.fh" 3526#include "stdio.fh" 3527c::functions 3528 logical bas_check_handle 3529 external bas_check_handle 3530c:: passed 3531 integer basisin !< [Input] the basis set handle 3532 integer nbf_max !< [Output] the maximum number of basis functions 3533 !< on any center 3534c:local 3535 integer basis, mynat, iat, mylo, myhi, mynbf 3536 integer myclo, mychi, mynumcont 3537c 3538#include "bas_ibs_sfn.fh" 3539c 3540 bas_nbf_ce_max = bas_check_handle(basisin,'bas_nbf_ce_max') 3541 if (.not. bas_nbf_ce_max ) then 3542 write(LuOut,*) 'bas_nbf_ce_max: basis handle not valid ' 3543 return 3544 endif 3545c 3546 basis = basisin + Basis_Handle_Offset 3547c 3548 bas_nbf_ce_max = geom_ncent(ibs_geom(basis),mynat) 3549 if (mynat.le.0.or. .not.bas_nbf_ce_max) then 3550 write(LuOut,*)' bas_nbf_ce_max: ERROR ' 3551 write(LuOut,*)' number of centers is zero or weird' 3552 write(LuOut,*)' nat = ',mynat 3553 return 3554 endif 3555c 3556 nbf_max = -1 3557 do 00100 iat = 1,mynat 3558 myclo = (sf_ibs_ce2cnr(1,iat,basis)) 3559 mychi = (sf_ibs_ce2cnr(2,iat,basis)) 3560 mynumcont = mychi - myclo + 1 3561* Bqs with no basis functions have mynumcont < 1 3562 if (mynumcont.gt.0) then 3563 mylo = sf_ibs_cn2bfr(1,myclo,basis) 3564 myhi = sf_ibs_cn2bfr(2,mychi,basis) 3565 mynbf = myhi - mylo + 1 3566 nbf_max = max(nbf_max, mynbf) 3567 endif 356800100 continue 3569c 3570 return 3571 end 3572*..................................................................... 3573c 3574C> \brief Retrieve the geometry that was specified when loading the 3575C> basis 3576C> 3577C> Simply look up the geometry that was stored with the basis set when 3578C> the basis set was loaded from the RTDB. 3579C> 3580C> \return Return .true. if successfull, and .false. otherwise. 3581c 3582 logical function bas_geom(basisin,geom) 3583 implicit none 3584#include "errquit.fh" 3585#include "nwc_const.fh" 3586#include "basP.fh" 3587#include "basdeclsP.fh" 3588#include "geobasmapP.fh" 3589c::functions 3590 logical bas_check_handle 3591 external bas_check_handle 3592c:: passed 3593 integer basisin !< [Input] the basis set handle 3594 integer basis 3595 integer geom !< [Output] the geometry used to load basis set 3596c 3597 if (.not.bas_check_handle(basisin,'bas_geom')) 3598 & call errquit('bas_geom: handle invalid',911, BASIS_ERR) 3599c 3600 basis = basisin + Basis_Handle_Offset 3601c 3602 geom = ibs_geom(basis) 3603 bas_geom = .true. 3604 end 3605*............................................... 3606c 3607C> \brief Retrieve the maximum number of contraction coefficients 3608C> associated with any shell 3609C> 3610C> This function scans through all shells in a basis set and checks 3611C> the number of contraction coefficients for each. It takes into 3612C> account segmented and generally contracted shells. The maximum 3613C> number of contraction coefficients for any shell is passed back to 3614C> the calling routine. 3615C> 3616C> \return Return .true. if successfull, .false. otherwise. 3617c 3618 logical function bas_ncoef_cn_max(basisin, ncoef_max) 3619 implicit none 3620#include "errquit.fh" 3621#include "nwc_const.fh" 3622#include "basP.fh" 3623#include "basdeclsP.fh" 3624 integer basisin !< [Input] the basis set handle 3625 integer ncoef_max !< [Output] the maximum number of contraction 3626 !< coefficients 3627c 3628 integer basis 3629 integer nu_cont, iu_cont, my_prim, my_gen 3630c 3631 basis = basisin + Basis_Handle_Offset 3632 if (.not.bsactive(basis)) call errquit 3633 & ('bas_ncoef_cn_max: basis handle invalid',911, BASIS_ERR) 3634 nu_cont = infbs_head(HEAD_NCONT,basis) 3635 ncoef_max = 0 3636 do 00100 iu_cont = 1,nu_cont 3637 my_prim = infbs_cont(CONT_NPRIM,iu_cont,basis) 3638 my_gen = infbs_cont(CONT_NGEN, iu_cont,basis) 3639 ncoef_max = max(ncoef_max,(my_prim*my_gen)) 364000100 continue 3641 bas_ncoef_cn_max = .true. 3642 end 3643*..................................................................... 3644c 3645C> \brief Retrieve the maximum number of primitive basis functions 3646C> associated with any shell 3647C> 3648C> This function scans through all shells in a basis set and checks 3649C> the number of primitive basis functions (or equivalently the 3650C> number of exponents) for each. The maximum number of primitive basis 3651C> basis functions for any shell is passed back to the calling routine. 3652C> 3653C> \return Return .true. if successfull, .false. otherwise. 3654c 3655 logical function bas_nprim_cn_max(basisin, nprim_max) 3656 implicit none 3657#include "errquit.fh" 3658#include "nwc_const.fh" 3659#include "basP.fh" 3660#include "basdeclsP.fh" 3661 integer basisin !< [Input] the basis set handle 3662 integer nprim_max !< [Output] the maximum number of primitives in 3663 !< any shell 3664c 3665 integer basis 3666 integer nu_cont, iu_cont, my_prim 3667c 3668 basis = basisin + Basis_Handle_Offset 3669 if (.not.bsactive(basis)) call errquit 3670 & ('bas_nprim_cn_max: basis handle invalid',911, BASIS_ERR) 3671 nu_cont = infbs_head(HEAD_NCONT,basis) 3672 nprim_max = 0 3673 do 00100 iu_cont = 1,nu_cont 3674 my_prim = infbs_cont(CONT_NPRIM,iu_cont,basis) 3675 nprim_max = max(nprim_max,my_prim) 367600100 continue 3677 bas_nprim_cn_max = .true. 3678 end 3679C> @} 3680*..................................................................... 3681 logical function bas_norm_get(basisin,norm_id) 3682 implicit none 3683#include "nwc_const.fh" 3684#include "basP.fh" 3685#include "basdeclsP.fh" 3686c 3687 integer basisin ! [input] basis set handle 3688 integer norm_id ! [output] Normalization id type 3689c 3690 integer basis 3691c 3692 basis = basisin + Basis_Handle_Offset 3693 norm_id = bas_norm_id(basis) 3694 bas_norm_get = norm_id.ge.BasNorm_lo.and.norm_id.le.BasNorm_hi 3695 end 3696*..................................................................... 3697 logical function bas_norm_set(basisin,norm_id) 3698 implicit none 3699#include "nwc_const.fh" 3700#include "basP.fh" 3701#include "basdeclsP.fh" 3702c 3703 integer basisin ! [input] basis set handle 3704 integer norm_id ! [input] Normalization id type 3705c 3706 integer basis 3707c 3708 basis = basisin + Basis_Handle_Offset 3709c 3710 if (norm_id.ge.BasNorm_lo.and.norm_id.le.BasNorm_hi) then 3711 bas_norm_id(basis) = norm_id 3712 bas_norm_set = .true. 3713 else 3714 bas_norm_set = .false. 3715 endif 3716 end 3717*..................................................................... 3718 logical function bas_norm_print(basisin) 3719 implicit none 3720#include "nwc_const.fh" 3721#include "basP.fh" 3722#include "basdeclsP.fh" 3723#include "stdio.fh" 3724c 3725 integer basisin ! [input] basis set handle 3726c 3727 integer basis, norm_id 3728c 3729 basis = basisin + Basis_Handle_Offset 3730 norm_id = bas_norm_id(basis) 3731 bas_norm_print = .true. 3732 if (norm_id.eq.BasNorm_UN) then 3733 write(luout,*)' basis is unnormalized' 3734 else if (norm_id.eq.BasNorm_STD) then 3735 write(luout,*)' basis has standard normalization via int_norm' 3736 else if (norm_id.eq.BasNorm_2c) then 3737 write(luout,*) 3738 & ' basis has dft/fitting normalization via int_norm_2c' 3739 else if (norm_id.eq.BasNorm_rel) then 3740 write(luout,*) 3741 & ' basis has relativistic normalization via int_norm' 3742 else 3743 write(luout,*)' basis handle: ',basisin 3744 write(luout,*)' norm_id = ',norm_id 3745 write(luout,*)' does not match a known normalization mode' 3746 bas_norm_print = .false. 3747 endif 3748 end 3749*..................................................................... 3750C> \ingroup bas 3751C> @{ 3752C> 3753C> \brief Print unique in core ECP information 3754C> 3755C> \return Return .true. if successfull, and .false. otherwise. 3756c 3757 logical function ecp_print(ecpidin) 3758c 3759c routine to print unique ecpid information that is in core 3760c 3761 implicit none 3762#include "mafdecls.fh" 3763#include "nwc_const.fh" 3764#include "basP.fh" 3765#include "basdeclsP.fh" 3766#include "inp.fh" 3767#include "geom.fh" 3768#include "stdio.fh" 3769#include "bas_starP.fh" 3770#include "errquit.fh" 3771c 3772c function declarations 3773c 3774 logical ecp_check_handle 3775 external ecp_check_handle 3776c 3777cc AJL/Begin/SPIN-POLARISED ECPs 3778 logical ecp_get_high_chan 3779 external ecp_get_high_chan 3780cc AJL/End 3781c 3782c:: passed 3783 integer ecpidin !< [Input] the basis set handle 3784c:: local 3785 integer mytags, myucont, myprim, mycoef, ecpid 3786 integer i,j,k,l, ifcont, mygen, mytype, irexptr, iexptr, icfptr 3787 integer atn, len_tag, len_ele 3788c 3789cc AJL/Begin/SPIN-POLARISED ECPs 3790 integer mychannel 3791 integer channels 3792 character*5 channel_type(0:2) 3793cc AJL/End 3794c 3795 character*2 symbol 3796 character*16 element 3797 character*9 cartorsph 3798 character*3 ctype(-1:6) 3799 character*3 shell_type 3800 character*60 buffer 3801c 3802#include "bas_exndcf.fh" 3803c 3804 ctype(-1)='U L' 3805 ctype(0) ='U-s' 3806 ctype(1) ='U-p' 3807 ctype(2) ='U-d' 3808 ctype(3) ='U-f' 3809 ctype(4) ='U-g' 3810 ctype(5) ='U-h' 3811 ctype(6) ='U-i' 3812c 3813cc AJL/Begin/SPIN-POLARISED ECPs 3814cc Aesthetic clean up for print. Is this necessary? 3815C Can I just set the array without the if statement? To test. 3816 channels = 1 3817 if (.not.ecp_get_high_chan(ecpidin,channels)) 3818 & call errquit('ecp_print error',911, BASIS_ERR) 3819 3820 channel_type(0) = 'Both' 3821 channel_type(1) = 'Alpha' 3822 channel_type(2) = 'Beta' 3823cc AJL/End 3824c 3825 ecp_print = .true. 3826 ecpid = ecpidin + Basis_Handle_Offset 3827c 3828 ecp_print = ecp_check_handle(ecpidin,'ecp_print') 3829 if (.not. ecp_print) return 3830c 3831c print basis set information 3832c 3833 if (bas_spherical(ecpid)) then 3834 cartorsph = 'spherical' 3835 else 3836 cartorsph = 'cartesian' 3837 endif 3838 write(LuOut,1)bs_name(ecpid)(1:inp_strlen(bs_name(ecpid))), 3839 $ bs_trans(ecpid)(1:inp_strlen(bs_trans(ecpid))),cartorsph 3840 1 format(' ECP "',a,'" -> "',a,'"',' (',a,')'/ 3841 $ ' -----') 3842 mytags = infbs_head(HEAD_NTAGS,ecpid) 3843 if (mytags.le.0) then 3844 write(LuOut,*)'No explicit ECP functions are defined !' 3845 write(LuOut,*) 3846c 3847c there could be star tags defined, so check that before returning 3848c 3849 goto 00010 3850 endif 3851c 3852 myucont = infbs_head(HEAD_NCONT,ecpid) 3853 myprim = infbs_head(HEAD_NPRIM,ecpid) 3854 mycoef = infbs_head(HEAD_NCOEF,ecpid) 3855c 3856 do 00100 i=1,mytags 3857 3858 if (geom_tag_to_element(bs_tags(i,ecpid), symbol, element, 3859 $ atn)) then 3860 len_tag = inp_strlen(bs_tags(i,ecpid)) 3861 len_ele = inp_strlen(element) 3862 write (buffer, 3863 & '(a,'' ('',a,'') Replaces '',i5,'' electrons'' )') 3864 & bs_tags(i,ecpid)(1:len_tag), element(1:len_ele), 3865 & infbs_tags(Tag_Nelec,i,ecpid) 3866 else 3867 buffer = bs_tags(i,ecpid) 3868 endif 3869 len_tag = inp_strlen(buffer) 3870 call util_print_centered(LuOut, buffer, len_tag/2 + 1, .true.) 3871 3872 myucont = infbs_tags(TAG_NCONT,i,ecpid) 3873c 3874 ifcont = infbs_tags(TAG_FCONT,i,ecpid) 3875c 3876 write(LuOut,6) 3877cc AJL/Begin/SPIN-POLARISED ECPs 3878c 6 format( 3879c $ ' R-exponent Exponent Coefficients '/ 3880 6 format(13(' '), 3881 $ 'Channel R-exponent Exponent Coefficients'/ 3882cc AJL/End 3883 $ ' ------------ ',57('-')) 3884 do 00200 j=1,myucont 3885 myprim = infbs_cont(CONT_NPRIM,ifcont,ecpid) 3886 mygen = infbs_cont(CONT_NGEN,ifcont,ecpid) 3887 3888 mytype = infbs_cont(CONT_TYPE, ifcont, ecpid) 3889 shell_type = ctype(mytype) 3890 iexptr = infbs_cont(CONT_IEXP, ifcont,ecpid) - 1 3891 icfptr = infbs_cont(CONT_ICFP, ifcont,ecpid) - 1 3892 irexptr = infbs_cont(CONT_IREXP,ifcont,ecpid) - 1 3893 mychannel = infbs_cont(CONT_CHANNEL,ifcont,ecpid) 3894c 3895 do 00300 k=1,myprim 3896 write(LuOut,7) j, shell_type, 3897 & channel_type(mychannel), 3898 & sf_exndcf((irexptr+k),ecpid), 3899 & sf_exndcf((iexptr+k),ecpid), 3900 & (sf_exndcf((icfptr+k+(l-1)*myprim),ecpid),l=1,mygen) 390100300 continue 3902 write(LuOut,*) 3903 ifcont = ifcont + 1 390400200 continue 390500100 continue 3906c 3907c Check if we have star tag definitions in the basis set 3908c 390900010 if (star_nr_tags .gt. 0) then 3910 write(LuOut,*) 3911 write(LuOut,8) 3912 8 format(' In addition, one or more string tags have been', 3913 & ' defined containing a * .'/' These tags, and ', 3914 & 'their exceptions list are printed below.'// 3915 & ' Tag ',12(' '),' Description ',18(' '), 3916 & ' Excluding '/' ',16('-'),' ',30('-'),' ',30('-')) 3917 do i=1,star_nr_tags 3918 k = 1 3919 if (i .gt. 1) k = star_nr_excpt(i-1) + 1 3920 write(LuOut,9) star_tag(i), star_bas_typ(i)(1:30), 3921 & (star_excpt(j)(1:inp_strlen(star_excpt(j))), 3922 & j=k,star_nr_excpt(i)) 3923 9 format(1x,a16,1x,a30,1x,10(a)) 3924 enddo 3925 write(LuOut,*) 3926 write(LuOut,*) 3927 endif 3928c 3929c If geom is set print out the info about total basis info 3930c associated with the geometry also 3931c 3932c ... not done yet 3933c 3934 return 3935cc AJL/Begin/SPIN-POLARISED ECPs 3936c 7 format(1x,i2,1x,a3,1x,f9.2,2x,f14.6,20f15.6) 3937 7 format(1x,i2,1x,a3,7x,a5,3x,f9.2,2x,f14.6,20f15.6) 3938cc AJL/End 3939 end 3940*..................................................................... 3941C> 3942C> \brief Print unique in core spin-orbit potential information 3943C> 3944C> \return Return .true. if successfull, and .false. otherwise. 3945c 3946 logical function so_print(soidin) 3947c 3948c routine to print unique soid information that is in core 3949c 3950 implicit none 3951#include "mafdecls.fh" 3952#include "nwc_const.fh" 3953#include "basP.fh" 3954#include "basdeclsP.fh" 3955#include "inp.fh" 3956#include "geom.fh" 3957#include "stdio.fh" 3958#include "bas_starP.fh" 3959c 3960c function declarations 3961c 3962 logical so_check_handle 3963 external so_check_handle 3964c:: passed 3965 integer soidin !< [Input] the basis set handle 3966c:: local 3967 integer mytags, myucont, myprim, mycoef, soid 3968 integer i,j,k,l, ifcont, mygen, mytype, irexptr, iexptr, icfptr 3969 integer atn, len_tag, len_ele 3970 character*2 symbol 3971 character*16 element 3972 character*3 ctype(-1:6) 3973 character*3 shell_type 3974 character*60 buffer 3975c 3976#include "bas_exndcf.fh" 3977c 3978 ctype(-1)='U L' 3979 ctype(0) ='U-s' 3980 ctype(1) ='U-p' 3981 ctype(2) ='U-d' 3982 ctype(3) ='U-f' 3983 ctype(4) ='U-g' 3984 ctype(5) ='U-h' 3985 ctype(6) ='U-i' 3986 so_print = .true. 3987 soid = soidin + Basis_Handle_Offset 3988c 3989 so_print = so_check_handle(soidin,'so_print') 3990 if (.not. so_print) return 3991c 3992c print basis set information 3993c 3994 write(LuOut,1)bs_name(soid)(1:inp_strlen(bs_name(soid))), 3995 $ bs_trans(soid)(1:inp_strlen(bs_trans(soid))) 3996 1 format(' SO Potential "',a,'" -> "',a,'"'/ 3997 $ ' -----') 3998 if (bas_spherical(soid)) write(LuOut,2) 3999 2 format(' SO Basis is spherical, 5d, 7f, 9g ... ') 4000 mytags = infbs_head(HEAD_NTAGS,soid) 4001 if (mytags.le.0) then 4002 write(LuOut,*)'No explicit SO ECP functions are defined !' 4003 write(LuOut,*) 4004c 4005c there could be star tags defined, so check that before returning 4006c 4007 goto 00010 4008 endif 4009c 4010 myucont = infbs_head(HEAD_NCONT,soid) 4011 myprim = infbs_head(HEAD_NPRIM,soid) 4012 mycoef = infbs_head(HEAD_NCOEF,soid) 4013c 4014 do 00100 i=1,mytags 4015 4016 if (geom_tag_to_element(bs_tags(i,soid), symbol, element, 4017 $ atn)) then 4018 len_tag = inp_strlen(bs_tags(i,soid)) 4019 len_ele = inp_strlen(element) 4020 write (buffer, 4021 & '(a,'' ('',a,'')'' )') 4022 & bs_tags(i,soid)(1:len_tag), element(1:len_ele) 4023 else 4024 buffer = bs_tags(i,soid) 4025 endif 4026 len_tag = inp_strlen(buffer) 4027 call util_print_centered(LuOut, buffer, len_tag/2 + 1, .true.) 4028 4029 myucont = infbs_tags(TAG_NCONT,i,soid) 4030c 4031 ifcont = infbs_tags(TAG_FCONT,i,soid) 4032c 4033 write(LuOut,6) 4034 6 format( 4035 $ ' R-exponent Exponent Coefficients '/ 4036 $ ' ------------ ',60('-')) 4037 do 00200 j=1,myucont 4038 myprim = infbs_cont(CONT_NPRIM,ifcont,soid) 4039 mygen = infbs_cont(CONT_NGEN,ifcont,soid) 4040 4041 mytype = infbs_cont(CONT_TYPE, ifcont, soid) 4042 shell_type = ctype(mytype) 4043 iexptr = infbs_cont(CONT_IEXP, ifcont,soid) - 1 4044 icfptr = infbs_cont(CONT_ICFP, ifcont,soid) - 1 4045 irexptr = infbs_cont(CONT_IREXP,ifcont,soid) - 1 4046 do 00300 k=1,myprim 4047 write(LuOut,7) j, shell_type, 4048 & sf_exndcf((irexptr+k),soid), 4049 & sf_exndcf((iexptr+k),soid), 4050 & (sf_exndcf((icfptr+k+(l-1)*myprim),soid),l=1,mygen) 405100300 continue 4052 write(LuOut,*) 4053 ifcont = ifcont + 1 405400200 continue 405500100 continue 4056c 4057c Check if we have star tag definitions in the basis set 4058c 405900010 if (star_nr_tags .gt. 0) then 4060 write(LuOut,*) 4061 write(LuOut,8) 4062 8 format(' In addition, one or more string tags have been', 4063 & ' defined containing a * .'/' These tags, and ', 4064 & 'their exceptions list are printed below.'// 4065 & ' Tag ',12(' '),' Description ',18(' '), 4066 & ' Excluding '/' ',16('-'),' ',30('-'),' ',30('-')) 4067 do i=1,star_nr_tags 4068 k = 1 4069 if (i .gt. 1) k = star_nr_excpt(i-1) + 1 4070 write(LuOut,9) star_tag(i), star_bas_typ(i)(1:30), 4071 & (star_excpt(j)(1:inp_strlen(star_excpt(j))), 4072 & j=k,star_nr_excpt(i)) 4073 9 format(1x,a16,1x,a30,1x,10(a)) 4074 enddo 4075 write(LuOut,*) 4076 write(LuOut,*) 4077 endif 4078c 4079c If geom is set print out the info about total basis info 4080c associated with the geometry also 4081c 4082c ... not done yet 4083c 4084 return 4085 7 format(1x,i2,1x,a3,1x,f9.2,2x,f14.6,20f15.6) 4086 end 4087*..................................................................... 4088C> 4089C> \brief Set the name of the ECP for a basis set 4090C> 4091C> \return Return .true. when successfull, and .false. otherwise 4092c 4093 logical function bas_set_ecp_name(basisin,ecp_name) 4094 implicit none 4095#include "nwc_const.fh" 4096#include "basP.fh" 4097 integer basisin !< [Input] the basis set handle 4098 character*(*) ecp_name !< [Input] the ECP name 4099c 4100 name_assoc(1,(basisin+Basis_Handle_Offset)) = 4101 & ecp_name 4102 bas_set_ecp_name = .true. 4103 end 4104*..................................................................... 4105C> 4106C> \brief Set the name of the spin-orbit potential for a basis set 4107C> 4108C> \return Return .true. when successfull, and .false. otherwise 4109c 4110 logical function bas_set_so_name(basisin,so_name) 4111 implicit none 4112#include "nwc_const.fh" 4113#include "basP.fh" 4114 integer basisin !< [Input] the basis set handle 4115 character*(*) so_name !< [Input] the spin-orbit potential name 4116c 4117 name_assoc(2,(basisin+Basis_Handle_Offset)) = 4118 & so_name 4119 bas_set_so_name = .true. 4120 end 4121*..................................................................... 4122C> 4123C> \brief Retrieve the name of the ECP of a basis set 4124C> 4125C> \return Return .true. when successfull, and .false. otherwise 4126c 4127 logical function bas_get_ecp_name(basisin,ecp_name) 4128 implicit none 4129#include "nwc_const.fh" 4130#include "basP.fh" 4131 integer basisin !< [Input] the basis set handle 4132 character*(*) ecp_name !< [Output] the ECP name 4133c 4134 ecp_name = 4135 & name_assoc(1,(basisin+Basis_Handle_Offset)) 4136 bas_get_ecp_name = .true. 4137 end 4138*..................................................................... 4139C> 4140C> \brief Retrieve the name of the spin-orbit potential of a basis set 4141C> 4142C> \return Return .true. when successfull, and .false. otherwise 4143c 4144 logical function bas_get_so_name(basisin,so_name) 4145 implicit none 4146#include "nwc_const.fh" 4147#include "basP.fh" 4148 integer basisin !< [Input] the basis set handle 4149 character*(*) so_name !< [Output] the spin-orbit potential name 4150c 4151 so_name = 4152 & name_assoc(2,(basisin+Basis_Handle_Offset)) 4153 bas_get_so_name = .true. 4154 end 4155*..................................................................... 4156C> 4157C> \brief Associate an ECP with a given basis set 4158C> 4159C> \return Return .true. if successfull, and .false. otherwise 4160C> 4161 logical function bas_set_ecp_handle(basisin,ecp_handle) 4162 implicit none 4163#include "nwc_const.fh" 4164#include "basP.fh" 4165 integer basisin !< [Input] the molecular basis set handle 4166 integer ecp_handle !< [Input] the EPC handle 4167c 4168 handle_assoc(1,(basisin+Basis_Handle_Offset)) = 4169 & ecp_handle 4170 bas_set_ecp_handle = .true. 4171 end 4172*..................................................................... 4173C> 4174C> \brief Retrieve an ECP of a given basis set 4175C> 4176C> \return Return .true. if successfull, and .false. otherwise 4177C> 4178 logical function bas_get_ecp_handle(basisin,ecp_handle) 4179 implicit none 4180#include "nwc_const.fh" 4181#include "basP.fh" 4182 integer basisin !< [Input] the molecular basis set handle 4183 integer ecp_handle !< [Output] the ECP handle 4184c 4185 ecp_handle = 4186 & handle_assoc(1,(basisin+Basis_Handle_Offset)) 4187 if (ecp_handle.ne.0) then 4188 bas_get_ecp_handle = .true. 4189 else 4190 bas_get_ecp_handle = .false. 4191 endif 4192 end 4193*..................................................................... 4194C> 4195C> \brief Set the molecular basis set for an ECP 4196C> 4197C> \return Return .true. if successfull, and .false. otherwise 4198C> 4199 logical function ecp_set_parent_handle(ecp_handle,basis_handle) 4200 implicit none 4201#include "errquit.fh" 4202#include "nwc_const.fh" 4203#include "basP.fh" 4204 integer ecp_handle !< [Input] the ECP handle 4205 integer basis_handle !< [Input] the molecular basis set handle 4206c 4207 logical bas_check_handle, ecp_check_handle 4208 external bas_check_handle, ecp_check_handle 4209c 4210 if (.not.bas_check_handle(basis_handle,'ecp_set_parent_handle')) 4211 & call errquit 4212 & ('ecp_set_parent_handle: basis_handle invalid',911, 4213 & BASIS_ERR) 4214 if (.not.ecp_check_handle(ecp_handle,'ecp_set_parent_handle')) 4215 & call errquit 4216 & ('ecp_set_parent_handle: not an ECP ...'// 4217 $ ' check for name conflict between basis and ECP',911, 4218 & BASIS_ERR) 4219 parent_assoc(1,(ecp_handle+Basis_Handle_Offset)) = 4220 & basis_handle 4221 bas_nassoc(basis_handle+Basis_Handle_Offset) = 4222 & bas_nassoc(basis_handle+Basis_Handle_Offset) + 1 4223 ecp_set_parent_handle = .true. 4224 end 4225*..................................................................... 4226C> 4227C> \brief Retrieve the molecular basis set of an ECP 4228C> 4229C> \return Return .true. if successfull, and .false. otherwise 4230C> 4231 logical function ecp_get_parent_handle(ecp_handle,basis_handle) 4232 implicit none 4233#include "errquit.fh" 4234#include "nwc_const.fh" 4235#include "basP.fh" 4236 integer ecp_handle !< [Input] the ECP handle 4237 integer basis_handle !< [Output] the molecular basis set handle 4238c 4239 logical bas_check_handle, ecp_check_handle 4240 external bas_check_handle, ecp_check_handle 4241c 4242 if (.not.ecp_check_handle(ecp_handle,'ecp_get_parent_handle')) 4243 & call errquit 4244 & ('ecp_get_parent_handle: ecp_handle invalid',911, 4245 & BASIS_ERR) 4246 basis_handle = parent_assoc(1,ecp_handle+Basis_Handle_Offset) 4247 if (basis_handle.ne.0) then 4248 ecp_get_parent_handle = .true. 4249 else 4250 ecp_get_parent_handle = .false. 4251 endif 4252 if (.not.bas_check_handle(basis_handle,'ecp_get_parent_handle')) 4253 & call errquit 4254 & ('ecp_get_parent_handle: stored basis_handle invalid',911, 4255 & BASIS_ERR) 4256 end 4257*..................................................................... 4258C> 4259C> \brief Set the spin-orbit potential for a molecular basis set 4260C> 4261C> \return Return .true. if successfull, and .false. otherwise 4262C> 4263 logical function bas_set_so_handle(basisin,so_handle) 4264 implicit none 4265#include "nwc_const.fh" 4266#include "basP.fh" 4267 integer basisin !< [Input] the molecular basis set handle 4268 integer so_handle !< [Input] the spin-orbit potential handle 4269c 4270 handle_assoc(2,(basisin+Basis_Handle_Offset)) = 4271 & so_handle 4272 bas_set_so_handle = .true. 4273 end 4274*..................................................................... 4275C> 4276C> \brief Retrieve the spin-orbit potential of a molecular basis set 4277C> 4278C> \return Return .true. if successfull, and .false. otherwise 4279C> 4280 logical function bas_get_so_handle(basisin,so_handle) 4281 implicit none 4282#include "nwc_const.fh" 4283#include "basP.fh" 4284 integer basisin !< [Input] the molecular basis set handle 4285 integer so_handle !< [Output] the spin-orbit potential handle 4286c 4287 so_handle = 4288 & handle_assoc(2,(basisin+Basis_Handle_Offset)) 4289 if (so_handle.ne.0) then 4290 bas_get_so_handle = .true. 4291 else 4292 bas_get_so_handle = .false. 4293 endif 4294 end 4295*..................................................................... 4296C> 4297C> \brief Set the molecular basis set for a spin-orbit potential 4298C> 4299C> \return Return .true. if successfull, and .false. otherwise 4300C> 4301 logical function so_set_parent_handle(so_handle,basis_handle) 4302 implicit none 4303#include "errquit.fh" 4304#include "nwc_const.fh" 4305#include "basP.fh" 4306 integer so_handle !< [Input] the spin-orbit potential handle 4307 integer basis_handle !< [Input] the molecular basis set handle 4308c 4309 logical bas_check_handle, so_check_handle 4310 external bas_check_handle, so_check_handle 4311c 4312 if (.not.bas_check_handle(basis_handle,'so_set_parent_handle')) 4313 & call errquit 4314 & ('so_set_parent_handle: basis_handle invalid',911, 4315 & BASIS_ERR) 4316 if (.not.so_check_handle(so_handle,'so_set_parent_handle')) 4317 & call errquit 4318 & ('so_set_parent_handle: so_handle invalid',911, 4319 & BASIS_ERR) 4320 parent_assoc(2,(so_handle+Basis_Handle_Offset)) = 4321 & basis_handle 4322 bas_nassoc(basis_handle+Basis_Handle_Offset) = 4323 & bas_nassoc(basis_handle+Basis_Handle_Offset) + 1 4324 so_set_parent_handle = .true. 4325 end 4326*..................................................................... 4327C> 4328C> \brief Retrieve the molecular basis set of a spin-orbit potential 4329C> 4330C> \return Return .true. if successfull, and .false. otherwise 4331C> 4332 logical function so_get_parent_handle(so_handle,basis_handle) 4333 implicit none 4334#include "errquit.fh" 4335#include "nwc_const.fh" 4336#include "basP.fh" 4337 integer so_handle !< [Input] the spin-orbit potential handle 4338 integer basis_handle !< [Output] the molecular basis set handle 4339c 4340 logical bas_check_handle, so_check_handle 4341 external bas_check_handle, so_check_handle 4342c 4343 if (.not.so_check_handle(so_handle,'so_get_parent_handle')) 4344 & call errquit 4345 & ('so_get_parent_handle: so_handle invalid',911, BASIS_ERR) 4346 basis_handle = 4347 & parent_assoc(2,(so_handle+Basis_Handle_Offset)) 4348 if (basis_handle.ne.0) then 4349 so_get_parent_handle = .true. 4350 else 4351 so_get_parent_handle = .false. 4352 endif 4353 if (.not.bas_check_handle(basis_handle,'so_get_parent_handle')) 4354 & call errquit 4355 & ('so_get_parent_handle: stored basis_handle invalid',911, 4356 & BASIS_ERR) 4357 end 4358C 4359C> \brief Summarize the internal data structures 4360C 4361 subroutine bas_print_allocated_info(msg) 4362 implicit none 4363#include "stdio.fh" 4364#include "mafdecls.fh" 4365#include "nwc_const.fh" 4366#include "basP.fh" 4367#include "geobasmapP.fh" 4368#include "global.fh" 4369#include "bas_ibs_dec.fh" 4370#include "bas_exndcf_dec.fh" 4371 character*(*) msg !< [Input] A message to identify the print 4372 integer basis 4373 integer me, nproc, inode 4374 integer inode_use 4375c 4376#include "bas_ibs_sfn.fh" 4377#include "bas_exndcf_sfn.fh" 4378c 4379 me = ga_nodeid() 4380 nproc = ga_nnodes() 4381* 4382 do inode = 0,(nproc-1) 4383 inode_use = inode 4384 call ga_sync() 4385 call ga_brdcst(1234,inode_use, 4386 & ma_sizeof(mt_int, 1, mt_byte),0) 4387 call ga_sync() 4388 if (inode_use.eq.me) then 4389 write(luout,*)' msg: ',msg,me 4390 write(luout,*)' basis data for node ',me 4391 write(luout,*)' number of possible basis sets',nbasis_bsmx 4392 do basis = 1,nbasis_bsmx 4393 write(luout,*)' active :',bsactive(basis) 4394 write(luout,*)' exndcf Handle/index/size :', 4395 & exndcf(H_exndcf,basis),exndcf(K_exndcf,basis), 4396 & exndcf(SZ_exndcf,basis) 4397 write(luout,*)' cn2ucn (H/K/SZ) :', 4398 & ibs_cn2ucn(H_ibs,basis), 4399 & ibs_cn2ucn(K_ibs,basis), 4400 & ibs_cn2ucn(SZ_ibs,basis) 4401 write(luout,*)' cn2ce (H/K/SZ) :', 4402 & ibs_cn2ce(H_ibs,basis), 4403 & ibs_cn2ce(K_ibs,basis), 4404 & ibs_cn2ce(SZ_ibs,basis) 4405 write(luout,*)' ce2uce (H/K/SZ) :', 4406 & ibs_ce2uce(H_ibs,basis), 4407 & ibs_ce2uce(K_ibs,basis), 4408 & ibs_ce2uce(SZ_ibs,basis) 4409 write(luout,*)' cn2bfr (H/K/SZ) :', 4410 & ibs_cn2bfr(H_ibs,basis), 4411 & ibs_cn2bfr(K_ibs,basis), 4412 & ibs_cn2bfr(SZ_ibs,basis) 4413 write(luout,*)' ce2cnr (H/K/SZ) :', 4414 & ibs_ce2cnr(H_ibs,basis), 4415 & ibs_ce2cnr(K_ibs,basis), 4416 & ibs_ce2cnr(SZ_ibs,basis) 4417 enddo 4418 call util_flush(luout) 4419 endif 4420 enddo 4421 end 4422c 4423C> \brief Matches a tag from the geometry to a tag in the basis set 4424C> 4425C> Associated with a shell of basis functions is a tag identifying the 4426C> kind of center the shell should reside on. These centers should of 4427C> course appear in the geometry. This routine takes a tag from the 4428C> geometry and matches it to a tag from the basis set. The index of 4429C> the first tag in the basis set that matches the geometry tag is 4430C> retrieved. 4431C> 4432C> \return Return .true. when a match was found, and .false. otherwise 4433c 4434 logical function bas_match_tags(tag_from_geom,basisin,btag) 4435 implicit none 4436#include "errquit.fh" 4437* 4438* This routine matches geometry tags to basis set tags in such a way 4439* that "H34" matches "H" or "hydrogen" or "h" if "H34" does not 4440* exist in the basis set specification. 4441* "H" will however not match "h" i.e., case sensitivity is 4442* preserved. 4443* 4444c::includes 4445#include "stdio.fh" 4446#include "inp.fh" 4447#include "nwc_const.fh" 4448#include "basP.fh" 4449#include "basdeclsP.fh" 4450c::functions 4451 logical geom_tag_to_element 4452 external geom_tag_to_element 4453c::passed 4454 character*(*) tag_from_geom !< [Input] the geometry tag 4455 integer basisin !< [Input] the basis set handle 4456 integer btag !< [Output] lexical basis set tag index 4457* bas_match_tags ! [output] true if they match 4458c::local 4459 integer basis ! lexical basis index 4460 integer nbtgs ! number of basis set tags 4461 character*16 gstring ! geometry or basis tags can only be *16 4462 character*16 bsmatch ! the matched basis set tag 4463 integer lgstring ! length of gstring 4464 integer lgstring_old ! length of gstring at first test 4465 integer lnotalpha 4466 integer i ! loop counter 4467 integer ind ! dummy counter 4468*rak: logical digits_found ! did geometry tag have any digits 4469 logical status ! dummy storage 4470 character*2 g_sym ! geometry tag -> symbol name 4471 character*16 g_elem ! geometry tag -> element name 4472 integer g_atn ! geometry tag -> atomic number 4473 character*2 b_sym ! basis set tag -> symbol name 4474 character*16 b_elem ! basis set tag -> element name 4475 integer b_atn ! basis set tag -> atomic number 4476 logical debug ! true for extra output 4477c 4478 integer na2z 4479 parameter (na2z = 26) 4480 character*1 a2z(na2z) 4481 data a2z /'a','b','c','d','e','f','g','h','i','j', 4482 & 'k','l','m','n','o','p','q','r','s','t', 4483 & 'u','v','w','x','y','z'/ 4484*rak: integer ndigits ! number of digits 4485*rak: parameter (ndigits=10) ! set number of digits 4486*rak: character*1 digits(ndigits) ! array of character digits 4487*rak:c 4488*rak: data digits /'0','1','2','3','4','5','6','7','8','9'/ 4489c 4490 debug = .false. 4491c 4492 bas_match_tags = .false. 4493 bsmatch = ' ' 4494c 4495 basis = basisin+Basis_Handle_Offset 4496c 4497 gstring = tag_from_geom ! copy tag to work on it 4498c 4499 if (gstring(1:4).eq.'bqgh') then 4500 gstring = gstring(5:) 4501 endif 4502c 4503c... first match full geometry tag to full basis tag list 4504c did the user specifically assign a tag ? 4505 4506 nbtgs = infbs_head(HEAD_NTAGS,basis) 4507 do i = 1, nbtgs 4508 if (gstring.eq.bs_tags(i,basis)) then 4509 bas_match_tags = .true. 4510 btag = i 4511 bsmatch = bs_tags(i,basis) 4512 goto 00009 4513 endif 4514 enddo 4515c 4516c... Now check to see if the tag has non alpha chracters. 4517c... the substring prior to the first non-alpha character is the users idea 4518c... of the "name" of the tag. 4519c 4520 lgstring = inp_strlen(gstring) 4521 lnotalpha = lgstring+1 4522 do i = 1,lgstring 4523 if (.not.(inp_match(na2z,.false., 4524 & gstring(i:i),a2z,ind))) then ! compare character to alpha (case-less test) 4525 lnotalpha = i 4526 goto 00001 4527 endif 4528 enddo 452900001 continue 4530 do i = lnotalpha,lgstring 4531 gstring(i:i) = ' ' 4532 enddo 4533 if (debug) write(luout,*)' the string:', tag_from_geom, 4534 & 'has the user substring of ',gstring 4535*rak: 4536*rak: digits_found = .false. 4537*rak:00001 continue 4538*rak: lgstring = inp_strlen(gstring) ! get the length of the geometry tag 4539*rak: if (lgstring.eq.0) return ! if empty string or string of digits return false 4540*rak: if (inp_match(ndigits,.false., 4541*rak: & gstring(lgstring:lgstring),digits,i)) then ! compare last character to a digit 4542*rak: gstring(lgstring:lgstring) = ' ' ! if a digit remove it 4543*rak: digits_found = .true. 4544*rak: goto 00001 4545*rak: endif 4546*rak:c 4547*rak:c if no digits then enforce case matching between say "H" and "h" by a false return 4548*rak:c 4549*rak: if (.not.digits_found) return 4550*rak:c 4551*rak:c... first match numberless geometry tag to full basis tag list 4552c 4553c... match user substring to basis set tags (only if it gstring has a different length) 4554 lgstring_old = lgstring 4555 lgstring = inp_strlen(gstring) 4556 if (lgstring_old.gt.lgstring) then 4557 do i = 1, nbtgs 4558 if (gstring.eq.bs_tags(i,basis)) then 4559 bas_match_tags = .true. 4560 btag = i 4561 bsmatch = bs_tags(i,basis) 4562 goto 00009 4563 endif 4564 enddo 4565 endif 4566c 4567c... now get symbol and element names for each tag and match those 4568c geometry tag is based on users substring 4569c basis set tag is based on user input to the basis object 4570c... bq's with basis functions must match from the above rules! 4571c 4572 status = geom_tag_to_element(gstring,g_sym,g_elem,g_atn) 4573 if (.not.status)then 4574 if (.not.(g_sym.eq.'bq')) then 4575 write(luout,*)'geometry tag <',gstring, 4576 & '> could not be matched to an element symbol' 4577 call errquit('bas_match_tags: fatal error ',911, BASIS_ERR) 4578 endif 4579 endif 4580 if (g_sym.eq.'bq') goto 00009 ! bq labes with basis functions must match from above 4581 do i = 1, nbtgs 4582 status = 4583 & geom_tag_to_element(bs_tags(i,basis),b_sym,b_elem,b_atn) 4584 if (.not.status) then 4585 if (.not.(b_sym.eq.'bq')) then 4586 write(luout,*)'basis tag',bs_tags(i,basis), 4587 & ' could not be matched to an element symbol' 4588 call errquit('bas_match_tags: fatal error ',911, BASIS_ERR) 4589 endif 4590 endif 4591 if (g_elem.eq.b_elem) then 4592 bas_match_tags = .true. 4593 btag = i 4594 bsmatch = bs_tags(i,basis) 4595 goto 00009 4596 else if (g_sym.eq.b_sym) then 4597 bas_match_tags = .true. 4598 btag = i 4599 bsmatch = bs_tags(i,basis) 4600 goto 00009 4601 endif 4602 enddo 4603 if (debug) 4604 & write(luout,*)'bas_match_tags:debug: no match for tag <', 4605 & tag_from_geom,'>' 4606 return ! no match 4607c 460800009 continue 4609 if (debug) then 4610 write(luout,10000)tag_from_geom,bsmatch 4611 endif 4612 return 461310000 format('bas_match_tags:debug: geometry tag ',a16, 4614 & ' matched basis set tag ',a16) 4615 end 4616c 4617C> \brief Work out whether multipoles can be calculated for a basis set 4618C> 4619C> The multipole code can expand the molecular potential in multipoles. 4620C> However, it cannot handle all basis sets. In particular it cannot 4621C> handle generally contracted or SP shells. 4622C> 4623C> \return Return .true. if this basis set is compatible with the 4624C> multipole code, and .false. otherwise 4625c 4626 logical function bas_cando_mpoles(basis) 4627 implicit none 4628#include "errquit.fh" 4629 integer basis !< [Input] the basis set handle 4630c 4631c Return true if it's possible to compute multipoles for this basis 4632c 4633c The multipole code cannot handle general contractions or SP shells. 4634c 4635 integer nshell, ishell, type, nprim, ngen, sphcart 4636 logical bas_numcont, bas_continfo 4637 external bas_numcont, bas_continfo 4638c 4639 bas_cando_mpoles = .false. 4640c 4641 if (.not. bas_numcont(basis, nshell)) call errquit 4642 $ ('multipole: basis bad?',0, BASIS_ERR) 4643 do ishell = 1, nshell 4644 if (.not. bas_continfo(basis, ishell, type, 4645 $ nprim, ngen, sphcart)) call errquit 4646 $ ('multipole: basis bad?',0, BASIS_ERR) 4647 if (type.lt.0 .or. ngen.gt.1) return 4648 enddo 4649c 4650 bas_cando_mpoles = .true. 4651c 4652 end 4653c 4654C> \brief Lookup whether a basis set consists of spherical harmonic 4655C> functions 4656C> 4657C> At present cartesian and spherical harmonic basis sets are supported. 4658C> This routine looks up whether the specified basis set is a spherical 4659C> harmonic basis set. 4660C> 4661C> \return Return .true. if the basis set consists of spherical 4662C> harmonic functions, and .false. otherwise 4663C 4664 logical function bas_is_spherical(basisin) 4665 implicit none 4666#include "errquit.fh" 4667#include "nwc_const.fh" 4668#include "basP.fh" 4669#include "stdio.fh" 4670*::passed 4671 integer basisin !< [Input] the basis set handle 4672*::local 4673 integer basis ! lexical index of basis set handle 4674*::functions 4675 logical bas_check_handle 4676 external bas_check_handle 4677*::code 4678 bas_is_spherical = bas_check_handle(basisin,'bas_is_spherical') 4679 if (.not.bas_is_spherical) then 4680 write(luout,*)'invalid/inactive basis handle:',basisin 4681 call errquit('bas_is_spherical: fatal error',911, BASIS_ERR) 4682 endif 4683 basis = basisin + Basis_Handle_Offset 4684 bas_is_spherical = bas_spherical(basis) 4685 return 4686* 4687 end 4688C> @} 4689 subroutine bas_dump_info(msg) 4690 implicit none 4691#include "stdio.fh" 4692#include "inp.fh" 4693#include "nwc_const.fh" 4694#include "basP.fh" 4695#include "basdeclsP.fh" 4696#include "geobasmapP.fh" 4697 character*(*) msg 4698c 4699 integer llen 4700 integer basis 4701 integer ntags 4702 llen = inp_strlen(msg) 4703 write(luout,00001) 4704 write(luout,*)'bas_dump_info: <',msg(1:llen),'>' 4705 write(luout,*)'bas_dump_info: bsversion:', bsversion 4706 do basis = 1, nbasis_bsmx 4707 if (bsactive(basis)) then 4708 write(luout,*)'---------------------------------------------' 4709 write(luout,*)'---------------------------------------------' 4710 write(luout,*)'handle(computed) :', 4711 & (basis-basis_handle_offset) 4712 write(luout,*)'geometry loader :', 4713 & ibs_geom(basis) 4714 write(luout,*)'bs_name :', 4715 & bs_name(basis)(1:len_bs_name(basis)) 4716 write(luout,*)'bs_trans :', 4717 & bs_trans(basis)(1:len_bs_trans(basis)) 4718 llen = inp_strlen(name_assoc(1,basis)) 4719 write(luout,*)'ecp_name_assoc :', 4720 & name_assoc(1,basis)(1:llen) 4721 llen = inp_strlen(name_assoc(2,basis)) 4722 write(luout,*)'so_name_assoc :', 4723 & name_assoc(2,basis)(1:llen) 4724 ntags = infbs_head(Head_Ntags,basis) 4725 write(luout,*)'head ntags :',ntags 4726 write(luout,*)'head ncont :', 4727 & infbs_head(Head_Ncont,basis) 4728 write(luout,*)'head nprim :', 4729 & infbs_head(Head_Nprim,basis) 4730 write(luout,*)'head ncoef :', 4731 & infbs_head(Head_Ncoef,basis) 4732 write(luout,*)'head excfptr :', 4733 & infbs_head(Head_Excfptr,basis) 4734 write(luout,*)'head sph :', 4735 & infbs_head(Head_Sph,basis) 4736 write(luout,*)'head ecp :', 4737 & infbs_head(Head_Ecp,basis) 4738 write(luout,*)'bas_spherical :', 4739 & bas_spherical(basis) 4740 write(luout,*)'bas_any_gc :', 4741 & bas_any_gc(basis) 4742 write(luout,*)'bas_any_sp_shell :', 4743 & bas_any_sp_shell(basis) 4744 write(luout,*)'bas_norm_id :', 4745 & bas_norm_id(basis) 4746 write(luout,*)'angular_bs :', 4747 & angular_bs(basis) 4748 write(luout,*)'nbfmax_bs :', 4749 & nbfmax_bs(basis) 4750 write(luout,*)'ecp_handle_assoc :', 4751 & handle_assoc(1,basis) 4752 write(luout,*)'ecp_parent_assoc :', 4753 & parent_assoc(1,basis) 4754 write(luout,*)'so_handle_assoc :', 4755 & handle_assoc(2,basis) 4756 write(luout,*)'so_parent_assoc :', 4757 & parent_assoc(2,basis) 4758 write(luout,*)'num of assoc :', 4759 & bas_nassoc(basis) 4760 write(luout,*)'^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^' 4761 write(luout,*)'---------------------------------------------' 4762 endif 4763 enddo 4764 write(luout,00001) 476500001 format(1x,80('=')) 4766 end 4767*..................................................................... 4768C> \ingroup bas 4769C> @{ 4770C> 4771C> \brief Find the length of the longest contraction in a basis set 4772C> 4773C> Retrieves the length of the longest contraction (i.e. the maximum 4774C> number of Gaussian primitive terms) in the specified basis set, 4775C> for any contraction/shell, atom, etc. 4776C> 4777C> \return Return .true. if successfull, and .false. otherwise 4778C 4779 Logical function bas_ncontr_cn_max( basis_hand, ncontr_max ) 4780c 4781c returns the largest number of contractions, i.e., columns of 4782c contraction coefficients for a given contraction set, 4783c in the basis input, for any contraction set in the basis 4784c 4785 implicit none 4786 logical bas_numcont, bas_continfo 4787 external bas_numcont, bas_continfo 4788c::args 4789 integer basis_hand !< [Input] the basis set handle 4790 integer ncontr_max !< [Output] the maximum contraction length 4791c::locals 4792 logical LResult 4793 integer itype, nprimo, isphcart 4794 integer ncontr, icontset, numcontset 4795c::exec 4796 LResult = .true. 4797 LResult = LResult .and. bas_numcont(basis_hand, numcontset) 4798c 4799c loop over contraction sets; find largest number of contractions 4800c (i.e., columns of contraction coefficients) 4801c 4802 ncontr_max = 0 4803 do icontset = 1,numcontset 4804 LResult = LResult .and. 4805 & bas_continfo( basis_hand, icontset, itype, nprimo, 4806 & ncontr, isphcart) 4807 4808 ncontr_max = max( ncontr_max, ncontr ) 4809 enddo 4810 4811 bas_ncontr_cn_max = LResult 4812 4813 return 4814 end 4815*---------------------------------------------------------------------- 4816C> \brief Retrieve the number of centers with either ECPs or spin-orbit 4817C> potentials 4818C> 4819C> \return Return .true. if successfull, and .false. otherwise 4820c 4821 logical function ecpso_ncent(geom,soidin,ecpidin,num_cent) 4822 implicit none 4823#include "errquit.fh" 4824* 4825* return the combined number of ecp/so centers 4826* 4827*::includes 4828#include "geom.fh" 4829*::functions 4830 logical bas_match_tags 4831 external bas_match_tags 4832*::passed 4833 integer geom !< [Input] the geometry handle 4834 integer soidin !< [Input] the spin-orbit potential (so) handle 4835 integer ecpidin !< [Input] the ECP handle 4836 integer num_cent !< [Output] the number of ECP/so centers 4837*::local 4838 integer ntotal ! total number of centers in goemetry 4839 integer ic ! loop counter for centers 4840 integer tag_indx ! basis set unique center index for tag 4841 character*16 gtag ! geometry tag 4842 double precision cxyz(3) ! coordinates 4843 double precision chg ! charge 4844 logical inthe_ecp ! is tag in the ecp? 4845 logical inthe_so ! is the tag in the so pot.? 4846* 4847 ecpso_ncent = geom_ncent(geom,ntotal) 4848 if (.not.ecpso_ncent) call errquit 4849 & ('ecpso_ncent: geom_ncent failed?',911, GEOM_ERR) 4850 num_cent = 0 4851 do ic = 1,ntotal 4852 if (.not.geom_cent_get(geom,ic,gtag,cxyz,chg)) call errquit 4853 & ('ecpso_ncent: geom_cent_get failed',911, GEOM_ERR) 4854 inthe_ecp = bas_match_tags(gtag,ecpidin,tag_indx) 4855 inthe_so = bas_match_tags(gtag,soidin,tag_indx) 4856 if (inthe_ecp.or.inthe_so) num_cent = num_cent + 1 4857 enddo 4858 end 4859*---------------------------------------------------------------------- 4860C> \brief Retrieve the number of centers with spin-orbit potentials 4861C> \return Return .true. if successfull, and .false. otherwise 4862 logical function so_ncent(geom,soidin,num_cent) 4863 implicit none 4864#include "errquit.fh" 4865* 4866* return the number of so centers 4867* 4868*::includes 4869#include "geom.fh" 4870*::functions 4871 logical bas_match_tags 4872 external bas_match_tags 4873*::passed 4874 integer geom !< [Input] the geometry handle 4875 integer soidin !< [Input] the spin-orbit potential (so) handle 4876 integer num_cent !< [Output] the number of so centers 4877*::local 4878 integer ntotal ! total number of centers in goemetry 4879 integer ic ! loop counter for centers 4880 integer tag_indx ! basis set unique center index for tag 4881 character*16 gtag ! geometry tag 4882 double precision cxyz(3) ! coordinates 4883 double precision chg ! charge 4884 logical inthe_so ! is the tag in the so pot.? 4885* 4886 so_ncent = geom_ncent(geom,ntotal) 4887 if (.not.so_ncent) call errquit 4888 & ('so_ncent: geom_ncent failed?',911, GEOM_ERR) 4889 num_cent = 0 4890 do ic = 1,ntotal 4891 if (.not.geom_cent_get(geom,ic,gtag,cxyz,chg)) call errquit 4892 & ('so_ncent: geom_cent_get failed',911, GEOM_ERR) 4893 inthe_so = bas_match_tags(gtag,soidin,tag_indx) 4894 if (inthe_so) num_cent = num_cent + 1 4895 enddo 4896 end 4897*---------------------------------------------------------------------- 4898C> \brief Retrieve the number of centers with ECPs 4899C> \return Return .true. if successfull, and .false. otherwise 4900 logical function ecp_ncent(geom,ecpidin,num_cent) 4901 implicit none 4902#include "errquit.fh" 4903* 4904* return the number of ecp centers 4905* 4906*::includes 4907#include "geom.fh" 4908*::functions 4909 logical bas_match_tags 4910 external bas_match_tags 4911*::passed 4912 integer geom !< [Input] the geometry handle 4913 integer ecpidin !< [Input] the ECP handle 4914 integer num_cent !< [Output] the number of ECP centers 4915*::local 4916 integer ntotal ! total number of centers in goemetry 4917 integer ic ! loop counter for centers 4918 integer tag_indx ! basis set unique center index for tag 4919 character*16 gtag ! geometry tag 4920 double precision cxyz(3) ! coordinates 4921 double precision chg ! charge 4922 logical inthe_ecp ! is the tag in the ecp pot.? 4923* 4924 ecp_ncent = geom_ncent(geom,ntotal) 4925 if (.not.ecp_ncent) call errquit 4926 & ('ecp_ncent: geom_ncent failed?',911, GEOM_ERR) 4927 num_cent = 0 4928 do ic = 1,ntotal 4929 if (.not.geom_cent_get(geom,ic,gtag,cxyz,chg)) call errquit 4930 & ('ecp_ncent: geom_cent_get failed',911, GEOM_ERR) 4931 inthe_ecp = bas_match_tags(gtag,ecpidin,tag_indx) 4932 if (inthe_ecp) num_cent = num_cent + 1 4933 enddo 4934 end 4935*---------------------------------------------------------------------- 4936C> \brief Retrieve the list of centers with either ECP or spin-orbit 4937C> potentials 4938C> 4939C> \return Return .true. if successfull, and .false. otherwise 4940c 4941 logical function ecpso_list_ncent(geom,soidin,ecpidin, 4942 & num_cent,ecpso_list) 4943 implicit none 4944#include "errquit.fh" 4945* 4946* return the combined number of ecp/so centers and the list of centers 4947* 4948*::includes 4949#include "geom.fh" 4950*::functions 4951 logical bas_match_tags 4952 external bas_match_tags 4953*::passed 4954 integer geom !< [Input] the geometry handle 4955 integer soidin !< [Input] the spin-orbit potential handle 4956 integer ecpidin !< [Input] the ECP handle 4957 integer num_cent !< [Output] the number of ECP/so centers 4958 integer ecpso_list(*) !< [Output] the list of ECP/so centers 4959*::local 4960 integer ntotal ! total number of centers in goemetry 4961 integer ic ! loop counter for centers 4962 integer tag_indx ! basis set unique center index for tag 4963 character*16 gtag ! geometry tag 4964 double precision cxyz(3) ! coordinates 4965 double precision chg ! charge 4966 logical inthe_ecp ! is tag in the ecp? 4967 logical inthe_so ! is the tag in the so pot.? 4968* 4969 ecpso_list_ncent = geom_ncent(geom,ntotal) 4970 if (.not.ecpso_list_ncent) call errquit 4971 & ('ecpso_list_ncent: geom_ncent failed?',911, GEOM_ERR) 4972 num_cent = 0 4973 do ic = 1,ntotal 4974 if (.not.geom_cent_get(geom,ic,gtag,cxyz,chg)) call errquit 4975 & ('ecpso_list_ncent: geom_cent_get failed',911, GEOM_ERR) 4976 inthe_ecp = bas_match_tags(gtag,ecpidin,tag_indx) 4977 inthe_so = bas_match_tags(gtag,soidin,tag_indx) 4978 if (inthe_ecp.or.inthe_so) then 4979 num_cent = num_cent + 1 4980 ecpso_list(num_cent) = ic 4981 endif 4982 enddo 4983 end 4984*---------------------------------------------------------------------- 4985C> \brief Retrieve the list of centers with spin-orbit 4986C> potentials 4987C> 4988C> \return Return .true. if successfull, and .false. otherwise 4989c 4990 logical function so_list_ncent(geom,soidin, 4991 & num_cent,so_list) 4992 implicit none 4993#include "errquit.fh" 4994* 4995* return the number of so centers and the list of centers 4996* 4997*::includes 4998#include "geom.fh" 4999*::functions 5000 logical bas_match_tags 5001 external bas_match_tags 5002*::passed 5003 integer geom !< [Input] the geometry handle 5004 integer soidin !< [Input] the spin-orbit potential (so) handle 5005 integer num_cent !< [Output] the number of so centers 5006 integer so_list(*) !< [Output] the list of so centers 5007*::local 5008 integer ntotal ! total number of centers in goemetry 5009 integer ic ! loop counter for centers 5010 integer tag_indx ! basis set unique center index for tag 5011 character*16 gtag ! geometry tag 5012 double precision cxyz(3) ! coordinates 5013 double precision chg ! charge 5014 logical inthe_so ! is the tag in the so pot.? 5015* 5016 so_list_ncent = geom_ncent(geom,ntotal) 5017 if (.not.so_list_ncent) call errquit 5018 & ('so_list_ncent: geom_ncent failed?',911, GEOM_ERR) 5019 num_cent = 0 5020 do ic = 1,ntotal 5021 if (.not.geom_cent_get(geom,ic,gtag,cxyz,chg)) call errquit 5022 & ('so_list_ncent: geom_cent_get failed',911, GEOM_ERR) 5023 inthe_so = bas_match_tags(gtag,soidin,tag_indx) 5024 if (inthe_so) then 5025 num_cent = num_cent + 1 5026 so_list(num_cent) = ic 5027 endif 5028 enddo 5029 end 5030*---------------------------------------------------------------------- 5031C> \brief Retrieve the list of centers with ECP 5032C> potentials 5033C> 5034C> \return Return .true. if successfull, and .false. otherwise 5035c 5036 logical function ecp_list_ncent(geom,ecpidin, 5037 & num_cent,ecp_list) 5038 implicit none 5039#include "errquit.fh" 5040* 5041* return the number of ecp centers and the list of centers 5042* 5043*::includes 5044#include "geom.fh" 5045*::functions 5046 logical bas_match_tags 5047 external bas_match_tags 5048*::passed 5049 integer geom !< [Input] the geometry handle 5050 integer ecpidin !< [Input] the ECP handle 5051 integer num_cent !< [Output] the number of ECP centers 5052 integer ecp_list(*) !< [Output] the list of ECP centers 5053*::local 5054 integer ntotal ! total number of centers in goemetry 5055 integer ic ! loop counter for centers 5056 integer tag_indx ! basis set unique center index for tag 5057 character*16 gtag ! geometry tag 5058 double precision cxyz(3) ! coordinates 5059 double precision chg ! charge 5060 logical inthe_ecp ! is tag in the ecp? 5061* 5062 ecp_list_ncent = geom_ncent(geom,ntotal) 5063 if (.not.ecp_list_ncent) call errquit 5064 & ('ecp_list_ncent: geom_ncent failed?',911, GEOM_ERR) 5065 num_cent = 0 5066 do ic = 1,ntotal 5067 if (.not.geom_cent_get(geom,ic,gtag,cxyz,chg)) call errquit 5068 & ('ecp_list_ncent: geom_cent_get failed',911, GEOM_ERR) 5069 inthe_ecp = bas_match_tags(gtag,ecpidin,tag_indx) 5070 if (inthe_ecp) then 5071 num_cent = num_cent + 1 5072 ecp_list(num_cent) = ic 5073 endif 5074 enddo 5075 end 5076*---------------------------------------------------------------------- 5077C> \brief Retrieve the coordinates of all centers that have either 5078C> an ECP or spin-orbit potential 5079C> 5080C> \return Return .true. if successfull, and .false. otherwise 5081C 5082 logical function ecpso_get_coords(geom,soidin,ecpidin, 5083 & nxyz,xyz) 5084 implicit none 5085#include "errquit.fh" 5086* 5087* get the coordinates for the ecp/so centers 5088* 5089*::includes 5090#include "geom.fh" 5091*::functions 5092 logical bas_match_tags 5093 external bas_match_tags 5094*::passed 5095 integer geom !< [Input] the geometry handle 5096 integer soidin !< [Input] the spin-orbit potential handle 5097 integer ecpidin !< [Input] the ecp handle 5098 integer nxyz !< [Input] the size of xyz array 5099 double precision xyz(3,nxyz) !< [Output] the coordinates 5100*::local 5101 integer ntotal ! total number of centers in goemetry 5102 integer ic ! loop counter for centers 5103 integer num_cent ! number of centers returned 5104 integer tag_indx ! basis set unique center index for tag 5105 character*16 gtag ! geometry tag 5106 double precision cxyz(3) ! coordinates 5107 double precision chg ! charge 5108 logical inthe_ecp ! is tag in the ecp? 5109 logical inthe_so ! is the tag in the so pot.? 5110* 5111 ecpso_get_coords = geom_ncent(geom,ntotal) 5112 if (.not.ecpso_get_coords) call errquit 5113 & ('ecpso_get_coords: geom_ncent failed?',911, GEOM_ERR) 5114 num_cent = 0 5115 do ic = 1,ntotal 5116 if (.not.geom_cent_get(geom,ic,gtag,cxyz,chg)) call errquit 5117 & ('ecpso_get_coords: geom_cent_get failed',911, GEOM_ERR) 5118 inthe_ecp = bas_match_tags(gtag,ecpidin,tag_indx) 5119 inthe_so = bas_match_tags(gtag,soidin,tag_indx) 5120 if (inthe_ecp.or.inthe_so) then 5121 num_cent = num_cent + 1 5122 if (num_cent.gt.nxyz) call errquit 5123 & ('ecpso_get_coords: array passed in was too small',911, 5124 & GEOM_ERR) 5125 xyz(1,num_cent) = cxyz(1) 5126 xyz(2,num_cent) = cxyz(2) 5127 xyz(3,num_cent) = cxyz(3) 5128 endif 5129 enddo 5130 end 5131*---------------------------------------------------------------------- 5132C> \brief Retrieve the coordinates of all centers that have 5133C> a spin-orbit potential 5134C> 5135C> \return Return .true. if successfull, and .false. otherwise 5136C 5137 logical function so_get_coords(geom,soidin, 5138 & nxyz,xyz) 5139 implicit none 5140#include "errquit.fh" 5141* 5142* get the coordinates for the so centers 5143* 5144*::includes 5145#include "geom.fh" 5146*::functions 5147 logical bas_match_tags 5148 external bas_match_tags 5149*::passed 5150 integer geom !< [Input] the geometry handle 5151 integer soidin !< [Input] the spin-orbit potential handle 5152 integer nxyz !< [Input] the size of xyz array 5153 double precision xyz(3,nxyz) !< [Output] the coordinates 5154*::local 5155 integer ntotal ! total number of centers in goemetry 5156 integer ic ! loop counter for centers 5157 integer num_cent ! number of centers returned 5158 integer tag_indx ! basis set unique center index for tag 5159 character*16 gtag ! geometry tag 5160 double precision cxyz(3) ! coordinates 5161 double precision chg ! charge 5162 logical inthe_so ! is the tag in the so pot.? 5163* 5164 so_get_coords = geom_ncent(geom,ntotal) 5165 if (.not.so_get_coords) call errquit 5166 & ('so_get_coords: geom_ncent failed?',911, GEOM_ERR) 5167 num_cent = 0 5168 do ic = 1,ntotal 5169 if (.not.geom_cent_get(geom,ic,gtag,cxyz,chg)) call errquit 5170 & ('so_get_coords: geom_cent_get failed',911, GEOM_ERR) 5171 inthe_so = bas_match_tags(gtag,soidin,tag_indx) 5172 if (inthe_so) then 5173 num_cent = num_cent + 1 5174 if (num_cent.gt.nxyz) call errquit 5175 & ('so_get_coords: array passed in was too small',911, 5176 & GEOM_ERR) 5177 xyz(1,num_cent) = cxyz(1) 5178 xyz(2,num_cent) = cxyz(2) 5179 xyz(3,num_cent) = cxyz(3) 5180 endif 5181 enddo 5182 end 5183*---------------------------------------------------------------------- 5184C> \brief Retrieve the coordinates of all centers that have 5185C> an ECP 5186C> 5187C> \return Return .true. if successfull, and .false. otherwise 5188C 5189 logical function ecp_get_coords(geom,ecpidin, 5190 & nxyz,xyz) 5191 implicit none 5192#include "errquit.fh" 5193* 5194* get the coordinates for the ecp centers 5195* 5196*::includes 5197#include "geom.fh" 5198*::functions 5199 logical bas_match_tags 5200 external bas_match_tags 5201*::passed 5202 integer geom !< [Input] the geometry handle 5203 integer ecpidin !< [Input] the ecp handle 5204 integer nxyz !< [Input] the size of xyz array 5205 double precision xyz(3,nxyz) !< [Output] the coordinates 5206*::local 5207 integer ntotal ! total number of centers in goemetry 5208 integer ic ! loop counter for centers 5209 integer num_cent ! number of centers returned 5210 integer tag_indx ! basis set unique center index for tag 5211 character*16 gtag ! geometry tag 5212 double precision cxyz(3) ! coordinates 5213 double precision chg ! charge 5214 logical inthe_ecp ! is tag in the ecp? 5215* 5216 ecp_get_coords = geom_ncent(geom,ntotal) 5217 if (.not.ecp_get_coords) call errquit 5218 & ('ecp_get_coords: geom_ncent failed?',911, GEOM_ERR) 5219 num_cent = 0 5220 do ic = 1,ntotal 5221 if (.not.geom_cent_get(geom,ic,gtag,cxyz,chg)) call errquit 5222 & ('ecp_get_coords: geom_cent_get failed',911, GEOM_ERR) 5223 inthe_ecp = bas_match_tags(gtag,ecpidin,tag_indx) 5224 if (inthe_ecp) then 5225 num_cent = num_cent + 1 5226 if (num_cent.gt.nxyz) call errquit 5227 & ('ecp_get_coords: array passed in was too small',911, 5228 & GEOM_ERR) 5229 xyz(1,num_cent) = cxyz(1) 5230 xyz(2,num_cent) = cxyz(2) 5231 xyz(3,num_cent) = cxyz(3) 5232 endif 5233 enddo 5234 end 5235*---------------------------------------------------------------------- 5236C> \brief Check whether a given center from the geometry is known to 5237C> the basis set by that name 5238C> 5239C> This routine looks the tag of the specified center up in the 5240C> geometry and queries the basis set with that tag. 5241C> 5242C> \return Return .true. if the basis set knows of such centers, and 5243C> .false. otherwise 5244C 5245 logical function bas_isce(geom,basisin,center) 5246 implicit none 5247#include "errquit.fh" 5248* 5249* return true if the given center is in the designated basis set 5250* 5251*::includes 5252#include "geom.fh" 5253#include "stdio.fh" 5254*::functions 5255 logical bas_match_tags 5256 external bas_match_tags 5257*::passed 5258 integer geom !< [Input] the geometry handle 5259 integer basisin !< [Input] the basis set handle 5260 integer center !< [Input] the lexical geometry index of center 5261*::local 5262 integer tag_indx ! basis set unique center index for tag 5263 character*16 gtag ! geometry tag 5264 double precision cxyz(3) ! coordinates 5265 double precision chg ! charge 5266* 5267 bas_isce = geom_cent_get(geom,center,gtag,cxyz,chg) 5268 if (.not.bas_isce) call errquit 5269 & ('bas_iscet: geom_cent_get failed',911, GEOM_ERR) 5270* write(LuOut,*)' center = ',center,' gtag:<',gtag,'>' 5271 bas_isce = bas_match_tags(gtag,basisin,tag_indx) 5272 end 5273C 5274C> \brief Retrieve the largest exponent any primitive function in the 5275C> basis set 5276C> 5277C> \return Return .true. if successfull, and .false. otherwise 5278C 5279 logical function bas_large_exponent(basisin,exponent) 5280 implicit none 5281#include "errquit.fh" 5282c 5283c returns the largest primitive exponent in the given basis set 5284c 5285#include "stdio.fh" 5286#include "mafdecls.fh" 5287#include "nwc_const.fh" 5288#include "basP.fh" 5289#include "geobasmapP.fh" 5290#include "basdeclsP.fh" 5291#include "bas_exndcf_dec.fh" 5292#include "bas_ibs_dec.fh" 5293c::functions 5294 logical bas_check_handle 5295 external bas_check_handle 5296c::passed 5297 integer basisin !< [Input] the basis set handle 5298 double precision exponent !< [Output] the maximum exponent 5299c::local 5300 integer basis 5301 integer num_cont 5302 integer icont, ucont 5303 integer myexptr 5304 integer mynprim 5305 integer iprim 5306 double precision test_exponent 5307c 5308#include "bas_exndcf_sfn.fh" 5309#include "bas_ibs_sfn.fh" 5310c 5311 bas_large_exponent = 5312 & bas_check_handle(basisin,'bas_large_exponent') 5313 if (.not.bas_large_exponent) then 5314 write(luout,*)'invalid/inactive basis handle',basisin 5315 call errquit('bas_large_exponent: fatal error',911, BASIS_ERR) 5316 endif 5317 basis = basisin + Basis_Handle_Offset 5318 num_cont = ncont_tot_gb(basis) 5319 exponent = -565.6589d00 5320 do icont = 1,num_cont 5321 ucont = sf_ibs_cn2ucn(icont,basis) 5322 myexptr = infbs_cont(CONT_IEXP,ucont,basis) 5323 mynprim = infbs_cont(CONT_NPRIM,ucont,basis) 5324 do iprim = 1, mynprim 5325 test_exponent = sf_exndcf((myexptr-1+iprim),basis) 5326 exponent = max(exponent,test_exponent) 5327* write(luout,*)' exponent ',exponent,' test ', test_exponent 5328 enddo 5329 enddo 5330* write(luout,*)' largest exponent ',exponent 5331 end 5332C 5333C> \brief Checks whether the basis set contains any generally contracted 5334C> or SP shells 5335C> 5336C> \return Return .true. if any generally contracted or SP shells are 5337C> present, and .false. otherwise 5338C 5339 logical function bas_any_gcorsp(basisin) 5340 implicit none 5341#include "errquit.fh" 5342#include "stdio.fh" 5343#include "nwc_const.fh" 5344#include "basP.fh" 5345c 5346c returns true if the basis set has any sp functions or general contraction 5347c functions. 5348c 5349c::functions 5350 logical bas_check_handle 5351 external bas_check_handle 5352c::passed 5353 integer basisin !< [Input] the basis set handle 5354c::local 5355 integer basis 5356 bas_any_gcorsp = 5357 & bas_check_handle(basisin,'bas_any_gcorsp') 5358 if (.not.bas_any_gcorsp) then 5359 write(luout,*)'invalid/inactive basis handle',basisin 5360 call errquit('bas_any_gcorsp: fatal error',911, BASIS_ERR) 5361 endif 5362 basis = basisin + Basis_Handle_Offset 5363 bas_any_gcorsp = bas_any_gc(basis).or.bas_any_sp_shell(basis) 5364 end 5365C 5366C> \brief Retrieves the largest exponent in the specified shell in the 5367C> basis set 5368C> 5369C> \return Return .true. if successfull, and .false. otherwise 5370C 5371 logical function bas_cont_large_exponent(basisin,incont,exponent) 5372 implicit none 5373#include "errquit.fh" 5374#include "stdio.fh" 5375#include "mafdecls.fh" 5376#include "nwc_const.fh" 5377#include "basP.fh" 5378#include "geobasmapP.fh" 5379#include "basdeclsP.fh" 5380#include "bas_exndcf_dec.fh" 5381#include "bas_ibs_dec.fh" 5382c::functions 5383 logical bas_check_handle 5384 external bas_check_handle 5385c::passed 5386 integer basisin !< [Input] the basis set handle 5387 integer incont !< [Input] the shell rank 5388 double precision exponent !< [Output] the maximum exponent 5389c::local 5390 integer basis 5391 integer num_cont 5392 integer ucont 5393 integer myexptr 5394 integer mynprim 5395 integer iprim 5396 double precision test_exponent 5397c 5398#include "bas_exndcf_sfn.fh" 5399#include "bas_ibs_sfn.fh" 5400c 5401 bas_cont_large_exponent = 5402 & bas_check_handle(basisin,'bas_cont_large_exponent') 5403 if (.not.bas_cont_large_exponent) then 5404 write(luout,*)'invalid/inactive basis handle',basisin 5405 call errquit('bas_cont_large_exponent: fatal error',911, 5406 & BASIS_ERR) 5407 endif 5408 basis = basisin + Basis_Handle_Offset 5409 num_cont = ncont_tot_gb(basis) 5410 exponent = -565.6589d00 5411 ucont = sf_ibs_cn2ucn(incont,basis) 5412 myexptr = infbs_cont(CONT_IEXP,ucont,basis) 5413 mynprim = infbs_cont(CONT_NPRIM,ucont,basis) 5414 do iprim = 1, mynprim 5415 test_exponent = sf_exndcf((myexptr-1+iprim),basis) 5416 exponent = max(exponent,test_exponent) 5417* write(luout,*)' exponent ',exponent,' test ', test_exponent 5418 enddo 5419* write(luout,*)' largest exponent ',exponent 5420 end 5421*..................................................................... 5422C 5423C> \brief Checks whether a relativistic basis set exists on the RTDB 5424C> 5425C> This routine checks for a relativistic basis set by looking for 5426C> the name of the small component basis functions on the RTDB. 5427C> The name of the small component basis is pulled from a common block 5428C> and hence not visible to the caller. 5429C> 5430C> \return Return .true. if the basis set exists on the RTDB, and 5431C> .false. otherwise 5432c 5433 logical function bas_rel_exists (rtdb) 5434 implicit none 5435 integer rtdb !< [Input] the RTDB handle 5436#include "nwc_const.fh" 5437#include "rel_nwc.fh" 5438 logical bas_name_exist_rtdb 5439 external bas_name_exist_rtdb 5440 bas_rel_exists = bas_name_exist_rtdb(rtdb,small_cmpt_name) 5441 & .or. bas_name_exist_rtdb(rtdb,auto_small_cmpt_name) 5442 end 5443*..................................................................... 5444C> \brief Check whether a basis set is a (part of) a relativistic 5445C> basis set 5446C> 5447C> \return Return .true. if the basis set is a relativistic basis set, 5448C> and .false. otherwise 5449C 5450 logical function basis_is_rel (basisin) 5451 implicit none 5452 integer basisin !< [Input] the basis set handle 5453#include "nwc_const.fh" 5454#include "rel_nwc.fh" 5455 character*32 basis_name,trans_name 5456 logical bas_name 5457 external bas_name 5458 basis_is_rel = bas_name(basisin,basis_name,trans_name) 5459 if (basis_is_rel) then 5460 basis_is_rel = (basis_name .eq. small_cmpt_name) 5461 & .or. (basis_name .eq. large_cmpt_name) 5462 & .or. (basis_name .eq. auto_small_cmpt_name) 5463 & .or. (basis_name .eq. auto_large_cmpt_name) 5464 end if 5465 end 5466*..................................................................... 5467C> \brief Retrieves the handle of the basis set with the name 5468C> "ao basis" 5469C> 5470C> \return Return .true. if successfull, and .false. otherwise 5471C 5472 logical function bas_get_ao_handle (basis,nbas,ao_handle) 5473 implicit none 5474 integer nbas !< [Input] the number of basis sets 5475 integer basis(nbas) !< [Input] the array of basis set handles 5476 integer ao_handle !< [Output] the ao basis set handle 5477c 5478c Returns handle of "ao basis" 5479c 5480 integer ibas 5481 logical odum 5482 character*255 basis_name,trans_name 5483* 5484 logical bas_name 5485 external bas_name 5486c 5487 do ibas = 1,nbas 5488 odum = bas_name(basis(ibas),basis_name,trans_name) 5489 if (basis_name(1:8) .eq. 'ao basis') then 5490 bas_get_ao_handle = .true. 5491 ao_handle = basis(ibas) 5492 return 5493 end if 5494 end do 5495 bas_get_ao_handle = .false. 5496c 5497 return 5498 end 5499*..................................................................... 5500C> \brief Checks whether the basis set contains any relativistic shells 5501C> 5502C> \return Return .true. if there are any relativistic shell present, 5503C> and .false. otherwise 5504C 5505 logical function bas_any_rel_shells (basisin) 5506 implicit none 5507#include "errquit.fh" 5508c 5509 integer basisin !< [Input] the basis set handle 5510#include "nwc_const.fh" 5511#include "basdeclsP.fh" 5512#include "basP.fh" 5513c 5514c Checks basis set to see if there are any relativistic shells 5515c 5516 integer ibas ! basis set 5517 integer nucont ! number of unique contractions 5518 integer i, isum 5519c 5520 logical bas_check_handle 5521 external bas_check_handle 5522c 5523 if (.not. bas_check_handle(basisin,'bas_any_rel_shells')) 5524 & call errquit('bas_any_rel_shells: invalid handle',99, 5525 & BASIS_ERR) 5526 ibas = basisin+BASIS_HANDLE_OFFSET 5527 nucont = infbs_head(HEAD_NCONT,ibas) 5528 isum = 0 5529 do i = 1,nucont 5530 isum = isum+infbs_cont(CONT_RELLS,i,ibas) 5531 end do 5532 bas_any_rel_shells = isum .ne. 0 5533 return 5534 end 5535*..................................................................... 5536C> \brief Retrieve the highest angular momentum of all relativistic 5537C> shells in a basis set 5538C> 5539C> \return Return .true. if successfull, and .false. otherwise 5540C 5541 logical function bas_rel_high_ang(basisin,high_angular) 5542 implicit none 5543c 5544c calculate and return highest angular momentem 5545c for relativistic shells in given basis. 5546c 5547#include "stdio.fh" 5548#include "nwc_const.fh" 5549#include "basP.fh" 5550#include "basdeclsP.fh" 5551c::functions 5552 logical bas_check_handle 5553 external bas_check_handle 5554c:: passed 5555 integer basisin !< [Input] the basis set handle 5556 integer high_angular !< [Output] the highest angular momentum 5557c:local 5558 integer basis, i, myang 5559 integer ucont 5560c 5561 bas_rel_high_ang = bas_check_handle(basisin,'bas_high_angular') 5562 if (.not. bas_rel_high_ang ) then 5563 write(LuOut,*) 'bas_rel_high_ang: basis handle not valid ' 5564 return 5565 endif 5566c 5567 basis = basisin + Basis_Handle_Offset 5568 ucont = infbs_head(head_ncont,basis) 5569 high_angular = -565 5570 do i = 1,ucont 5571 if (infbs_cont(cont_rells,i,basis) .ne. 0) then 5572 myang = abs(infbs_cont(cont_type,i,basis)) 5573 high_angular = max(high_angular,myang) 5574 end if 5575 enddo 5576c 5577 bas_rel_high_ang = .true. 5578 return 5579 end 5580C 5581C> \brief Checks whether two tags match 5582C> 5583C> This routine matches two input tags using the geometry/basis set 5584C> tags matching mechanism, in such a way that "H34" matches "H" or 5585C> "hydrogen" or "h" if "H34" is one of the tags. "H" will however 5586C> not match "h" i.e., case sensitivity is preserved. 5587C> 5588C> \return Return .true. if the tags match, and .false. otherwise 5589C 5590 logical function bas_do_tags_match(tag_one,tag_two) 5591 implicit none 5592#include "errquit.fh" 5593* 5594* This routine matches two input tags using the geometry/basis set 5595* tags matching mechanism, in such a way that "H34" matches "H" or 5596* "hydrogen" or "h" if "H34" is one of the tags. "H" will however 5597* not match "h" i.e., case sensitivity is preserved. 5598* 5599c::includes 5600#include "stdio.fh" 5601#include "inp.fh" 5602#include "nwc_const.fh" 5603c::functions 5604 logical geom_tag_to_element 5605 external geom_tag_to_element 5606c::passed 5607 character*(*) tag_one !< [Input] a geometry/basis tag 5608 character*(*) tag_two !< [Input] another geometry/basis tag 5609* bas_do_tags_match ! [output] true if they match 5610c::local 5611 character*16 one16 ! geometry or basis tags can only be *16 5612 character*16 two16 ! geometry or basis tags can only be *16 5613 integer lnotalpha ! string index to non alpha character 5614 integer i ! loop counter 5615 integer ind ! dummy counter 5616 logical status ! dummy storage 5617 integer len_one ! length of tag one or substring of tag 5618 integer len_one_old ! length of tag one 5619 integer len_two ! length of tag two or substring of tag 5620 integer len_two_old ! length of tag two 5621 character*2 one_sym ! tag one -> symbol name 5622 character*16 one_elem ! tag one -> element name 5623 integer one_atn ! tag one -> atomic number 5624 character*2 two_sym ! tag two -> symbol name 5625 character*16 two_elem ! tag two -> element name 5626 integer two_atn ! tag two -> atomic number 5627 logical debug ! true for extra output 5628c 5629 integer na2z 5630 parameter (na2z = 26) 5631 character*1 a2z(na2z) 5632 data a2z /'a','b','c','d','e','f','g','h','i','j', 5633 & 'k','l','m','n','o','p','q','r','s','t', 5634 & 'u','v','w','x','y','z'/ 5635c 5636 debug = .true. 5637c 5638 bas_do_tags_match = .false. 5639c 5640 one16 = tag_one ! copy tag to work on it 5641 two16 = tag_two ! copy tag to work on it 5642c 5643* remove bq ghost info if it exists 5644c 5645 if (one16(1:4).eq.'bqgh') then 5646 one16 = one16(5:) 5647 endif 5648 if (two16(1:4).eq.'bqgh') then 5649 two16 = two16(5:) 5650 endif 5651c 5652c... first match full geometry tag to full basis tag list 5653c did the user specifically assign a tag ? 5654 5655 if (one16.eq.two16) then 5656 bas_do_tags_match = .true. 5657 goto 00009 5658 endif 5659c 5660c... Now check to see if either tag has non alpha chracters. 5661c... the substring prior to the first non-alpha character is the users idea 5662c... of the "name" of the tag. 5663c 5664 len_one = inp_strlen(one16) 5665 len_one_old = len_one 5666 lnotalpha = len_one+1 5667 do i = 1,len_one 5668 if (.not.(inp_match(na2z,.false., 5669 & one16(i:i),a2z,ind))) then ! compare character to alpha (case-less test) 5670 lnotalpha = i 5671 goto 00001 5672 endif 5673 enddo 567400001 continue 5675 do i = lnotalpha,len_one 5676 one16(i:i) = ' ' 5677 enddo 5678 if (debug) write(luout,*)' the string:', tag_one, 5679 & 'has the user substring of ',one16 5680 5681 len_two = inp_strlen(two16) 5682 len_two_old = len_two 5683 lnotalpha = len_two+1 5684 do i = 1,len_two 5685 if (.not.(inp_match(na2z,.false., 5686 & two16(i:i),a2z,ind))) then ! compare character to alpha (case-less test) 5687 lnotalpha = i 5688 goto 00002 5689 endif 5690 enddo 569100002 continue 5692 do i = lnotalpha,len_two 5693 two16(i:i) = ' ' 5694 enddo 5695 if (debug) write(luout,*)' the string:', tag_two, 5696 & 'has the user substring of ',two16 5697c 5698c... match user substrings 5699 len_one = inp_strlen(one16) 5700 len_two = inp_strlen(two16) 5701 if ((len_one_old.gt.len_one).or.(len_two_old.gt.len_two)) then 5702 if (one16.eq.two16) then 5703 bas_do_tags_match = .true. 5704 goto 00009 5705 endif 5706 endif 5707c 5708c... now get symbol and element names for each tag and match those 5709c use the substring for each tag to match these 5710c... bq's with basis functions must match from the above rules! 5711c 5712 status = geom_tag_to_element(one16,one_sym,one_elem,one_atn) 5713 if (.not.status)then 5714 if (.not.(one_sym.eq.'bq')) then 5715 write(luout,*)'tag <',one16, 5716 & '> could not be matched to an element symbol' 5717 call errquit('bas_do_tags_match: fatal error ',911, BASIS_ERR) 5718 endif 5719 endif 5720 if (one_sym.eq.'bq') goto 00009 ! bq labels with basis functions must match from above 5721 status = geom_tag_to_element(two16,two_sym,two_elem,two_atn) 5722 if (.not.status)then 5723 if (.not.(two_sym.eq.'bq')) then 5724 write(luout,*)'tag <',two16, 5725 & '> could not be matched to an element symbol' 5726 call errquit('bas_do_tags_match: fatal error ',911, BASIS_ERR) 5727 endif 5728 endif 5729 if (two_sym.eq.'bq') goto 00009 ! bq labels with basis functions must match from above 5730* 5731 if (one_elem.eq.two_elem) then 5732 bas_do_tags_match = .true. 5733 goto 00009 5734 else if (one_sym.eq.two_sym) then 5735 bas_do_tags_match = .true. 5736 goto 00009 5737 endif 5738 if (debug) 5739 & write(luout,*)'bas_do_tags_match:debug: no match for tags <', 5740 & tag_one,'> and <',tag_two,'>' 5741 return ! no match 5742c 574300009 continue 5744 if (debug) then 5745 write(luout,10000)tag_one,tag_two 5746 endif 5747 return 574810000 format('bas_do_tags_match:debug: tag ',a16, 5749 & ' matched this tag ',a16) 5750 end 5751C> @} 5752