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