1 subroutine texas_hf2_m(ij_basis,ish,jsh,kl_basis,ksh,lsh,nquart, 2 * q4, use_q4, 3 * ra,rb,rc,rd,use_r, 4 * eri,leri,icf,jcf,kcf,lcf,integ_n0,labels, 5 * more_int,blscr,l_blscr,zerotol,int_type) 6* 7* $Id$ 8* 9 implicit none 10#include "mafdecls.fh" 11#include "bas.fh" 12c 13 character*8 int_type 14 integer ij_basis, ish(*), jsh(*), kl_basis, ksh(*), lsh(*) 15 integer nquart 16 double precision q4(*) 17 logical use_q4 18 double precision ra(*), rb(*), rc(*), rd(*) 19 logical use_r 20 double precision eri(*) 21 integer leri 22 integer icf(*), jcf(*), kcf(*), lcf(*) 23 integer integ_n0 24 logical labels 25 logical more_int 26 double precision blscr(*) 27 integer l_blscr 28 double precision zerotol 29c 30 double precision time_beg,time_end 31 double precision time4texas_hf2_m,time4mul_quart 32 common /pnl_time/ time4texas_hf2_m,time4mul_quart 33c------------------------------------------------------------------ 34 call txs_second(time_beg) 35c------------------------------------------------------------------ 36c send info about type of integrals requested : 37c 38 call requested_task(int_type) 39c------------------------------------------------------------------ 40 call texas_hf2_m2(ij_basis,ish,jsh,kl_basis,ksh,lsh,nquart, 41 * q4, use_q4, 42 * ra,rb,rc,rd,use_r, 43 * eri,leri,icf,jcf,kcf,lcf,integ_n0,labels, 44 * more_int,blscr,l_blscr,zerotol) 45c------------------------------------------------------------------ 46 call txs_second(time_end) 47 time4texas_hf2_m=time4texas_hf2_m + (time_end-time_beg) 48c------------------------------------------------------------------ 49 end 50c================================================================= 51 subroutine texas_hf2_m2(ij_basis,ish,jsh,kl_basis,ksh,lsh,nquart, 52 * q4, use_q4, 53 * ra,rb,rc,rd,use_r, 54 * eri,leri,icf,jcf,kcf,lcf,integ_n0,labels, 55 * more_int,blscr,l_blscr,zerotol) 56c-------------------------------------------------------- 57c This subroutine delivers two-el.integrals 58c 59c (1) all of them (zeros also) WITHOUT labels 60c (2) only non-zero integrals WITH labels 61c 62c NQUART quartets of contr. shells is requested 63c ish(nquart), ..,lsh(nquart) 64c Integrals return in ERI(leri) without labels or with labels 65c in icf(leri),..,lcf(leri) . 66c-------------------------------------------------------- 67 implicit none 68#include "errquit.fh" 69#include "mafdecls.fh" 70 logical use_r,labels, more_int 71c PNL scratch for calculations replacing our bl() and its size : 72 integer l_blscr 73 double precision blscr(l_blscr) 74c PNL requested set of contracted shell quartets : 75 integer ij_basis, kl_basis 76 integer nquart,ish(nquart),jsh(nquart),ksh(nquart),lsh(nquart) 77 logical use_q4 78 double precision q4(nquart) 79c returning integral's indeces : 80 integer icf(*),jcf(*),kcf(*),lcf(*), integ_n0 81 integer leri 82 double precision eri(leri) 83 double precision ra(3),rb(3),rc(3),rd(3) 84c screening threshold for output integrals 85 double precision zerotol 86c---------------------------------------------------------------------- 87 double precision savezerotol 88 common /csavezerotol/ savezerotol ! Used in detbul 89c---------------------------------------------------------------------- 90 integer ntxs_bl_scr 91 common /bl_txs_add/ ntxs_bl_scr 92c---------------------------------------------------------------------- 93 integer basis_init 94 common /c_basis_init/ basis_init 95c---------------------------------------------------------------------- 96 integer num_bas_1,num_bas_2,num_bas_3 97 integer ncs_bas_1,ncs_bas_2,ncs_bas_3 98 integer nps_bas_1,nps_bas_2,nps_bas_3 99 integer nat_bas_1,nat_bas_2,nat_bas_3 100 integer ncf_bas_1,ncf_bas_2,ncf_bas_3 101 common /multi_basis/ num_bas_1,num_bas_2,num_bas_3, ! Basis set handle 102 * ncs_bas_1,ncs_bas_2,ncs_bas_3, ! Cummulative #shells 103 * nps_bas_1,nps_bas_2,nps_bas_3, ! Cummulative #prims unused 104 * nat_bas_1,nat_bas_2,nat_bas_3, ! Cummulative #atoms unused 105 * ncf_bas_1,ncf_bas_2,ncf_bas_3 ! Cummulative #basis functions 106c 107c to see if any of these first shells from this request are ZERO 108c 109 integer ish_first,jsh_first,ksh_first,lsh_first 110c---------------------------------------------------------------------- 111c... check for consistency the input basis is what was initialized. 112c 113 savezerotol = zerotol 114c 115 if(.not.more_int) then 116c ---------------------------------------------- 117 if(ij_basis.eq.num_bas_1) go to 1001 118 if(ij_basis.eq.num_bas_2) go to 1001 119 if(ij_basis.eq.num_bas_3) go to 1001 120c 121 write(6,*)' basis sets initialized :' 122 if(num_bas_1.gt.0) write(6,*) num_bas_1 123 if(num_bas_2.gt.0) write(6,*) num_bas_2 124 if(num_bas_3.gt.0) write(6,*) num_bas_3 125 write(6,*)' ij_basis handle :',ij_basis 126 call errquit 127 & ('texas_hf2_m: called with different basis set handle',911, 128 & BASIS_ERR) 129 1001 continue 130c ---------------------------------------------- 131 if(kl_basis.eq.num_bas_1) go to 1002 132 if(kl_basis.eq.num_bas_2) go to 1002 133 if(kl_basis.eq.num_bas_3) go to 1002 134c 135 write(6,*)' basis sets initialized :' 136 if(num_bas_1.gt.0) write(6,*) num_bas_1 137 if(num_bas_2.gt.0) write(6,*) num_bas_2 138 if(num_bas_3.gt.0) write(6,*) num_bas_3 139 write(6,*)' kl_basis handle :',kl_basis 140 call errquit 141 & ('texas_hf2_m: called with different basis set handle',911, 142 & BASIS_ERR) 143 1002 continue 144c -------------------------------------------------------------- 145c remember the first shells in this request(if ZERO or not) 146 ish_first=ish(1) 147 jsh_first=jsh(1) 148 ksh_first=ksh(1) 149 lsh_first=lsh(1) 150c -------------------------------------------------------------- 151c Check what basis sets are involved and what type of integrals 152c is requested (4-c , 3-c or 2-c ; c=center) 153c 154 call request_update(ish,jsh,ksh,lsh,nquart,ij_basis,kl_basis ) 155c -------------------------------------------------------------- 156c switch from txs bl() to pnl blscr() 157c 158 call switch_scr(dbl_mb(ntxs_bl_scr),blscr,l_blscr) 159c -------------------------------------------------------------- 160 endif 161c---------------------------------------------------------------------- 162 call mul_quart(ish,jsh,ksh,lsh,nquart, q4,use_q4, 163 * ra,rb,rc,rd,use_r, blscr,l_blscr,eri,leri, 164 * icf,jcf,kcf,lcf,integ_n0,labels, more_int) 165c---------------------------------------------------------------------- 166c In the case of different basis sets update the label's arrays 167c 168 if(labels) then 169 call labels_update(icf,jcf,kcf,lcf,integ_n0,ij_basis,kl_basis, 170 * ish_first,jsh_first,ksh_first,lsh_first) 171 endif 172c---------------------------------------------------------------------- 173 end 174c================================================================= 175 subroutine mul_quart(icspnl,jcspnl,kcspnl,lcspnl,nquart,q4,use_q4, 176 * ra,rb,rc,rd,use_r, bl ,l_blscr,eri,leri, 177 * icf,jcf,kcf,lcf,integ_nx,labels, more_int) 178 implicit real*8 (a-h,o-z) 179#include "mafdecls.fh" 180 logical use_q4,use_r,labels,more_int 181 double precision ra(3), rb(3), rc(3), rd(3) 182 common /ctxs_index/ maxsh,ifp,inx(1) 183 common /ganz/ lcore,iov,last,lflag(4),inuc,ibas,na,nbf,nsh,ncf,ncs 184 1,nsy(4),nsym,nganz(35),lopt(30) 185c---------------------------------------------------------------- 186cccc common /memmax/ ispblx, maxme1,iforwhat 187 common /mem_max_min/ ispblx,maxme1,max_111,iforwhat 188c---------------------------------------------------------------- 189 common /pnl001/ ispec,ijpres2,klpres2, ijblock,klblock,iqorder 190 common /pnl002/ ncshell,ncfunct,nblock2,integ_n0 191c---------------------------------------------------------------- 192 common /pnl006/ nsplit,isplit,isbl_split, isbl_part 193c---------------------------------------------------------------- 194 common /pnl_time/ time4texas_hf2_m,time4mul_quart 195c---------------------------------------------------------------- 196 common /pnl_nqrt/ ncall_pnl,ncall_41,nqrts_pnl,npart_pnl,nsize_pnl 197c---------------------------------------------------------------- 198cpnl common /big/ bl(1) 199 dimension bl(l_blscr), q4(nquart) 200c 201c requested quartets of shells : 202 dimension icspnl(nquart),jcspnl(nquart), 203 * kcspnl(nquart),lcspnl(nquart) 204c returning integrals and indeces : 205 dimension eri(leri) 206 dimension icf(*),jcf(*),kcf(*),lcf(*) 207c---------------------------------------------------------------- 208c eri - PNL buffer for returning integrals 209c---------------------------------------------------------------- 210 call txs_second(time_beg) 211c---------------------------------------------------------------- 212c memory checking ; 213c 214c call getmem(0,last11) 215c call retmem(1) 216c---------------------------------------------------------------- 217c MORE_INT=true means that it is called for the same request again 218c 219 999 continue 220c 221 if(more_int) then 222 nsplit=nsplit+1 223 go to 1000 224 endif 225c---------------------------------------------------------------- 226 ncall_pnl=ncall_pnl+1 227 nqrts_pnl=nqrts_pnl+nquart 228 if(nquart.eq.1) ncall_41=ncall_41 + 1 229ctest 230c write(6,*) 231c * ' call no=',ncall_pnl,' nquart=',nquart,' nq_tot=',nqrts_pnl 232c---------------------------------------------------------------- 233c if use_r is true replace current coordinates by Ra-Rd 234c 235c if (use_r) then 236c call replacer(bl(inuc+1),ra,rb,rc,rd,inx,nquart, 237c * icspnl,jcspnl,kcspnl,lcspnl,'forw') 238c endif 239c---------------------------------------------------------------- 240c find texas-shells corresponding to the given set of pnl-shells. 241c re-oredr TXS-requested shells in quartets to TXS-work ordering : 242c find wihch TXS pairs are present in a request (ijpres2,klpres2) 243c find all pair-blocks to which requested shells belong : 244c 245 call txs_setup(bl,inx,ncs,leri, 246 * nquart,icspnl,jcspnl,kcspnl,lcspnl,labels) 247ctest 248c call texas_terminate() 249c stop 'stopped by kw after txs_setup' 250ctest 251c 252c 11 calls to getmem 253c 12 calls to getmem 254c if(.not.labels) +1 call to getmem 255c---------------------------------------------------------------- 256c if this pnl-request is too big (too many integrals at once) 257c 258 if( isplit.gt.0 ) then 259 more_int=.true. 260 nsplit=0 261 go to 999 262 endif 263c---------------------------------------------------------------- 264 1000 continue 265c---------------------------------------------------------------- 266c calculate two-electron integrals : 267c 268 npart_pnl=npart_pnl+1 269c 270 if(.not.labels) then 271 ispec=1 272 integ_n0=0 273c 274c calculate integrals without labels : 275c 276c ncfunct is not used on this call 277 call calcint2(bl,inx, q4,use_q4, icf,jcf,kcf,lcf,eri, 278 * dbl_mb(ncfunct)) 279 integ_nx=integ_n0 280c 281 else 282c 283 ispec=2 284 integ_n0=0 285c 286c calculate integrals with labels : 287c 288 call calcint2(bl,inx, q4,use_q4, icf,jcf,kcf,lcf,eri, 289 * dbl_mb(ncfunct)) 290 integ_nx=integ_n0 291c 292 endif 293c 294 if(isplit.gt.0) then 295 if(nsplit.eq.isplit+1) more_int=.false. 296 endif 297c 298 if(.not.more_int) then 299 if(.not.labels) then 300c10 call retmem(13) 301 call retmem(11) 302 else 303c10 call retmem(12) 304 call retmem(10) 305 endif 306c---------------------------------------------------------------- 307c if use_r is true return original coordinates 308c Probably NOT needed ! 309c--- 310c if (use_r) then 311c call replacer(bl(inuc+1),ra,rb,rc,rd,inx,nquart, 312c * icspnl,jcspnl,kcspnl,lcspnl,'back') 313c endif 314c---------------------------------------------------------------- 315c check and update memory status in pnl_scratch 316c after each pnl request is done; 317c 318 call memstat_pnl 319c 320 endif 321c---------------------------------------------------------------- 322 call txs_second(time_end) 323 time4mul_quart=time4mul_quart + (time_end -time_beg) 324c---------------------------------------------------------------- 325 if(integ_nx.eq.0 .and. more_int) go to 999 326c---------------------------------------------------------------- 327c memory checking ; 328c 329c call getmem(0,last12) 330c call retmem(1) 331c if(last11.ne.last12) then 332c write(6,*)'** Memory alocations in MUL_QUART **' 333c write(6,*)' ncall_pnl=',ncall_pnl,'more_int=',more_int 334c write(6,*)' isplit=',isplit,' nsplit=',nsplit 335c write(6,*)' at the beginning of calc.=',last11 336c write(6,*)' at the end of cal. =',last12 337c endif 338c---------------------------------------------------------------- 339c 340 end 341c================================================================== 342 subroutine replacer(datnuc,ra,rb,rc,rd,inx,nquart, 343 * icspnl,jcspnl,kcspnl,lcspnl, direction) 344 implicit real*8 (a-h,o-z) 345 character*4 direction 346 common /keepr/ ri(3),rj(3),rk(3),rl(3) 347 dimension ra(3),rb(3),rc(3),rd(3) 348 dimension datnuc(5,*) 349 dimension inx(12,*) 350 dimension icspnl(*),jcspnl(*),kcspnl(*),lcspnl(*) ! dim=nquart 351c 352 if(direction.eq.'forw') then 353 do 100 iq=1,nquart 354 ics=icspnl(iq) 355 jcs=jcspnl(iq) 356 kcs=kcspnl(iq) 357 lcs=lcspnl(iq) 358c 359 iat=inx(2,ics) 360 jat=inx(2,jcs) 361 kat=inx(2,kcs) 362 lat=inx(2,lcs) 363c 364 do 10 ii=1,3 365 if(iat.gt.0) ri(ii)=datnuc(ii+1,iat) 366 if(jat.gt.0) rj(ii)=datnuc(ii+1,jat) 367 if(kat.gt.0) rk(ii)=datnuc(ii+1,kat) 368 if(lat.gt.0) rl(ii)=datnuc(ii+1,lat) 369 10 continue 370 do 20 ii=1,3 371 if(iat.gt.0) datnuc(ii+1,iat)=ra(ii) 372 if(jat.gt.0) datnuc(ii+1,jat)=rb(ii) 373 if(kat.gt.0) datnuc(ii+1,kat)=rc(ii) 374 if(lat.gt.0) datnuc(ii+1,lat)=rd(ii) 375 20 continue 376 100 continue 377 endif 378c 379 if(direction.eq.'back') then 380 do 200 iq=1,nquart 381 ics=icspnl(iq) 382 jcs=jcspnl(iq) 383 kcs=kcspnl(iq) 384 lcs=lcspnl(iq) 385c 386 iat=inx(2,ics) 387 jat=inx(2,jcs) 388 kat=inx(2,kcs) 389 lat=inx(2,lcs) 390 do 30 ii=1,3 391 if(iat.gt.0) datnuc(ii+1,iat)=ri(ii) 392 if(jat.gt.0) datnuc(ii+1,jat)=rj(ii) 393 if(kat.gt.0) datnuc(ii+1,kat)=rk(ii) 394 if(lat.gt.0) datnuc(ii+1,lat)=rl(ii) 395 30 continue 396 200 continue 397 endif 398c 399 end 400c================================================================== 401 subroutine txs_setup(bl,inx,ncs,leri, 402 * nquart,icspnl,jcspnl,kcspnl,lcspnl,labels) 403 implicit real*8 (a-h,o-z) 404#include "mafdecls.fh" 405 logical labels 406c 407 common /bl_txs_add/ ntxs_bl_scr 408c 409 common /memor1/ npard,mxsize,nblock1,nblock1_back 410 common /memor11/ mxpair 411 common /memor1b/ nbl2 412c 413 common /pnl001/ ispec,ijpres2,klpres2, ijblock,klblock,iqorder 414 common /pnl002/ ncshell,ncfunct,nblock2,integ_n0 415 common /pnl003/ nqrtpnl,icstxs,jcstxs,kcstxs,lcstxs 416 common /pnl004/ isize,jsize,ksize,lsize,itxspnl 417 common /pnl005/ isblsize,isblqrts,isblpoint 418 common /pnl006/ nsplit,isplit,isbl_split, isbl_part 419 dimension icspnl(*),jcspnl(*),kcspnl(*),lcspnl(*) ! dim=nquart 420 dimension inx(12,*) 421 dimension bl(*) 422c-------------------------------------------------------- 423 nqrtpnl=nquart 424 nsupblk=nbl2*(nbl2+1)/2 425c-------------------------------------------------------- 426c find texas-shells corresponding to the given set of pnl-shells. 427c 428 call memo1_int(nquart,icstxs) 429 call memo1_int(nquart,jcstxs) 430 call memo1_int(nquart,kcstxs) 431 call memo1_int(nquart,lcstxs) 432c 433C@@@@ call txs_pnl(nquart,bl(ncshell), 434 call txs_pnl(nquart,dbl_mb(ncshell), 435 * bl(icstxs),bl(jcstxs),bl(kcstxs),bl(lcstxs), 436 * icspnl,jcspnl,kcspnl,lcspnl) 437c 438c on return TXS shells in icspnl... again 439c 440 call retmem(4) 441c 442c-------------------------------------------------------- 443c reserve memory 444c for pairs ij, kl which are present in a request 445c 446 call memo1_int(nquart, ijpres2) 447 call memo1_int(nquart, klpres2) 448c- 449 call memo1_int(nquart, iqorder) 450c-------------------------------------------------------- 451c reserve memory 452c for pair-blocks ij, kl which are present in a request 453c 454 call memo1_int(nbl2, ijblock) 455 call memo1_int(nbl2, klblock) 456c-------------------------------------------------------- 457c reserve memory for present super-blocks 458c 459 call memo1_int(nquart,isblqrts) 460 call memo1_int(nsupblk,isblsize) 461c-------------------------------------------------------- 462c 7 calls to getmem 463c-------------------------------------------------------- 464 call memo1_int(nquart,isupblk) 465 call memo1_int(nsupblk,isblpoin) 466* call txs_second(time_beg) 467 call swap_shells(nquart,icspnl,jcspnl,kcspnl,lcspnl,nbl2, 468 * bl(ijblock),bl(klblock),bl(ijpres2),bl(klpres2), 469 * bl(isupblk),bl(isblsize),bl(isblqrts),bl(nblock1_back), 470 * bl(iqorder),bl(isblpoin) ) 471 call txs_second(time_end) 472* time_swap=time_swap +(time_end-time_beg) 473 call retmem(1) 474 call retmem(1) 475c 476 if(.not.labels) then 477 call get_sizes(inx,icspnl,jcspnl,kcspnl,lcspnl, 478 * isize,jsize,ksize,lsize,ngcont) 479ccc write(6,*)'general contraction deep=',ngcont 480c 481 call memo1_int(ngcont*isize*jsize*ksize*lsize,itxspnl) 482c 483 call ordr_shells(icspnl,jcspnl,kcspnl,lcspnl, 484 * isize,jsize,ksize,lsize,bl(itxspnl), 485 * bl(iqorder),ngcont) 486c 1 call to getmem 487 endif 488c--------------------------------------------------------- 489c calculate the pointer to the block's quartets 490c (to eliminate calcul. in the get_nbls_pnl routine) 491c 492 call memo1_int(nsupblk,isblpoint) 493 call supblks_point(nsupblk,bl(isblsize),bl(isblpoint)) 494c-------------------------------------------------------- 495c check current sizes of super-blocks against maximum sizes 496c and available size of integral's buffer (leri): 497c (maximum sizes calculated in blksizer (spec_block.f) ) 498c 499 call memo1_int(nsupblk,isbl_part) 500 call supblks_split(bl,nquart,bl(mxsize),bl(mxpair), 501 * bl(ijblock),bl(klblock), 502 * leri,inx, bl(nblock1),nbl2, 503 * bl(isblqrts),bl(isblpoint), 504 * bl(ijpres2),bl(klpres2),bl(isblsize),bl(isbl_part)) 505c--------------------------------------------------------- 506c calculate the size of the request and split it if needed 507c 508 call memo1_int(nsupblk+2,isbl_split) 509 call request_split(leri,inx, nbl2,bl(ijblock),bl(klblock), 510 * bl(ijpres2),bl(klpres2),bl(isblsize),bl(isblqrts), 511 * bl(isbl_part), 512 * bl(isbl_split),isplit ) 513c 514 call retmem(1) 515 call memo1_int(isplit+2, isbl_split) 516c-------------------------------------------------------- 517c2002 : efficeiency of using blocking capabilities of Texas: 518c 519 nquartot=ncs*(ncs+1)/2 520 nquartot=nquartot*(nquartot+1)/2 521 call blk_capab(nbl2,bl(npard),bl(isblsize),nquart,nquartot) 522c-------------------------------------------------------- 523 end 524c================================================================== 525 subroutine get_sizes(inx, icstxs,jcstxs,kcstxs,lcstxs, 526 * isize,jsize,ksize,lsize,ngcont) 527 dimension inx(12,*),icstxs(*),jcstxs(*), kcstxs(*),lcstxs(*) 528c 529 ics=icstxs(1) 530 jcs=jcstxs(1) 531 ijcs=(ics-1)*ics/2 +jcs 532 kcs=kcstxs(1) 533 lcs=lcstxs(1) 534 klcs=(kcs-1)*kcs/2 +lcs 535c 536 igcon=inx(4, ics) + 1 537 jgcon=inx(4, jcs) + 1 538 kgcon=inx(4, kcs) + 1 539 lgcon=inx(4, lcs) + 1 540c 541 ngcont=igcon*jgcon*kgcon*lgcon 542c 543 isize=inx(3, ics) 544 jsize=inx(3, jcs) 545 ksize=inx(3, kcs) 546 lsize=inx(3, lcs) 547c 548 end 549c================================================================== 550 subroutine swap_shells(nquart,icstxs,jcstxs,kcstxs,lcstxs, 551 * nbl2, 552 * ijblock,klblock,ijpres2,klpres2, 553 * isupblk, isbl_s,isbl_q, nblock1_back, 554 * iqorder,isbl_p) 555 dimension icstxs(nquart),jcstxs(nquart), 556 * kcstxs(nquart),lcstxs(nquart) 557 dimension ijblock(nbl2),klblock(nbl2) 558 dimension ijpres2(nquart),klpres2(nquart) 559 dimension iqorder(nquart) 560 dimension isupblk(nquart) 561 dimension isbl_s(*) ! dimension nbl2*(nbl2+1)/2 num of super-blks 562 dimension isbl_q(nquart) 563 dimension nblock1_back(*) ! dimension ncs 564 dimension ipnl(4) 565 dimension isbl_p(*) ! dimension nbl2*(nbl2+1)/2 num of super-blks 566c 567c 568 isbl=0 569 do 05 ibl=1,nbl2 570 ijblock(ibl)=0 571 klblock(ibl)=0 572 do 05 kbl=1,ibl 573 isbl=isbl+1 574 isbl_s(isbl)=0 575 05 continue 576c 577 do 10 iq=1,nquart 578 ipnl(1)=1 579 ipnl(2)=2 580 ipnl(3)=3 581 ipnl(4)=4 582c 583 ics=icstxs(iq) 584 jcs=jcstxs(iq) 585 kcs=kcstxs(iq) 586 lcs=lcstxs(iq) 587c 588 if(ics.lt.jcs) then 589 ii=ics 590 ics=jcs 591 jcs=ii 592 ipnl(1)=2 593 ipnl(2)=1 594 endif 595 if(kcs.lt.lcs) then 596 kk=kcs 597 kcs=lcs 598 lcs=kk 599 ipnl(3)=4 600 ipnl(4)=3 601 endif 602c 603c now find ijcs and klcs pairs and find to which 604c pair-blocks they belong : 605c 606 ijcs=ics*(ics-1)/2+jcs 607 klcs=kcs*(kcs-1)/2+lcs 608cnew 609 i_block1=nblock1_back(ics) 610 j_block1=nblock1_back(jcs) 611 ijblock1=i_block1*(i_block1-1)/2 + j_block1 612 k_block1=nblock1_back(kcs) 613 l_block1=nblock1_back(lcs) 614 klblock1=k_block1*(k_block1-1)/2 + l_block1 615cnew 616c 617 iswitch=0 618 if(ijblock1.lt.klblock1) iswitch=1 619 if(ijblock1.eq.klblock1 .and. ijcs.lt.klcs) iswitch=1 620c 621 if( iswitch.eq.1 ) then 622 ii=ics 623 ics=kcs 624 kcs=ii 625 jj=jcs 626 jcs=lcs 627 lcs=jj 628 ij=ijcs 629 ijcs=klcs 630 klcs=ij 631 i1=ipnl(1) 632 j1=ipnl(2) 633 ipnl(1)=ipnl(3) 634 ipnl(2)=ipnl(4) 635 ipnl(3)=i1 636 ipnl(4)=j1 637 ijblk=ijblock1 638 ijblock1=klblock1 639 klblock1=ijblk 640 endif 641c 642 ijblock(ijblock1)=1 643 klblock(klblock1)=1 644c 645 icstxs(iq)=ics 646 jcstxs(iq)=jcs 647 kcstxs(iq)=kcs 648 lcstxs(iq)=lcs 649c 650 ijpres2(iq)=ijcs 651 klpres2(iq)=klcs 652c 653 isupblk(iq)=ijblock1*(ijblock1-1)/2+klblock1 654c 655 iorder=1000*ipnl(1)+100*ipnl(2)+10*ipnl(3)+ipnl(4) 656 iqorder(iq)=iorder 657c 658 10 continue 659c 660c--------------------------------------------------------- 661c calculate how many super-blocks are matched by a given 662c PNL request and to what degree : 663c 664c 665 do 30 iq=1,nquart 666 isblq=isupblk(iq) 667 isbl_s(isblq)=isbl_s(isblq)+1 668 30 continue 669c 670 ipoint=0 671 isbl_p(1)=0 672 do 40 ikbl=2,nbl2*(nbl2+1)/2 673 ipoint=ipoint+isbl_s(ikbl-1) 674 isbl_p(ikbl)=ipoint 675 40 continue 676c 677 do 50 iq=1,nquart 678 isblq=isupblk(iq) 679 isbl_p(isblq)=isbl_p(isblq)+1 680 iquart=isbl_p(isblq) 681 isbl_q(iquart)=iq 682 50 continue 683c--------------------------------------------------------- 684c 685 end 686c================================================================== 687 subroutine ordr_shells(icstxs,jcstxs,kcstxs,lcstxs, 688 * isize,jsize,ksize,lsize,itxspnl, 689 * iqorder,ngcq) 690 dimension icstxs(*),jcstxs(*),kcstxs(*),lcstxs(*) 691 dimension iqorder(*),itxspnl(*) 692c 693c all dimensions are = nquart 694c--------------------------------------------------------- 695c Only for a run without labels : 696c 697 iorder=iqorder(1) 698c 699 if(iorder.eq.1234) then 700c ijkl->ijkl 701 ls=0 702 ks=lsize 703 js=ksize*ks 704 is=jsize*js 705 im=0 706 jm=0 707 km=0 708 lm=1 709 call txspnl(is,im, js,jm, ks,km ,ls,lm, itxspnl,ngcq) 710 endif 711 if(iorder.eq.1243) then 712c ijkl->ijlk 713 ls=ksize 714 ks=0 715 js=lsize*ls 716 is=jsize*js 717 im=0 718 jm=0 719 km=1 720 lm=0 721 call txspnl(is,im, js,jm, ks,km ,ls,lm, itxspnl,ngcq) 722 endif 723 if(iorder.eq.2134) then 724c ijkl->jikl 725 ls=0 726 ks=lsize 727 is=ksize*ks 728 js=isize*is 729 im=0 730 jm=0 731 km=0 732 lm=1 733 call txspnl(is,im, js,jm, ks,km ,ls,lm, itxspnl,ngcq) 734 endif 735 if(iorder.eq.2143) then 736c ijkl->jilk 737 ls=ksize 738 is=lsize*ls 739 js=isize*is 740 ks=0 741 im=0 742 jm=0 743 km=1 744 lm=0 745 call txspnl(is,im, js,jm, ks,km ,ls,lm, itxspnl,ngcq) 746 endif 747c 748 if(iorder.eq.3412) then 749c ijkl->klij 750 is=jsize 751 ls=isize*is 752 ks=lsize*ls 753 js=0 754 im=0 755 jm=1 756 km=0 757 lm=0 758 call txspnl(is,im, js,jm, ks,km ,ls,lm, itxspnl,ngcq) 759 endif 760 if(iorder.eq.4312) then 761c ijkl->klji 762 is=0 763 js=isize 764 ls=jsize*js 765 ks=lsize*ls 766 im=1 767 jm=0 768 km=0 769 lm=0 770 call txspnl(is,im, js,jm, ks,km ,ls,lm, itxspnl,ngcq) 771 endif 772 if(iorder.eq.3421) then 773c ijkl->lkij 774 js=0 775 is=jsize 776 ks=isize*is 777 ls=ksize*ks 778 im=0 779 jm=1 780 km=0 781 lm=0 782 call txspnl(is,im, js,jm, ks,km ,ls,lm, itxspnl,ngcq) 783 endif 784 if(iorder.eq.4321) then 785c ijkl->lkji 786 is=0 787 js=isize 788 ks=jsize*js 789 ls=ksize*ks 790 im=1 791 jm=0 792 km=0 793 lm=0 794 call txspnl(is,im, js,jm, ks,km ,ls,lm, itxspnl,ngcq) 795 endif 796c 797 end 798c====================================================================== 799 subroutine txs_pnl(nquart,ncshell, 800 * icstxs,jcstxs,kcstxs,lcstxs, 801 * icspnl,jcspnl,kcspnl,lcspnl) 802c 803 dimension ncshell(*) 804 dimension 805 * icspnl(nquart),jcspnl(nquart),kcspnl(nquart),lcspnl(nquart), 806 * icstxs(nquart),jcstxs(nquart),kcstxs(nquart),lcstxs(nquart) 807c 808c txs-order pnl-order 809c 810 do 10 iq=1,nquart 811 icstxs(iq)=ncshell(icspnl(iq)) 812 jcstxs(iq)=ncshell(jcspnl(iq)) 813 kcstxs(iq)=ncshell(kcspnl(iq)) 814 lcstxs(iq)=ncshell(lcspnl(iq)) 815ctest98 816c itxs=icstxs(iq) 817c ipnl=icspnl(iq) 818c jtxs=jcstxs(iq) 819c jpnl=jcspnl(iq) 820c ktxs=kcstxs(iq) 821c kpnl=kcspnl(iq) 822c ltxs=lcstxs(iq) 823c lpnl=ncshell(lcspnl(iq)) 824c write(6,*)' PNL shells =',ipnl,jpnl,kpnl,lpnl 825c write(6,*)' TXS shells =',itxs,jtxs,ktxs,ltxs 826ctest98 827 10 continue 828c 829 do 20 iq=1,nquart 830 icspnl(iq)=icstxs(iq) 831 jcspnl(iq)=jcstxs(iq) 832 kcspnl(iq)=kcstxs(iq) 833 lcspnl(iq)=lcstxs(iq) 834 20 continue 835c 836 end 837c================================================================= 838 subroutine pnl_txs(ncfunct, icf,jcf,kcf,lcf,integ_n0 ) 839 dimension ncfunct(*),icf(*),jcf(*),kcf(*),lcf(*) 840c 841 do 10 integ=1,integ_n0 842 itxs=icf(integ) 843 jtxs=jcf(integ) 844 ktxs=kcf(integ) 845 ltxs=lcf(integ) 846c 847 ipnl=ncfunct(itxs) 848 jpnl=ncfunct(jtxs) 849 kpnl=ncfunct(ktxs) 850 lpnl=ncfunct(ltxs) 851c 852 icf(integ)=ipnl 853 jcf(integ)=jpnl 854 kcf(integ)=kpnl 855 lcf(integ)=lpnl 856 10 continue 857c 858 end 859c================================================================= 860 subroutine txspnl(is,im, js,jm, ks,km ,ls,lm,itxspnl,ngcq) 861 common /pnl004/ isize,jsize,ksize,lsize,iiiiiii 862 dimension itxspnl(*) 863c 864c Temprorarly re-order integrals txs-pnl 865c 866 increm=isize*jsize*ksize*lsize 867 itxs=0 868 do 10 iqu=1,ngcq 869 integ=(iqu-1)*increm 870 do 20 i=1,isize 871 ii=(i-1)*is +i*im 872 do 20 j=1,jsize 873 jj=(j-1)*js +j*jm 874 ij=ii+jj 875 do 20 k=1,ksize 876 kk=(k-1)*ks +k*km 877 ijk=ij+kk 878 do 20 l=1,lsize 879 ll=(l-1)*ls +l*lm 880 itxs=itxs+1 881ccccc ipnl=ijk+ll 882 ipnl=ijk+ll + integ 883 itxspnl(itxs)=ipnl 884 20 continue 885 10 continue 886c 887c write(6 ,*)'itxs , ipnl ( itxspnl(itxs))' 888c do 1111 ii=1,isize*jsize*ksize*lsize 889c write(6 ,66)ii,itxspnl(ii) 890c1111 continue 891c 66 format('ii=',i3,1x,'itxspnl(ii)=',i3) 892c 893 end 894c===================================================================== 895 subroutine switch_scr(bltxs, blscr,l_blscr) 896 implicit real*8 (a-h,o-z) 897c 898c copy from local bl() as defined in texas_face : 899c 900 common /memor1_R/ npard_R,mxsize_R,nblock1_R,nblock1_back_R 901 common /memor11_R/ mxpair_R 902c 903c to pnl : 904c 905 common /memor1/ npard,mxsize,nblock1,nblock1_back 906 common /memor11/ mxpair 907c 908c and for corresp. sizes : 909c 910 common /memor1_S/ npard_S,mxsize_S,nblock1_S,nblock1_back_S 911 common /memor11_S/ mxpair_S 912c 913 common /ganz/ lcore,iov,last,lflag(4),inuc,ibas,na,nbf,nsh,ncf,ncs 914 1,nsy(4),nXXX,nganz(35),lopt(30) 915c 916 dimension bltxs(*) 917 dimension blscr(l_blscr) 918c------------------------------- 919cccc write(6,*)' in switch_scr ; l_blscr=',l_blscr 920c------------------------------- 921 lcore=l_blscr 922c 923c now get/ret operate on blscr . 924c 925 call retall 926c------------------------------- 927c 928 call memo1_int(npard_S ,npard) 929 call memo1_int(mxsize_S,mxsize) 930 call memo1_int(mxpair_S,mxpair) 931 call memo1_int(nblock1_S,nblock1) 932 call memo1_int(nblock1_back_S,nblock1_back) 933c------------------------------- 934c copy data : 935c 936 call tfer_i(bltxs(npard_R), blscr(npard), npard_S ) 937 call tfer_i(bltxs(mxsize_R),blscr(mxsize), mxsize_S ) 938 call tfer_i(bltxs(mxpair_R),blscr(mxpair), mxpair_S ) 939 call tfer_i(bltxs(nblock1_R),blscr(nblock1), nblock1_S ) 940 call tfer_i 941 * (bltxs(nblock1_back_R),blscr(nblock1_back), nblock1_back_S ) 942c------------------------------- 943 end 944c========== 945 subroutine tfer_i(ia,ib,n) 946 dimension ia(n),ib(n) 947c 948 do 10 ii=1,n 949 ib(ii)=ia(ii) ! let the compiler do the unrolling. 950 10 continue 951c 952c$$$ m = mod(n,7) 953c$$$ if( m .eq. 0 ) go to 20 954c$$$ do 10 ii=1,n 955c$$$ ib(ii)=ia(ii) 956c$$$ 10 continue 957c$$$ if( n .lt. 7) return 958c$$$ 20 mp1=m+1 959c$$$ do 30 i=mp1,n,7 960c$$$ ib(i ) = ia(i ) 961c$$$ ib(i+1) = ia(i+1) 962c$$$ ib(i+2) = ia(i+2) 963c$$$ ib(i+3) = ia(i+3) 964c$$$ ib(i+4) = ia(i+4) 965c$$$ ib(i+5) = ia(i+5) 966c$$$ ib(i+6) = ia(i+6) 967c$$$ 30 continue 968c 969 end 970c===================================================================== 971 subroutine supblks_split(bl,nquart,maxsize,maxpair, 972 * ijblock,klblock, 973 * leri,inx, nblock1,nbl2, 974 * isbl_q,isbl_point, 975 * ijpres2,klpres2,isbl_s,isbl_part) 976c-------------------------------------------------------------------- 977c this routine is called for every new PNL request i.e when 978c more_int=.false. 979c-------------------------------------------------------------------- 980 implicit real*8 (a-h,o-z) 981#include "errquit.fh" 982cccc common /intlim/ limxmem,limblks,limpair 983c 984 common /pnl_nqrt/ ncall_pnl,ncall_41,nqrts_pnl,npart_pnl,nsize_pnl 985c 986c output from uniq_pairs_1 call: 987c 988 common /map4uniq/ nij_uniqe, ij_uniqe_p, map_ij, 989 * nkl_uniqe, kl_uniqe_p, map_kl 990c 991 dimension maxsize(*),maxpair(*) 992 dimension inx(12,*),nblock1(*) 993 dimension isbl_s(*) ! dimension nbl2*(nbl2+1)/2 num of super-blks 994 dimension ijpres2(*),klpres2(*) ! dimension nquart 995 dimension isbl_part(*) ! dimension nbl2*(nbl2+1)/2 996 dimension isbl_q(*), isbl_point(*) 997c 998 dimension ijblock(nbl2),klblock(nbl2) 999 dimension bl(*) 1000c--------------------------------------------------------------------- 1001c for derivatives of all kinds : 1002cccc common /memmax/ ispblx, maxme1,iforwhat 1003 common /mem_max_min/ ispblx,maxme1,max_111,iforwhat 1004c--------------------------------------------------------------------- 1005 n_times=1 1006 if(iforwhat.eq.2) n_times= 6 ! giao Ist derivatives 1007 if(iforwhat.eq.3) n_times=12 ! gradient derivatives 1008 if(iforwhat.eq.4) n_times=78 ! hessian derivatives 1009c--------------------------------------------------------------------- 1010ctest 1011c write(6,*)'call=',ncall_pnl,' leri=',leri,' nquart=',nquart 1012c--------------------------------------------------------------------- 1013 i_blks=0 1014 isupb=0 1015 do 10 ibl=1,nbl2 1016 if(ijblock(ibl).eq.0) then 1017 iibl=ibl*(ibl-1)/2 1018 iibl1=iibl+1 1019 iibl2=iibl+ibl 1020 do iiii=iibl1,iibl2 1021 isbl_part(iiii)=1 1022 enddo 1023 isupb=isupb+ibl 1024 go to 10 1025 endif 1026 call get_ics_jcs(nblock1,ibl,1,ics1,jcs1) 1027 icont_ij=(inx(4,ics1)+1)*(inx(4,jcs1)+1) 1028 integ_ij=inx(3,ics1)*inx(3,jcs1)*icont_ij 1029 do 20 kbl=1,ibl 1030 isupb=isupb+1 1031 isbl_part(isupb)=1 1032 isbl_size=isbl_s(isupb) 1033 if(isbl_size.gt.0) then 1034c 1035 call get_ics_jcs(nblock1,kbl,1,kcs1,lcs1) 1036 icont_kl=(inx(4,kcs1)+1)*(inx(4,lcs1)+1) 1037 integ_kl=inx(3,kcs1)*inx(3,lcs1)*icont_kl 1038 integ1=integ_ij*integ_kl 1039c 1040ccccccc integ1=integ1*n_times ! derivatives 1041c 1042c first re-define maximum size according to leri : 1043c 1044 maxm_size=maxsize(isupb) 1045 leri_size=leri/integ1 1046c 1047 if(leri_size.eq.0) then 1048c 1049 write(6,*)' buffer size too small; leri=',leri, 1050 * ' needed=',integ1,' =',n_times,'*',integ1/n_times 1051 write(6,*)' shell sizes=', 1052 * inx(3,ics1),inx(3,jcs1),inx(3,kcs1),inx(3,lcs1), 1053 * ', shell g.con=',inx(4,ics1)+1,inx(4,jcs1)+1, 1054 * inx(4,kcs1)+1,inx(4,lcs1)+1 1055 call errquit('texas: supblks_split',1016, INT_ERR) 1056 endif 1057c 1058 maxsize(isupb)=min(leri_size,maxm_size) 1059c 1060c split each super-block in parts if its current pnl-size is greater 1061c than maximum-size allowed or number of pairs exceeds allowed limit 1062c 1063 maxm_size=maxsize(isupb) 1064 maxm_pair=maxpair(isupb) 1065c 1066 ipoint=isbl_point(isupb) 1067c 1068 call uniq_pairs_1('ks_split',ibl,kbl,bl,ijpres2,klpres2, 1069 * isbl_size,isbl_q,ipoint) 1070 call retmem(2) ! release allocations in uniq_pairs_1 1071c output nij_uniqe & nkl_uniqe in common /map4uniq/ 1072c 1073 if(nij_uniqe.le.maxm_pair.and.nkl_uniqe.le.maxm_pair) then 1074 if(isbl_size.gt.maxm_size) then 1075 nparts=isbl_size/maxm_size 1076 nrem =mod(isbl_size,maxm_size) 1077 if(nrem.gt.0) nparts=nparts+1 1078 isbl_part(isupb)=nparts 1079 maxsize(isupb)=maxm_size 1080 endif 1081 else 1082 maxm_size=min(maxm_size,maxm_pair) 1083 nparts=isbl_size/maxm_size 1084 nrem =mod(isbl_size,maxm_size) 1085 if(nrem.gt.0) nparts=nparts+1 1086 isbl_part(isupb)=nparts 1087 maxsize(isupb)=maxm_size 1088 endif 1089c 1090 i_blks=i_blks+isbl_part(isupb) 1091c 1092 endif 1093 20 continue 1094 10 continue 1095c--------------------------------------------------------- 1096 if(i_blks.gt.0) nsize_pnl=nsize_pnl+ nquart/i_blks 1097c--------------------------------------------------------- 1098 end 1099c===================================================================== 1100 subroutine request_split(leri,inx,nbl2, ijblock,klblock, 1101 * ijpres2,klpres2, 1102 * isbl_s,isbl_q, isbl_part, 1103 * isbl_split,isplit) 1104c-------------------------------------------------------------------- 1105c if a given pnl-request is too big then it is divided into parts: 1106c-------------------------------------------------------------------- 1107 dimension inx(12,*) 1108 dimension isbl_s(*) ! dimension nbl2*(nbl2+1)/2 1109 dimension isbl_q(*) ! dimension nquart 1110 dimension ijpres2(*),klpres2(*) ! dimension nquart 1111 dimension ijblock(nbl2),klblock(nbl2) 1112 dimension isbl_part(*) ! dim. nbl2*(nbl2+1)/2 1113 dimension isbl_split(0:*) 1114c--------------------------------------------------------------------- 1115c for derivatives of all kinds : 1116ccc common /memmax/ ispblx, maxme1,iforwhat 1117 common /mem_max_min/ ispblx,maxme1,max_111,iforwhat 1118c--------------------------------------------------------------------- 1119 n_times=1 1120 if(iforwhat.eq.2) n_times= 6 ! giao Ist derivatives 1121 if(iforwhat.eq.3) n_times=12 ! gradient derivatives 1122 if(iforwhat.eq.4) n_times=78 ! hessian derivatives 1123c--------------------------------------------------------------------- 1124c 1125 isplit=0 1126 integrals=0 1127 ipoint=0 1128c 1129 do 50 ibl=1,nbl2 1130 if(ijblock(ibl).eq.0) go to 50 1131 iibl=ibl*(ibl-1)/2 1132 do 60 kbl=1,ibl 1133 isupb=iibl+kbl 1134 isbl_size=isbl_s(isupb) 1135 nparts=isbl_part(isupb) 1136 if(isbl_size.eq.0) go to 60 1137ctest 1138ccc write(6,*)'block=',isupb,' size=',isbl_size,' nparts=',nparts 1139ctest 1140c 1141 iq1=isbl_q(ipoint+1) 1142 ijcs1=ijpres2(iq1) 1143 klcs1=klpres2(iq1) 1144 call get_ij_half(ijcs1,ics1,jcs1) 1145 icont_ij=(inx(4,ics1)+1)*(inx(4,jcs1)+1) 1146 call get_ij_half(klcs1,kcs1,lcs1) 1147 icont_kl=(inx(4,kcs1)+1)*(inx(4,lcs1)+1) 1148 integ1=inx(3,ics1)*inx(3,jcs1)*inx(3,kcs1)*inx(3,lcs1) 1149 integ1=integ1*icont_ij*icont_kl 1150c 1151cccccc integ1=integ1*n_times ! derivatives 1152c 1153c splitting among super-blocks : 1154c 1155 if(nparts.eq.1) then 1156 integ=integ1*isbl_size 1157 if((integrals+integ).gt.leri) then 1158 isplit=isplit+1 1159 integrals=integ 1160 isbl_split(isplit)=isupb 1161 else 1162 integrals=integrals+integ 1163 endif 1164 else 1165c splitting inside of one super-block : force splitting : 1166 integrals=leri 1167 isplit=isplit+1 1168 isbl_split(isplit)=isupb 1169 endif 1170c 1171 ipoint=ipoint+isbl_size 1172 60 continue 1173 50 continue 1174c 1175 isblast=nbl2*(nbl2+1)/2 1176 isbl_split(0)=1 1177 isbl_split(isplit+1)=isblast+1 1178c--------------------------------------------------------- 1179 end 1180c===================================================================== 1181 subroutine supblks_point(nsupblk,isbl_size,isbl_point) 1182 dimension isbl_size(*), isbl_point(*) 1183c 1184c calculates the pointer to block's quartets : iq=isbl_q(IPOINT+iqp) 1185c 1186 ipoint=0 1187 isbl_point(1)=0 1188 do 10 isbl=2,nsupblk 1189 ipoint=ipoint+isbl_size(isbl-1) 1190 isbl_point(isbl)=ipoint 1191 10 continue 1192c 1193 end 1194c===================================================================== 1195 subroutine memstat_pnl 1196 common /mem_pnl_scr/ nall_peak,mark_peak,mem_peak 1197c 1198 call memstat(nalloc_bl,nmark_bl,maxmem_bl,memtot_bl) 1199c 1200 if(nalloc_bl.gt.nall_peak) nall_peak=nalloc_bl 1201 if(nmark_bl .gt.mark_peak) mark_peak=nmark_bl 1202 if(maxmem_bl.gt. mem_peak) mem_peak=maxmem_bl 1203c 1204 end 1205c===================================================================== 1206 subroutine request_update(icspnl,jcspnl,kcspnl,lcspnl,nquart, 1207 * num_bas_ij,num_bas_kl) 1208c 1209c--------------------------------------------------------------------- 1210c I,J shells might belong to different basis set than K,L shells 1211c--------------------------------------------------------------------- 1212c num_bas_ij , num_bas_kl = basis sets for IJ,KL . 1213c For 4-center (ordinary) two-electron integrals ALL four vectors are 1214c non-zero (just numbers of contracted shells). 1215c For 3-center two-el. integrals ONE of these vectors should be ZERO. 1216c For 2-center two-el. integrals TWO of these vectors should be ZERO. 1217c 1218c If ZERO is found then it is replaced by NCS-the LAST contracted shell 1219c (in PNL order). This is the S-type shell (uncontracted & exp.=zero) 1220c which was added to the basis set. 1221c--------------------------------------------------------------------- 1222c 1223 logical cent2,cent3,cent4 1224 common /what_was_calc/ cent2,cent3,cent4 1225c--------------------------------------------------------------------- 1226 common /ganz/ lcore,iov,last,lflag(4),inuc,ibas,na,nbf,nsh,ncf,ncs 1227 1,nsy(4),nsym,nganz(35),lopt(30) 1228c--------------------------------------------------------------------- 1229 common /multi_basis/ num_bas_1,num_bas_2,num_bas_3, 1230 * ncs_bas_1,ncs_bas_2,ncs_bas_3, 1231 * nps_bas_1,nps_bas_2,nps_bas_3, 1232 * nat_bas_1,nat_bas_2,nat_bas_3, 1233 * ncf_bas_1,ncf_bas_2,ncf_bas_3 1234c 1235 dimension icspnl(nquart),jcspnl(nquart), 1236 * kcspnl(nquart),lcspnl(nquart) 1237c--------------------------------------------------------------------- 1238ctest 1239c write(6,*)'ncs_bas_1-3=',ncs_bas_1,ncs_bas_2,ncs_bas_3 1240c write(6,*)' entr iscpnl :', icspnl 1241c write(6,*)' entr jscpnl :', jcspnl 1242c write(6,*)' entr kscpnl :', kcspnl 1243c write(6,*)' entr lscpnl :', lcspnl 1244ctest 1245c 1246 ics=icspnl(1) 1247 jcs=jcspnl(1) 1248 kcs=kcspnl(1) 1249 lcs=lcspnl(1) 1250c 1251 ncenters=0 1252 if(ics.gt.0) ncenters=ncenters+1 1253 if(jcs.gt.0) ncenters=ncenters+1 1254 if(kcs.gt.0) ncenters=ncenters+1 1255 if(lcs.gt.0) ncenters=ncenters+1 1256c 1257c write(6,*) ncenters,'-center integrals are requested' 1258c 1259 if(ncenters.eq.4) then 1260 cent4=.true. 1261c 1262c all possible cases 1263c ----------------- 1264c if(num_bas_ij.eq.num_bas_1) then 1265c endif 1266 if(num_bas_ij.eq.num_bas_2) then 1267 call ncshl_update(icspnl,nquart,ncs_bas_1) 1268 call ncshl_update(jcspnl,nquart,ncs_bas_1) 1269 endif 1270 if(num_bas_ij.eq.num_bas_3) then 1271 call ncshl_update(icspnl,nquart,ncs_bas_2) 1272 call ncshl_update(jcspnl,nquart,ncs_bas_2) 1273 endif 1274c 1275c if(num_bas_kl.eq.num_bas_1) then 1276c endif 1277 if(num_bas_kl.eq.num_bas_2) then 1278 call ncshl_update(kcspnl,nquart,ncs_bas_1) 1279 call ncshl_update(lcspnl,nquart,ncs_bas_1) 1280 endif 1281 if(num_bas_kl.eq.num_bas_3) then 1282 call ncshl_update(kcspnl,nquart,ncs_bas_2) 1283 call ncshl_update(lcspnl,nquart,ncs_bas_2) 1284 endif 1285 endif 1286c 1287 if(ncenters.eq.3) then 1288 cent3=.true. 1289c 1290c all possible cases 1291c ----------------- 1292 if(ics.eq.0) then 1293 do ii=1,nquart 1294 icspnl(ii)=ncs 1295 enddo 1296 if(num_bas_ij.eq.num_bas_2) then 1297 call ncshl_update(jcspnl,nquart,ncs_bas_1) 1298 endif 1299 if(num_bas_ij.eq.num_bas_3) then 1300 call ncshl_update(jcspnl,nquart,ncs_bas_2) 1301 endif 1302 if(num_bas_kl.eq.num_bas_2) then 1303 call ncshl_update(kcspnl,nquart,ncs_bas_1) 1304 call ncshl_update(lcspnl,nquart,ncs_bas_1) 1305 endif 1306 if(num_bas_kl.eq.num_bas_3) then 1307 call ncshl_update(kcspnl,nquart,ncs_bas_2) 1308 call ncshl_update(lcspnl,nquart,ncs_bas_2) 1309 endif 1310 endif 1311c ----------------- 1312 if(jcs.eq.0) then 1313 do ii=1,nquart 1314 jcspnl(ii)=ncs 1315 enddo 1316 if(num_bas_ij.eq.num_bas_2) then 1317 call ncshl_update(icspnl,nquart,ncs_bas_1) 1318 endif 1319 if(num_bas_ij.eq.num_bas_3) then 1320 call ncshl_update(icspnl,nquart,ncs_bas_2) 1321 endif 1322 if(num_bas_kl.eq.num_bas_2) then 1323 call ncshl_update(kcspnl,nquart,ncs_bas_1) 1324 call ncshl_update(lcspnl,nquart,ncs_bas_1) 1325 endif 1326 if(num_bas_kl.eq.num_bas_3) then 1327 call ncshl_update(kcspnl,nquart,ncs_bas_2) 1328 call ncshl_update(lcspnl,nquart,ncs_bas_2) 1329 endif 1330 endif 1331c ----------------- 1332 if(kcs.eq.0) then 1333 do ii=1,nquart 1334 kcspnl(ii)=ncs 1335 enddo 1336 if(num_bas_ij.eq.num_bas_2) then 1337 call ncshl_update(icspnl,nquart,ncs_bas_1) 1338 call ncshl_update(jcspnl,nquart,ncs_bas_1) 1339 endif 1340 if(num_bas_ij.eq.num_bas_3) then 1341 call ncshl_update(icspnl,nquart,ncs_bas_2) 1342 call ncshl_update(jcspnl,nquart,ncs_bas_2) 1343 endif 1344 if(num_bas_kl.eq.num_bas_2) then 1345 call ncshl_update(lcspnl,nquart,ncs_bas_1) 1346 endif 1347 if(num_bas_kl.eq.num_bas_3) then 1348 call ncshl_update(lcspnl,nquart,ncs_bas_2) 1349 endif 1350 endif 1351c ----------------- 1352 if(lcs.eq.0) then 1353 do ii=1,nquart 1354 lcspnl(ii)=ncs 1355 enddo 1356 if(num_bas_ij.eq.num_bas_2) then 1357 call ncshl_update(icspnl,nquart,ncs_bas_1) 1358 call ncshl_update(jcspnl,nquart,ncs_bas_1) 1359 endif 1360 if(num_bas_ij.eq.num_bas_3) then 1361 call ncshl_update(icspnl,nquart,ncs_bas_2) 1362 call ncshl_update(jcspnl,nquart,ncs_bas_2) 1363 endif 1364 if(num_bas_kl.eq.num_bas_2) then 1365 call ncshl_update(kcspnl,nquart,ncs_bas_1) 1366 endif 1367 if(num_bas_kl.eq.num_bas_3) then 1368 call ncshl_update(kcspnl,nquart,ncs_bas_2) 1369 endif 1370 endif 1371 endif 1372c 1373 if(ncenters.eq.2) then 1374 cent2=.true. 1375c 1376c only j=0 & l=0 case !!! 1377c ----------------- 1378 do ii=1,nquart 1379 jcspnl(ii)=ncs 1380 lcspnl(ii)=ncs 1381 enddo 1382 if(num_bas_ij.eq.num_bas_2) then 1383 call ncshl_update(icspnl,nquart,ncs_bas_1) 1384 endif 1385 if(num_bas_ij.eq.num_bas_3) then 1386 call ncshl_update(icspnl,nquart,ncs_bas_2) 1387 endif 1388 if(num_bas_kl.eq.num_bas_2) then 1389 call ncshl_update(kcspnl,nquart,ncs_bas_1) 1390 endif 1391 if(num_bas_kl.eq.num_bas_3) then 1392 call ncshl_update(kcspnl,nquart,ncs_bas_2) 1393 endif 1394ctest 1395c write(6,*)' exit iscpnl :', icspnl 1396c write(6,*)' exit jscpnl :', jcspnl 1397c write(6,*)' exit kscpnl :', kcspnl 1398c write(6,*)' exit lscpnl :', lcspnl 1399ctest 1400 endif 1401c--------------------------------------------------------------------- 1402 end 1403c===================================================================== 1404 subroutine ncshl_update(icspnl,nquart,ncs_bas_1) 1405 dimension icspnl(nquart) 1406c 1407 do ii=1,nquart 1408 icspnl(ii)=icspnl(ii)+ncs_bas_1 1409 enddo 1410c 1411 end 1412c===================================================================== 1413 subroutine labels_update(icf,jcf,kcf,lcf,integ,ij_basis,kl_basis, 1414 * ish_first,jsh_first,ksh_first,lsh_first) 1415c--------------------------------------------------------------------- 1416 common /multi_basis/ num_bas_1,num_bas_2,num_bas_3, 1417 * ncs_bas_1,ncs_bas_2,ncs_bas_3, 1418 * nps_bas_1,nps_bas_2,nps_bas_3, 1419 * nat_bas_1,nat_bas_2,nat_bas_3, 1420 * ncf_bas_1,ncf_bas_2,ncf_bas_3 1421 dimension icf(*),jcf(*),kcf(*),lcf(*) 1422c--------------------------------------------------------------------- 1423ctest 1424c write(6,*)'ncf_bas_1-3=',ncf_bas_1,ncf_bas_2,ncf_bas_3 1425c write(6,*)' entr icf1-9 :',(icf(ii),ii=1,9) 1426c write(6,*)' entr jcf1-9 :',(jcf(ii),ii=1,9) 1427c write(6,*)' entr kcf1-9 :',(kcf(ii),ii=1,9) 1428c write(6,*)' entr lcf1-9 :',(lcf(ii),ii=1,9) 1429ctest 1430c 1431 if(ij_basis.eq.num_bas_1 .and. kl_basis.eq.num_bas_1) RETURN 1432c 1433c ij 1434c 1435 if(ij_basis.eq.num_bas_2) then 1436 if(ish_first.gt.0) call ncfun_update(icf,integ,ncf_bas_1) 1437 if(jsh_first.gt.0) call ncfun_update(jcf,integ,ncf_bas_1) 1438 endif 1439c 1440 if(ij_basis.eq.num_bas_3) then 1441 if(ish_first.gt.0) call ncfun_update(icf,integ,ncf_bas_2) 1442 if(jsh_first.gt.0) call ncfun_update(jcf,integ,ncf_bas_2) 1443 endif 1444c 1445c kl 1446c 1447 if(kl_basis.eq.num_bas_2) then 1448 if(ksh_first.gt.0) call ncfun_update(kcf,integ,ncf_bas_1) 1449 if(lsh_first.gt.0) call ncfun_update(lcf,integ,ncf_bas_1) 1450 endif 1451c 1452 if(kl_basis.eq.num_bas_3) then 1453 if(ksh_first.gt.0) call ncfun_update(kcf,integ,ncf_bas_2) 1454 if(lsh_first.gt.0) call ncfun_update(lcf,integ,ncf_bas_2) 1455 endif 1456ctest 1457c write(6,*)' exit icf1-9 :',(icf(ii),ii=1,9) 1458c write(6,*)' exit jcf1-9 :',(jcf(ii),ii=1,9) 1459c write(6,*)' exit kcf1-9 :',(kcf(ii),ii=1,9) 1460c write(6,*)' exit lcf1-9 :',(lcf(ii),ii=1,9) 1461ctest 1462c 1463 end 1464c===================================================================== 1465 subroutine ncfun_update(icf,integ,ncf_bas_1) 1466 dimension icf(*) 1467c 1468 do ii=1,integ 1469 icf(ii)=icf(ii)-ncf_bas_1 1470 enddo 1471c 1472 end 1473c===================================================================== 1474 subroutine requested_task(int_type) 1475 character*8 int_type 1476 character*27 type_i,type_r 1477 character*11 scftype 1478 character*8 where 1479 common /runtype/ scftype,where 1480 common /type_inited/ iforinit ! use here & in texas_face.F 1481 common /mem_max_min/ ispblx,maxme1,max_111,iforwhat 1482c 1483 irequest = 1 ! take care of compiler warnings 1484 if(int_type.eq.'scfd_int') then 1485 where='buff' 1486 irequest=1 1487 endif 1488 if(int_type.eq.'giao_int') then 1489 where='shif' 1490 irequest=2 1491 endif 1492 if(int_type.eq.'der1_int') then 1493 where='forc' 1494 irequest=3 1495 endif 1496 if(int_type.eq.'der2_int') then 1497 where='hess' 1498 irequest=4 1499 endif 1500c 1501c check if the texas integral program has been initiated 1502c for up to this requested task : 1503c 1504 if(irequest.gt.iforinit) then 1505 type_i = ' texas not initialized ' 1506 if(iforinit.eq.1) type_i='ordinary 2-el.integrals ' 1507 if(iforinit.eq.2) type_i='giao integral derivatives ' 1508 if(iforinit.eq.3) type_i='first geometry derivatives' 1509 if(iforinit.eq.4) type_i='second geometry derivatives' 1510c 1511 type_r = ' texas error on irequest' 1512 if(irequest.eq.1) type_r='ordinary 2-el.integrals ' 1513 if(irequest.eq.2) type_r='giao integral derivatives ' 1514 if(irequest.eq.3) type_r='first geometry derivatives' 1515 if(irequest.eq.4) type_r='second geometry derivatives' 1516c 1517 write(6,66) type_i,type_r 1518 66 format(' Texas integral program '/ 1519 & ' has been initialized only for :',a27,/, 1520 & ' and can not be executed for :',a27) 1521 call errquit('texas program stopped in requested_task',0, 1522 & INT_ERR) 1523 else 1524 iforwhat=irequest 1525 endif 1526c 1527 end 1528c================================================================== 1529 subroutine get_ics_jcs(nblock1,ijblock,ijpar,ics,jcs) 1530c-------------------------------------------------- 1531c ijblock - pair-block ![input] 1532c ijpar - a pair in that block ![input] 1533c ics,jcs : shells ![output] 1534c-------------------------------------------------- 1535 dimension nblock1(0:* ) ! nblock1(0:nbl1) 1536c 1537 call get_ij_half(ijblock, iblock,jblock) 1538c 1539 ics_b= nblock1(iblock-1)+1 ! first shell 1540 ics_e= nblock1(iblock) ! last shell 1541 jcs_b= nblock1(jblock-1)+1 ! first shell 1542 jcs_e= nblock1(jblock) ! last shell 1543c 1544c increm=ics_e-ics_b+1 1545c jncrem=jcs_e-jcs_b+1 1546c 1547 if(jblock.eq.iblock) then 1548 call get_ij_half(ijpar,ishell,jshell) 1549 else 1550 jncrem=jcs_e-jcs_b+1 1551 call get_ij_full(ijpar,jncrem,ishell,jshell) 1552 endif 1553c 1554c ishell,jshell - shells in the ijblock (numbered from 1 to increm 1555c 1 to jncrem) 1556c 1557 ics=ishell+ics_b -1 1558 jcs=jshell+jcs_b -1 1559c 1560 end 1561c================================================================== 1562 subroutine blk_capab(nbl2,npar,isbl_s,nquart_pnl,nquart_txs) 1563 implicit real*8 (a-h,o-z) 1564 common /pnl000/ xbluse,nbluse 1565 dimension isbl_s(*) , npar(*) 1566c 1567 xquart_pnl=dble(nquart_pnl) 1568 xquart_txs=dble(nquart_txs) 1569c 1570 ibluse=0 1571 ikbl=0 1572 do ibl=1,nbl2 1573 do kbl=1,ibl 1574 ikbl=ikbl+1 1575 iquart_pnl=isbl_s(ikbl) 1576 if(iquart_pnl.gt.0) ibluse=ibluse+1 1577 enddo 1578 enddo 1579c 1580 nbluse=nbluse+ibluse 1581c 1582 xblk_pnl=dble(ibluse) 1583 xblk_txs=dble(ikbl) 1584c 1585 size_pnl=xquart_pnl/xblk_pnl 1586 size_txs=xquart_txs/xblk_txs 1587 bluse=size_pnl/size_txs 1588 1589 xbluse=xbluse+bluse 1590c 1591 end 1592c================================================================== 1593