1* $Id$ 2c=================================================================== 3 subroutine calcint2(bl,inx, q4,use_q4, 4 * lab1,lab2,lab3,lab4,fockmat, 5 $ ncfunct) 6 implicit real*8 (a-h,o-z) 7 integer ncfunct(*) ! Map from txs to pnl basis functions 8 logical firstd 9 logical use_q4 10 character*11 scftype 11 character*8 where 12c---------------------------------------------------------------- 13c The variable where shows the type of calculated integrals: 14c where.eq.'buff' - for ordinary two-el. integrals 15c where.eq.'shif' - for the GIAO/NMR integral derivativs 16c where.eq.'forc' - for first geometry derivatives 17c where.eq.'hess' - for second geometry derivatives 18c It can be extanded for other targets. 19c---------------------------------------------------------------- 20 common /runtype/ scftype,where 21 common /txs_time/ timprep2,timcalc2,timblok2 22 common /number/ acc 23 common /infor/ icheck,firstd,ndirect,nprint,iblok,nbeg,nend 24 common /memor1/ npard,mxsize,nblock1,nblock1_back 25 common /memor1b/ nbl2 26c 27 dimension bl(*) 28 dimension inx(12,*) 29 dimension fockmat(*) 30 dimension lab1(*),lab2(*),lab3(*),lab4(*) 31 dimension q4(*) ! pnl-symmetry factors 32c---------------------------------------------------------------- 33 call txs_second(time_beg) 34c---------------------------------------------------------------- 35c calculations of two-electron integrals in blocks 36c Nbl2 - number of blocks of pairs 37c 38 call blockint(bl,nbl2,inx,bl(nblock1),bl(npard), 39 * bl(mxsize), 40 * q4,use_q4, 41 * lab1,lab2,lab3,lab4,fockmat,ncfunct) 42c 43c---------------------------------------------------------------- 44 call txs_second(time_end) 45 timcalc2=timcalc2 + time_end-time_beg 46c---------------------------------------------------------------- 47 end 48c=================================================================== 49 subroutine blockint(bl,nbl2,inx,nblock1,npar, 50 * mxsize, 51 * q4,use_q4, 52 * lab1,lab2,lab3,lab4,fockmat, ncfunct) 53c 54 implicit real*8 (a-h,o-z) 55#include "errquit.fh" 56 integer ncfunct(*) ! txs to pnl basis map 57 logical firstd 58 logical use_q4 59 character*11 scftype 60 character*8 where 61c---------------------------------------------------------------- 62c where is 'disk' or 'fock' for SCF : non -,semi- and full-direct 63c It can be extanded for special purposes : 64c it is 'shif' or 'sf10' for nmr shift calculations 65c it is 'forc' for gradient 66c---------------------------------------------------------------- 67c ispec - indicates a special use of the two-el.int.program (pnl) 68 common /pnl000/ xbluse,nbluse 69 common /pnl001/ ispec,ijpres2,klpres2, ijblock,klblock,iqorder 70 common /pnl005/ isblsize,isblqrts,isblpoint 71 common /pnl006/ nsplit,isplit,isbl_split, isbl_part 72 common /map4uniq/ nij_uniqe, ij_uniqe_p, map_ij, 73 * nkl_uniqe, kl_uniqe_p, map_kl 74c------ 75 common /superbl/ isupb,ibl,kbl 76 common /runtype/ scftype,where 77 common /infor/ icheck,firstd,ndirect,nprint,iblok,nbeg,nend 78 common /howmany/ ntotal,noijkl,nopres,nospec,nrimtot,nrimret 79 dimension bl(*), inx(12,*) 80cnew dimension ncost(*),mxsize(*) 81 dimension mxsize(*) 82 dimension nblock1(*),npar(*) 83 dimension fockmat(*) 84 dimension lab1(*),lab2(*),lab3(*),lab4(*) 85 dimension q4(*) 86c---------------------------------------------------------------- 87c setup info needed if SPEC (spec.use of two-el.prog) is ON 88c 89 if(ispec.ne.0) scftype='** PNL **' 90c---------------------------------------------------------------- 91c for checking only : 92 call getmem(1,last0) 93c write(8,*)'From Blockint (entry): first free address in BL=',last0 94 call retmem(1) 95c------------------------------------------------------------------c 96c begining of the loop over blocks of contracted quartets c 97c------------------------------------------------------------------c 98c loop over quartet-blocks : nbl4 =nbl2*(nbl2+1)/2 99c 100 isupb=0 101 do 100 ibl=1,nbl2 102 call doblock2(ibl,bl(ijblock),idoit) 103 if(idoit.eq.0) go to 100 104 ibl12=ibl*(ibl-1)/2 105 do 200 kbl=1,ibl 106 call doblock2(kbl,bl(klblock),jdoit) 107 if(jdoit.eq.0) go to 200 108 isupb=ibl12+kbl 109c 110 call get_nbls_pnl(isupb,bl(isblsize),bl(isblpoint), 111 * isbl_size,isbl_point) 112c 113 call dosupblk(isupb,bl(isblsize),bl(isblqrts),isbl_point, 114 * isbdo,bl(isbl_split) ) 115 if(isbdo.eq.0) go to 200 116c 117 if (where.eq.'buff' .or. 118 * where.eq.'shif' .or. 119 * where.eq.'forc' .or. 120 * where.eq.'hess' ) then 121c 122 maxsize=mxsize(isupb) 123c 124c for each super-block do update to the current pnl-sizes : 125c 126 call get_nparts(isupb,bl(isbl_part),nparts) 127c 128c check in quartets for a super-block which is too big : 129c 130 if(nparts.gt.1) then 131c-----------> call getmem(isbl_size,isbl_copy) ! <-----------------c 132 call memo1_int(isbl_size,isbl_copy) 133 call copy_block(isbl_size,bl(isblqrts),isbl_point, 134 * bl(isbl_copy) ) 135 call check_in(isupb,isbl_size,bl(isblqrts),isbl_point, 136 * bl(isbl_copy),maxsize) 137 else 138 isbl_copy = 0 139 endif 140c 141 call uniq_pairs_1('blockint', 142 * ibl,kbl,bl,bl(ijpres2),bl(klpres2), 143 * isbl_size,bl(isblqrts),isbl_point) 144ctest 145c WRITE(6,*) 146c * 'SUPER=',isupb,' IBL,KBL=',ibl,kbl,' PNL-SIZE=',isbl_size, 147c * 'MAX-SIZE=',maxsize,' PARTS=',nparts 148ctest 149c ! two calls to getmem 150 call prec2ij(ibl,nblock1,bl,inx) ! one call to getmem 151 call prec2kl(kbl,nblock1,bl,inx) ! one call to getmem 152c 153c 154c make map from present quartets to unique pairs in a given block 155c 156 call make_map2uniq(bl,ibl,kbl, 157 * isbl_size,bl(isblqrts),isbl_point, 158 * bl(ijpres2),bl(klpres2) ) 159c 160c ! two calls to getmem 161c 162c in above call there are two memory allocations (getmem=memo1_int) 163c output from this call is in the map4uniq common block 164c 165c Both, isbl_size and maxsize are needed in the following call 166c for the case when a block is split into parts. Memory in memo5c should 167c be reserved according to min(isbl_size,maxsize) while in memo3 ALWAYS 168c according to isbl_size : 169c 170 call onesuper(isupb,ibl,kbl, isbl_size,maxsize, 171 * bl,inx,npar,nbl2, where, q4,use_q4, 172 * lab1,lab2,lab3,lab4 ,fockmat, ncfunct) 173c 174 call retmem(2) ! release mem.alloc. in make_map2uniq 175c 176c check out calculated quartests : 177c 178 call check_out(isupb,bl(isblsize),bl(isblqrts),isbl_point, 179 * bl(isbl_split),bl(isbl_copy),nparts,maxsize, 180 * ibl,bl(ijblock),igo_back ) 181c 182 call retmem(1) ! from prec2kl 183 call retmem(1) ! from prec2ij 184 call retmem(2) ! from uniq_pairs_1 185c 186 if(nparts.gt.1) call retmem(1) ! release isbl_copy 187c 188 if(igo_back.eq.1) go to 300 189c 190 endif 191 200 continue 192 100 continue 193c 194 300 continue 195c 196c------------------------------------------------------------------c 197c******* end of the loop over blocks of contracted quartets *******c 198c------------------------------------------------------------------c 199c 200c for checking only : 201 call getmem(1,last1) 202 call retmem(1) 203 if(last1.ne.last0) then 204 write(8,*)'** WRONG memory alocations **' 205 write(8,*)' at the beginning of calc.=',last0 206 write(8,*)' at the end of cal. =',last1 207 write(8,*)'** Execution has been stoped in BLOCKINT **' 208* stop100 209 call errquit('texas:blockint',0, INT_ERR) 210 endif 211c* 212 end 213c=================================================================== 214 subroutine onesuper(isupb,ibl,kbl, nbls_pnl,maxsize, 215 * bl,inx,npar,nbl2, where, q4,use_q4, 216 * lab1,lab2,lab3,lab4 ,fockmat,ncfunct) 217 implicit real*8 (a-h,o-z) 218#include "errquit.fh" 219 integer ncfunct(*) ! txs to pnl basis map 220 logical use_q4 221 character*8 where 222c---------------------------------------------------------------- 223c nbls_size = full pnl-size 224c maxsize = maximum size (can be smaller than above) 225c---------------------------------------------------------------- 226 common /memor2/ nqrtd, nibld,nkbld, nijbd,nijed, nklbd,nkled 227 common /memor3/ nblok1d 228 dimension bl(*), inx(12,*) 229 dimension npar(nbl2) 230 dimension fockmat(*) 231 dimension lab1(*),lab2(*),lab3(*),lab4(*) 232 dimension q4(*) ! pnl-symmetry factors 233c---------------------------------------------------------------- 234 call mmark 235c---------------------------------------------------------------- 236c Construct blocks of contracted shell quartets for a given 237c super-block (isupb) ; set up common /memor2/ and /memor3/ . 238c 239c---- 240 call blockin4(isupb,ibl,kbl, nbls_pnl,bl,npar) 241c---- 242c---------------------------------------------------------------- 243c memory checking : 244 call getmem(0,last11) 245 call retmem(1) 246c---------------------------------------------------------------- 247c 248ccccc do 200 ikbl=1,nblokx ! always =1 for pnl 249c 250 ikbl=1 251 call oneblock(isupb,ikbl,bl, bl(nibld),bl(nkbld), 252 * bl(nijbd),bl(nijed),bl(nklbd),bl(nkled), 253 * bl(nblok1d),bl(nqrtd), nbl2,inx, 254 * where, q4,use_q4,maxsize, 255 * lab1,lab2,lab3,lab4,fockmat,ncfunct) 256c 257c---------------------------------------------------------------- 258c memory checking : 259 call getmem(0,last12) 260 call retmem(1) 261 if(last11.ne.last12) then 262 write(6,*)'** WRONG memory alocations **' 263 write(6,*)' sup-block no=',isupb,'small-block=',ikbl 264 write(6,*)' on entry =',last11 265 write(6,*)' on exit =',last12 266 write(6,*)'** Execution has been stopped in ONESUPER **' 267* stop200 268 call errquit('texas:onsuper',0, INT_ERR) 269 endif 270c 271cc200 continue 272c---------------------------------------------------------------- 273 call retmark 274c---------------------------------------------------------------- 275 end 276c=================================================================== 277 subroutine oneblock(isupb,ikbl,bl,nibl,nkbl, nijb,nije,nklb,nkle, 278 * nblok1,nqrt, nbl2,inx, 279 * where, q4,use_q4,maxsize, 280 * lab1,lab2,lab3,lab4,fockmat,map_txs_pnl) 281c 282 implicit real*8 (a-h,o-z) 283 integer map_txs_pnl(*) ! txs to pnl basis map = ncfunct 284 logical first,firstd 285 logical use_q4 286 character*8 where 287c------------------------------------------------------ 288 common /route/ iroute 289c------------------------------------------------------ 290 common /pnl000/ xbluse,nbluse 291 common /pnl001/ ispec,ijpres2,klpres2, ijblock,klblock,iqorder 292 common /pnl002/ ncshell,ncfunct,nblock2,integ_n0 293 common /pnl003/ nqrtpnl,icstx,jcstx,kcstx,lcstx 294 common /pnl004/ isize,jsize,ksize,lsize,itxspnl 295 common /pnl005/ isblsize,isblqrts,isblpoint 296 common /map4uniq/ nij_uniqe, ij_uniqe_p, map_ij, 297 * nkl_uniqe, kl_uniqe_p, map_kl 298c------------------------------------------------------ 299 common /howmany/ ntotal,noijkl,nopres,nospec,nrimtot,nrimret 300 common /infor/ icheck,firstd,ndirect,nprint,iblok,nbeg,nend 301 common /infob/ inuc,ibas,na,nbf,nsh,ncf,ncs 302 common /number/ acc 303c 304ctime 305 common /timex/ tconv1,tconv2,ttrobs 306 common /time0/ tprec2 307 common /time1/ tpre4n,txwpq ,tassem,tamshf,tdesti 308 common /time2/ terint,ttrans 309 common /time3/ tzeroi,tspeci 310ctime 311 common /types/itype,jtype,ktype,ltype,itype1,jtype1,ktype1,ltype1 312 common /contr/ ngci,ngcj,ngck,ngcl,lci,lcj,lck,lcl,lcij,lckl 313 common /lengt/ ilen,jlen,klen,llen, ilen1,jlen1,klen1,llen1 314 common /gcont/ ngci1,ngcj1,ngck1,ngcl1,ngcd 315c 316#include "texas_lpar.fh" 317 common/obarai/ 318 * lni,lnj,lnk,lnl,lnij,lnkl,lnijkl,mmax, 319 * nqi,nqj,nqk,nql,nsij,nskl, 320 * nqij,nqij1,nsij1,nqkl,nqkl1,nskl1,ijbeg,klbeg 321c 322 common/shell/lshellt,lshelij,lshelkl,lhelp,lcas2(4),lcas3(4) 323c 324 common /memor4/ iwt0,iwt1,iwt2,ibuf,ibuf2, 325 * ibfij1,ibfij2,ibfkl1,ibfkl2, 326 * ibf2l1,ibf2l2,ibf2l3,ibf2l4,ibfij3,ibfkl3, 327 * ibf3l,issss, 328 * ix2l1,ix2l2,ix2l3,ix2l4,ix3l1,ix3l2,ix3l3,ix3l4, 329 * ixij,iyij,izij, iwij,ivij,iuij,isij 330c 331 common /memor5x/ ieab,iecd 332 common /memor5a/ iaa,ibb,icc,idd,icis,icjs,icks,icls, 333 * ixab,ixp,ixpn,ixpp,iabnia,iapb,i1apb,ifij,icij,isab, 334 * ixcd,ixq,ixqn,ixqq,icdnia,icpd,i1cpd,ifkl,ickl,iscd 335 common /memor5b/ irppq, 336 * irho,irr1,irys,irhoapb,irhocpd,iconst,ixwp,ixwq,ip1234, 337 * idx1,idx2,indx 338 common /memor5c/ itxab,itxcd,iabcd,ihabcd 339 common /memor5d/ iabnix,icdnix,ixpnx,ixqnx,ihabcdx 340 common /memor5e/ igci,igcj,igck,igcl,indgc,igcoef, 341 * icfg,jcfg,kcfg,lcfg, igcij,igckl 342 common /memor5f/ indxp 343c 344c nmr derivatives : 345 common /memor6/ ixyab,ixycd 346c 347 common /memors/ nsym,ijshp,isymm 348 common /memor1/ npard,mxsize,nblock1,nblock1_back 349c 350 dimension bl(*), inx(12,*) 351 dimension nibl(*),nkbl(*),nijb(*),nije(*),nklb(*),nkle(*) 352 dimension nblok1(2,*),nqrt(*) 353c 354 dimension fockmat(*) 355 dimension lab1(*),lab2(*),lab3(*),lab4(*) 356 dimension q4(*) 357c----------------------------------------------------------- 358c if(icheck.gt.0) return 359c----------------------------------------------------------- 360c number of ij and kl pairs (npij,npkl) 361c (zero for npkl means that the block is diagonal) 362 ibl=nibl(ikbl) 363 kbl=nkbl(ikbl) 364c 365 nijbeg=nijb(ikbl) 366 nijend=nije(ikbl) 367 npij=nijend-nijbeg+1 368c 369 nklbeg=nklb(ikbl) 370 nklend=nkle(ikbl) 371 npkl=nklend-nklbeg+1 372 if(nklend.eq.0) npkl=0 373 npklx=npkl 374 if(npkl.eq.0) npklx=npij 375c----------------------------------------------------------- 376c Get pnl_size of the ISUPB block and a pointer to its qrts. 377c 378 call get_nbls_pnl(isupb,bl(isblsize),bl(isblpoint), 379 * isbl_size,isbl_point) 380c----------------------------------------------------------- 381c get the NBLOK1 array and pnl-block-size nbls_pnl=nqrt(ikbl) 382c (for pnl the ikbl is ALWAYS equal 1). The ISYMM(ii) is zero out 383c which is needed if a sup-block is split 384c 385 call block1_pnl(isupb,isbl_size,bl(isblqrts),isbl_point,maxsize, 386 * bl(ijpres2),bl(klpres2),bl(isymm),nblok1,nqrt) 387ccc nbls_pnl=nqrt(ikbl) 388 nbls_pnl=nqrt( 1 ) 389c 390c 391 nbls=nbls_pnl 392c----------------------------------------------------------- 393c When the SPEC option is on 394c 395 nblspec=0 396 if(ispec.ne.0) then 397c 398c2002 ntotal=ncs*(ncs+1)/2 399c ntotal=ntotal*(ntotal+1)/2 400c 401c call txs_second(tim_1doq) 402c 403c needed when a given super-block is split in parts : 404c 405 call doquarts(isupb,isbl_size,bl(isblqrts),isbl_point, 406 * bl(isymm),nblspec) 407ccc * bl(ijpres2),bl(klpres2),bl(isymm),nblspec) 408c 409c nblspec= could be the same or smaller than nbls(pnl) 410c is used ONLY for performance estimation 411c----------------------------------------------------------- 412c 413c call txs_second(tim_2doq) 414c time_doqrt=time_doqrt + tim_2doq-tim_1doq 415c 416c efficiency : 417c 418 nospec=nospec+nblspec 419c 420c2002 421c xbluse=xbluse+dble(nblspec)/dble(nbls) 422c nbluse=nbluse+1 423c 424 endif 425c----------------------------------------------------------- 426c return PNL-check-run 427c if(icheck.gt.0) return 428c----------------------------------------------------------- 429c first quartet of shells : 430cnew 431 call get_ics_jcs(bl(nblock1),ibl, 1 ,ics1,jcs1) 432 call get_ics_jcs(bl(nblock1),kbl, 1 ,kcs1,lcs1) 433cnew 434 ijcs1=ics1*(ics1-1)/2+jcs1 435 klcs1=kcs1*(kcs1-1)/2+lcs1 436cnew 437c----------------------------------------------------------- 438c set up the type (common types) and length of shells and 439c information concerning contractions (commons contr,gcont) : 440c 441 call txs_shells(inx,ics1,jcs1,kcs1,lcs1) 442c----------------------------------------------------------- 443c (setup the obarai and shell commons) : 444c 445c 446 call iobara(itype1,jtype1,ktype1,ltype1,where) 447c----------------------------------------------------------- 448c for memory handling 449c 450 nfumax=nfu(mmax) 451 mmax1=mmax-1 452 if(mmax1.eq.0) mmax1=1 453c 454c98 if(nsij.ge.nskl) then 455c98 nhabcd=nkl_uniqe 456c98 nfha=nfumax*nhabcd*lckl 457c98 else 458c98 nhabcd=nij_uniqe 459c98 nfha=nfumax*nhabcd*lcij 460c98 endif 461c98 462c 463c this is needed ONLY for TRACY's recursive 464c and NOW do it ALWAYS as follows : 465c 466 nhabcd=nkl_uniqe 467 nfha=nfumax*nhabcd*lckl 468c----------------------------------------------------------- 469c memory handling 470c 471c 1: for pairs precalculations (20 quantities) 472c (for whole block of contracted quartets and 473c all primitive quartets belonging to this block) 474c 475c 2: and quartets precalculations (12 quantities) 476c (for whole block of contracted quartets and 477c one primitive quartet ) 478c 479c 3: and for individual shells (8 quantities) 480c ( aa, bb, cc, dd - exponents ) 481c ( cis,cjs,cks,cls -contr coef.) 482c 483c nbls = a full size of the pnl-requested super-block 484c----------------------------------------------------------- 485 IF( iroute.eq.1 ) THEN 486 call memo5a_1(nij_uniqe ,mmax1) 487 call memo5b_1(nkl_uniqe,mmax1) 488 call memo5c_1(bl, nbls,mmax1,nij_uniqe,nkl_uniqe,nfha,nfumax) 489c 490c 50 or 56 calls of getmem 491 ncalls=50 492 if(ngcd.gt.1) ncalls=56 493 ELSE 494 call memo5a_2(nij_uniqe ,mmax1) 495 call memo5b_2(nkl_uniqe,mmax1) 496 call memo5c_2(nbls,mmax1,nij_uniqe,nkl_uniqe,nfumax) 497c 498c 43 or 47 calls of getmem plus 3 calls for where='forc'or'hess' 499c 500 ncalls=43 501 if(ngcd.gt.1) ncalls=47 502 if(where.eq.'forc' .or. where.eq.'hess') ncalls=ncalls+3 503 ENDIF 504c----------------------------------------------------------- 505c Perform pre-calculations for pairs of contracted sells 506c 507 call txs_second(tprec2b) 508c 509 IF( iroute.eq.1 ) THEN 510 call precalc2_1(isupb,bl,mmax,mmax1,nhabcd,nfumax,nbl2,nbls, 511 * inx,ibl,kbl,npkl) 512 ELSE 513 call precalc2_2(isupb,bl,mmax,mmax1, nfumax, nbl2,nbls, 514 * inx,ibl,kbl,npkl) 515 ENDIF 516c----------------------------------------------------------- 517c for NMR derivatives 518c 519c NOT READY for PNL : inuc, ibas refer to dbl_mb() 520c 521cnmr if(where.eq.'shif') then 522cnmr call memo6(npij,npklx) ! 2 calls of getmem 523cnmr call getmem(npij ,ijcent) 524cnmr call getmem(npklx,klcent) 525cnmr call precal2d(bl(inuc+1),iis,jjs,inx, npij,npklx,npkl, 526cnmr * ibl,kbl,ijbl,nbl2,nijbeg,nijend,nklbeg,nklend, 527cnmr * bl(ixyab),bl(ixycd), 528cnmr * bl(isymm),bl(ijcent),bl(klcent) ) 529cnmr call retmem(2) 530cnmr endif 531c----------------------------------------------------------- 532c neglect some of contracted quartets on the cont. shell level : 533c isym(ijkl)=0 or 1 , 0 means that ijkl quartet is neglected. 534c 535c (neglect on the contected shell level is not performed for PNL) 536c 537c NO call to the CSHNEG routine !!!! 538c----------------------------------------------------------- 539c At this moment some of c.s.quartets may not be present if this run 540c is with MORE_INT=.true. which means that this super-block was split 541c and some of its quartets were calculated in a previous call. 542c 543c Thus, reduce the block size from nbls ---> nblsp (present) 544c (only present c.s.quartets are considered) 545c 546c define idx1 (ijpar <--- ijkl quartets mapping) 547c and idx2 (klpar <--- ijkl quartets mapping) 548c 549 call indexp_pnl(isupb,isbl_size,bl(isblqrts),isbl_point, 550 * bl(map_ij),bl(map_kl), 551 * bl(idx1),bl(idx2),bl(isymm),bl(indxp),nblsp) 552c 553c----------------------------------------------------------- 554 call txs_second(tprec2e) 555 tprec2=tprec2+tprec2e-tprec2b 556c----------------------------------------------------------- 557c FROM hereafter the block-size is reduced nbls-->nblsp 558c 559 if(nblsp.eq.0) then 560c write(8,*)' From oneblock block=',ikbl,'ngcd=',ngcd 561c write(8,*)' nbls--->nblsp ',nbls,nblsp 562c ncalls-1 comes from the fact that no memory has been res. 563c for buf in this case : 564 ncalls=ncalls-1 565 go to 110 566 endif 567c 568 nopres=nopres+nblsp 569cccc nblsrem=nbls 570 nbls=nblsp 571c----------------------------------------------------------- 572c calculate special integrals with mmax<=2 573c (ssss), (psss),(lsss) etc. 574c 575 if(mmax.le.2) then 576 call txs_second(tspec1) 577 if(icheck.eq.0) then 578 IF( iroute.eq.1 ) THEN 579cnew call erintsp_1(isupb,bl,first,nbls,acc, ikbl,npij,npklx,npkl) 580 call erintsp_1(isupb,bl,first,nbls,acc, ikbl, npkl) 581 ELSE 582cnew call erintsp_2(isupb,bl,first,nbls,acc, ikbl,npij,npklx,npkl) 583 call erintsp_2(isupb,bl,first,nbls,acc, ikbl, npkl) 584 ENDIF 585 endif 586 call txs_second(tspec2) 587 tspeci=tspeci+tspec2-tspec1 588 if( first ) go to 110 589 go to 120 590 endif 591c 592c----------------------------------------------------------- 593c calculate integrals with MMAX > 2 594c 595ctime 596 call txs_second(teri) 597c 598 if(icheck.eq.0) then 599 IF( iroute.eq.1 ) THEN 600cnew call erinteg_1(isupb,bl,first,nbls,acc,ikbl,npij,npklx,npkl, 601c98 call erinteg_1(isupb,bl,first,nbls,acc,ikbl, npkl, 602c98 * lobsa,immax,kmmax, where) 603 call erinteg_1(isupb,bl,first,nbls,acc,ikbl, npkl, 604 * where) 605 ELSE 606cnew call erinteg_2(isupb,bl,first,nbls,acc,ikbl,npij,npklx,npkl, 607c98 call erinteg_2(isupb,bl,first,nbls,acc,ikbl, npkl, 608c98 * lobsa,immax,kmmax, where) 609 call erinteg_2(isupb,bl,first,nbls,acc,ikbl, npkl, 610 * where) 611 ENDIF 612 endif 613c 614c----------------------------------------------------------- 615 call txs_second(terie) 616 terint=terint+(terie-teri) 617c 618 if( first ) go to 110 619c 620c----------------------------------------------------------- 621c transformatin d6->d5, f10->f7, g15->g9, h21->h11, i28->i13 622c 623c do we need to transform ? 624c 625 lenall=ilen+jlen+klen+llen 626 lenall1=ilen1+jlen1+klen1+llen1 627 if(lenall1-lenall.ne.0 ) then 628c>>>> if( cart_2_sphe ) then 629 call txs_second(trans1) 630c ------------------------------------------------------ 631c get spherical harmonics PNL---types : ityps - ltyps : 632c using Texas cartesian types itype1-ltype1 633c 634c cartesian basis set shell sizes : ilen1 - llen1 : 635c spherical harmonics shell sizes : ilen - llen : 636c spherical harmonics PNL---types : ityps - ltyps : 637c 638 call get_spher_pnl_type(itype1,ityps) 639 call get_spher_pnl_type(jtype1,jtyps) 640 call get_spher_pnl_type(ktype1,ktyps) 641 call get_spher_pnl_type(ltype1,ltyps) 642c ------------------------------------------------------- 643c 644c 645 mnbls=nbls 646 nbuf=ibuf 647 if(where.eq.'shif') then 648c --- nmr derivatives ---- 649 mnbls=6*nbls 650 nbuf=ibuf+ngcd*nbls*lnijkl 651 endif 652 if(where.eq.'forc') then 653c --- gradient derivatives ---- 654 mnbls=9*nbls 655 nbuf=ibuf 656 endif 657 if(where.eq.'hess') then 658c --- hessian derivatives ---- 659 mnbls=45*nbls ! for second derivatives only 660cccc mnbls=54*nbls ! for second and first derivatives together 661 nbuf=ibuf 662 endif 663c 664 incrt=mnbls*lnijkl 665c 666 do 130 iqu=1,ngcd ! over gen.cont. 667 jbuf=nbuf+(iqu-1)*incrt 668 call transfor(bl,mnbls,jbuf, 669 * ityps , jtyps, ktyps, ltyps, 670 * ilen1 , jlen1, klen1, llen1, 671 * ilen , jlen , klen , llen ) 672 130 continue 673 125 continue 674 call txs_second(trans2) 675 ttrans=ttrans+trans2-trans1 676 endif 677c 678c end of transformation 679c----------------------------------------------------------- 680c 681 120 continue 682c----------------------------------------------------------- 683 if(icheck.gt.0) go to 110 684c----------------------------------------------------------- 685 call txs_second(tdest) 686c-------------------------------------------------------------- 687c Final use of integrals : 688c 689c For Battelle PNL : special request (SPEC option in the input) 690c 691c 5. put integrals from a selected quartet of shells into 692c a buffer 693c (when where='buff' - call destbuf or destbul) 694c (when where='forc' - call destduf or destdul) 695c 696c-------------------------------------------------------------- 697c 698 if ( where.eq.'buff') then 699c-5> put two-el.int. into the buffer - only a set of integrals 700c fockmat is now the buffer for integrals 701c 702 if(ispec.eq.1) then 703c integrals without labels (all of them): 704 call destbuf(ikbl,nbls,nblok1, ncs,inx, 705 * bl(ibuf),fockmat,bl(itxspnl), q4,use_q4, 706 * bl(icfg),bl(jcfg),bl(kcfg),bl(lcfg),ngcd,lnijkl, 707 * bl(indxp) ,bl(isymm)) 708 endif 709 if(ispec.eq.2) then 710c integrals (non-zeros) with labels 711 call destbul(ikbl,nbls,nblok1, ncs,inx, 712 * bl(ibuf),fockmat, lab1,lab2,lab3,lab4, q4,use_q4, 713 * bl(icfg),bl(jcfg),bl(kcfg),bl(lcfg),ngcd,lnijkl, 714 * bl(indxp) ,bl(isymm),bl(iqorder),map_txs_pnl) 715 endif 716 endif 717c 718 if ( where.eq.'forc') then 719c gradient derivatives : 720c 721 if(ispec.eq.1) then 722c integrals without labels (all of them): 723 call destduf(ikbl,nbls,nblok1, ncs,inx, 724 * bl(ibuf),fockmat,bl(itxspnl), q4,use_q4, 725 * bl(icfg),bl(jcfg),bl(kcfg),bl(lcfg),ngcd,lnijkl, 726 * bl(indxp) ,bl(isymm), bl(iqorder)) 727 endif 728 if(ispec.eq.2) then 729c integrals (non-zeros) with labels 730 call destdul(ikbl,nbls,nblok1, ncs,inx, 731 * bl(ibuf),fockmat, lab1,lab2,lab3,lab4, q4,use_q4, 732 * bl(icfg),bl(jcfg),bl(kcfg),bl(lcfg),ngcd,lnijkl, 733 * bl(indxp) ,bl(isymm),bl(iqorder),map_txs_pnl) 734 endif 735 endif 736c 737 if ( where.eq.'hess') then 738c hessian derivatives : 739c 740 if(ispec.eq.1) then 741c integrals without labels (all of them): 742 call desthuf(ikbl,nbls,nblok1, ncs,inx, 743 * bl(ibuf),fockmat,bl(itxspnl), q4,use_q4, 744 * bl(icfg),bl(jcfg),bl(kcfg),bl(lcfg),ngcd,lnijkl, 745 * bl(indxp) ,bl(isymm), bl(iqorder)) 746 endif 747 if(ispec.eq.2) then 748c integrals (non-zeros) with labels 749 call desthul(ikbl,nbls,nblok1, ncs,inx, 750 * bl(ibuf),fockmat, lab1,lab2,lab3,lab4, q4,use_q4, 751 * bl(icfg),bl(jcfg),bl(kcfg),bl(lcfg),ngcd,lnijkl, 752 * bl(indxp) ,bl(isymm),bl(iqorder),map_txs_pnl) 753 endif 754 endif 755c----------------------------------------------------------- 756c 757 call txs_second(tdeste) 758 tdesti=tdesti+(tdeste-tdest) 759c 760 110 continue 761c 762c----------------------------------------------------------- 763c release memory at the end of a given block : 764c 765c nmr deriv: 766c 767 if(where.eq.'shif') ncalls=ncalls+2 768c 769 if(icheck.eq.0) then 770 call retmem(ncalls+1) ! +1 to release ibuf 771 else 772 call retmem(ncalls) 773 endif 774c---------------------------------------- 775 end 776c=================================================================== 777 subroutine erinteg_1(isupb,bl,first,nbls,acc,ikbl,npkl, where) 778c 779c character variable WHERE is needed for derivatives only 780c 781 implicit real*8 (a-h,o-z) 782 logical first 783 character*8 where 784c 785 common /primij/ iabprim, ijdim ,ijpar1 786 common /primkl/ kabprim, kldim ,klpar1 787c 788 common /howmany/ ntotal,noijkl,nopres,nospec,nrimtot,nrimret 789ctime 790 common /timex/ tconv1,tconv2,ttrobs 791 common /time1/ tpre4n,txwpq ,tassem,tamshf,tdesti 792 common /time3/ tzeroi,tspeci 793 common /time4/ tderiv 794ctime 795c 796 common /types/itype,jtype,ktype,ltype,itype1,jtype1,ktype1,ltype1 797 common /contr/ ngci,ngcj,ngck,ngcl,lci,lcj,lck,lcl,lc12,lc34 798 common /gcont/ ngci1,ngcj1,ngck1,ngcl1,ngcd 799c 800#include "texas_lpar.fh" 801c 802 common/obarai/ 803 * lni,lnj,lnk,lnl,lnij,lnkl,lnijkl,mmax, 804 * nqi,nqj,nqk,nql,nsij,nskl, 805 * nqij,nqij1,nsij1,nqkl,nqkl1,nskl1,ijbeg,klbeg 806c 807 common /memor3/ nblok1d 808c 809 common /memor4/ iwt0,iwt1,iwt2,ibuf,ibuf2, 810 * ibfij1,ibfij2,ibfkl1,ibfkl2, 811 * ibf2l1,ibf2l2,ibf2l3,ibf2l4,ibfij3,ibfkl3, 812 * ibf3l,issss, 813 * ix2l1,ix2l2,ix2l3,ix2l4,ix3l1,ix3l2,ix3l3,ix3l4, 814 * ixij,iyij,izij, iwij,ivij,iuij,isij 815c 816 common /memor5x/ ieab,iecd 817 common /memor5a/ iaa,ibb,icc,idd,icis,icjs,icks,icls, 818 * ixab,ixp,ixpn,ixpp,iabnia,iapb,i1apb,ifij,icij,isab, 819 * ixcd,ixq,ixqn,ixqq,icdnia,icpd,i1cpd,ifkl,ickl,iscd 820 common /memor5b/ irppq, 821 * irho,irr1,irys,irhoapb,irhocpd,iconst,ixwp,ixwq,ip1234, 822 * idx1,idx2,indx 823c 824 common /memor5c/ itxab,itxcd,iabcd,ihabcd 825 common /memor5d/ iabnix,icdnix,ixpnx,ixqnx,ihabcdx 826c new for grad. derivatives: 827 common /memor5dd/ iaax,ibbx,iccx 828 common /memor5e/ igci,igcj,igck,igcl,indgc,igcoef, 829 * icfg,jcfg,kcfg,lcfg, igcij,igckl 830 common /memor5f/ indxp 831c nmr der. 832 common /memor6/ ixyab,ixycd 833c 834 common /memors/ nsym,ijshp,isymm 835c 836 common /pnl001/ ispec,ijpres2,klpres2, ijblock,klblock,iqorder 837 common /pnl005/ isblsize,isblqrts,isblpoint 838 common /map4uniq/ nij_uniqe, ij_uniqe_p, map_ij, 839 * nkl_uniqe, kl_uniqe_p, map_kl 840c 841 dimension bl(*) 842c----------------------------------------------------------------- 843c memory handling for a block 844c (independent of contractions) 845c reserves memory for trobsa and assemble 846c 847 call memo4a(bl, nbls, l11,l12,mem2,igmcnt) 848c 849 nfumax=nfu(mmax) 850 mmax1=mmax-1 851 nfumax1=nfumax 852 narray=1 853 nqijr=nqij 854 nqklr=nqkl 855 if(where.eq.'shif') nfumax1=nfu(mmax1) 856 if(where.eq.'forc') then 857 nfumax1=nfu(mmax1) 858 narray=4 859 nqij=nqij-1 860 nqkl=nqkl-1 861 if(nqij.le.0) nqij=1 862 if(nqkl.le.0) nqkl=1 863 endif 864 if(where.eq.'hess') then 865 nfumax1=nfu(mmax1-1) 866 narray=10 867 nqij=nqij-2 868 nqkl=nqkl-2 869 if(nqij.le.0) nqij=1 870 if(nqkl.le.0) nqkl=1 871 endif 872c----------------------------------------------------------------- 873 first=.true. 874c----------------------------------------------------------------- 875 igcoet=1 876 if(ngcd.gt.1) then 877 ngcij=ngci1*ngcj1 878 ngckl=ngck1*ngcl1 879 call getmem(ngcij*nbls,igcij) 880 call getmem(ngckl*nbls,igckl) 881 call getmem(ngcd*nbls ,igcoet) 882 endif 883c 884 if(where.eq.'forc' .or. where.eq.'hess') then 885 call getmem(nbls,iaax) 886 call getmem(nbls,ibbx) 887 call getmem(nbls,iccx) 888 endif 889c----------------------------------------------------------------- 890c get pnl-size of the super-block and the pointer to its quartets: 891c 892 call get_nbls_pnl(isupb,bl(isblsize),bl(isblpoint), 893 * isbl_size,isbl_point) 894c 895 call getmem(isbl_size,map_n0) 896 call get_nbls_n0(isbl_size,isbl_point,bl(isblqrts), 897 * bl(map_n0), nbls_n0) 898c----------------------------------------------------------------- 899c 900c loop over contraction : 901c 902 lcij=0 903 do 12 l1=1,lci 904 do 12 l2=1,lcj 905 if(ngcd.gt.1) then 906 call gcparij(nbls, bl(idx1),nij_uniqe, 907 * l1,l2,ngci1,ngcj1,ngcij, 908 * bl(igci),bl(igcj), bl(igcij)) 909 endif 910 lcij=lcij+1 911 lckl=0 912 nblsr=0 913 do 34 l3=1,lck 914 do 34 l4=1,lcl 915 if(ngcd.gt.1) then 916 call gcparkl(nbls, bl(idx2),nkl_uniqe, 917 * l3,l4,ngck1,ngcl1,ngckl, 918 * bl(igck),bl(igcl), bl(igckl)) 919 endif 920 lckl=lckl+1 921c 922c calculate : rppq,rho, rysx,rhoapb,rhocpd : 923c 924ctime 925 call txs_second(t4b) 926c 927c2002 call prec4neg_1(isbl_size,isbl_point,bl(isblqrts),npkl, 928c 929 call prec4neg_1(nbls_n0,bl(map_n0), lcij,lckl,lc12,lc34, 930 * bl(ieab),bl(iecd),bl(iapb),bl(icpd),bl(icij),bl(ickl), 931 * bl(ixp),bl(ixq), 932 * bl(map_ij),bl(map_kl),nij_uniqe,nkl_uniqe, 933 * bl(irppq),bl(irhoapb),bl(irhocpd),bl(irys),bl(iconst), 934 * nbls1,bl(indx) ) 935c 936 call txs_second(t4e) 937 tpre4n=tpre4n+(t4e-t4b) 938c 939 nrimret=nrimret+nbls1 940 nrimtot=nrimtot+nbls 941 if(nbls1.eq.0) go to 34 942c 943c note : 944c from here a given block (nbls) was reduced to (nbls1) *** 945c for given l1,l2,l3,l4 !!! 946c and 947c rho, rppq, rhoapb, rhocpd, rys, and const 948c have dimension nbls1 (not (nbls) ) 949c----------------------------------------------------------------- 950c substitute zero for all needed bufors for quartets 951c which do not appear first time after neglect 952c 953 if(first) then 954ctime 955 call txs_second(tzer) 956 nblsnot=nbls-nbls1 957 if(nblsnot.ne.0) then 958ckw call zeroint(bl,nbls,nbls1,nblsnot,lnij,lnkl,indx,ngcd) 959 call zeroint(bl,nbls,nbls1,nblsnot,lnij,lnkl,indx,ngcd,narray) 960 endif 961 call txs_second(tzero) 962 tzeroi=tzeroi+(tzero-tzer) 963 endif 964c----------------------------------------------------------------- 965c calculate products of contraction coefficients 966c for general contracted basis set 967c 968 if(ngcd.GT.1) then 969 call gcqijkl(nbls,nbls1, bl(indx), 970 * ngci1,ngcj1,ngck1,ngcl1,ngcd, 971 * bl(indgc),bl(igcoef), 972 * bl(igcij),ngcij, bl(igckl),ngckl) 973 endif 974c----------------------------------------------------------------- 975 call txs_second(txwpqb) 976c 977 call xwpq_1(nbls1,bl(ixwp),bl(ixwq),bl(ip1234), 978 * lcij,lckl,bl(idx1),bl(idx2),bl(indx), 979 * nij_uniqe,nkl_uniqe, 980 * bl(irppq),bl(ixp),bl(ixq),bl(ixpp),bl(ixqq), 981 * bl(itxab),bl(itxcd), 982 * bl(iapb),bl(i1cpd),bl(icpd),bl(i1apb) ) 983c 984c----------------------------------------------------------------- 985c calculate abcd(i)=apb/cpd or vice versa depending on stability : 986c 987 call abcd_1(nbls1, l1,l2, lcij,lckl, 988 * bl(idx1),bl(idx2),bl(indx),nij_uniqe,nkl_uniqe, 989 * bl(iabcd), bl(iapb),bl(i1cpd),bl(icpd),bl(i1apb), 990 * bl(iaa),bl(ibb),bl(iconst),bl(igcoef),ngcd) 991c----------------------------------------------------------------- 992c convert: abnia,cdnia, xpn,xqn, habcd 993c from pair to quartet (nbls1) quantities 994c for trobsa 995c 996cxxxxij 997ctime 998 call txs_second(tc1xb) 999 txwpq=txwpq+(tc1xb-txwpqb) 1000c 1001 if(nbls1.ne.nbls .or. nblsr.ne.nbls) then 1002 call conv1x_1(nbls1,mmax1,nij_uniqe ,lcij, bl(idx1),bl(indx), 1003 * bl(iabnia), bl(ixpn), bl(iabnix),bl(ixpnx) ) 1004 if(where.eq.'forc' .or. where.eq.'hess') then 1005 call conv1der_1(nbls1,nij_uniqe,l1,bl(idx1),bl(indx), 1006 * bl(iaa),bl(iaax) ) 1007 call conv1der_1(nbls1,nij_uniqe,l2,bl(idx1),bl(indx), 1008 * bl(ibb),bl(ibbx) ) 1009 endif 1010 endif 1011cxxxxkl 1012 call conv1x_1(nbls1,mmax1,nkl_uniqe,lckl, bl(idx2),bl(indx), 1013 * bl(icdnia), bl(ixqn), bl(icdnix),bl(ixqnx) ) 1014 if(where.eq.'forc' .or. where.eq.'hess') then 1015 call conv1der_1(nbls1,nkl_uniqe,l3,bl(idx2),bl(indx), 1016 * bl(icc),bl(iccx) ) 1017 endif 1018ctime 1019 call txs_second(tc1xe) 1020 tconv1=tconv1+(tc1xe-tc1xb) 1021cxxxx2 1022c98 if(nsij.ge.nskl) then 1023 call conv2x(nbls1,nfumax1,nkl_uniqe,lckl, bl(idx2),bl(indx), 1024 * bl(ihabcd),nfumax, bl(ihabcdx) ) 1025c98 else 1026c98 if(nbls1.ne.nbls .or. nblsr.ne.nbls) then 1027c98 call conv2x(nbls1,nfumax1,nij_uniqe,lcij, bl(idx1),bl(indx), 1028c98 * bl(ihabcd),nfumax, bl(ihabcdx) ) 1029c98 endif 1030c98 endif 1031ctime 1032 call txs_second(tc2xe) 1033 tconv2=tconv2+(tc2xe-tc1xe) 1034c 1035c98 call trobsa(bl,nbls1,l11,l12,mem2,immax,kmmax,lobsa) 1036 call trobsa_1(bl,nbls1,l11,l12,mem2) 1037c 1038ctime 1039 call txs_second(ttrob) 1040 ttrobs=ttrobs+(ttrob-tc2xe) 1041cxxxxxx 1042 if(ngcd.eq.1) then 1043 call assemblx(bl,first,nbls,nbls1,lnij,lnkl, 1044 * l1,l2,l3,l4,lcij,lckl,nij_uniqe,nkl_uniqe) 1045cnew * l1,l2,l3,l4,lcij,lckl,npij,npklx) 1046 else 1047c98 call gcqijkl(nbls,nbls1, bl(indx), 1048c98 * ngci1,ngcj1,ngck1,ngcl1,ngcd, 1049c98 * bl(indgc),bl(igcoef), 1050c98 * bl(igcij),ngcij, bl(igckl),ngckl) 1051c 1052c transpose gcoef(ngcd,nbls) into gcoet(nbls,ngcd) 1053 call trspmo(bl(igcoef),ngcd,bl(igcoet),nbls) 1054c 1055 call assemblg(bl,first,nbls,nbls1,lnij,lnkl,ngcd,igcoet) 1056 endif 1057ctime 1058 call txs_second(tasse) 1059 tassem=tassem+(tasse-ttrob) 1060cxxxxx 1061c 1062 nblsr=nbls1 1063c 1064 34 continue 1065 12 continue 1066c 1067c grad & hess derivatives: 1068c 1069 nqij=nqijr 1070 nqkl=nqklr 1071c 1072c----------------------------------------------------------------- 1073c release memory : 1074c 1075 call retmem(1) ! allocated for map_n0 1076c 1077 if(where.eq.'forc' .or. where.eq.'hess') call retmem(3) 1078c 1079 if(ngcd.gt.1) call retmem(3) 1080c 1081 if( first ) then 1082 call retmem(igmcnt) 1083 return 1084 endif 1085c 1086c release memory from trobsa 1087c 1088 call retmem(3) 1089c 1090c------------------------------------ 1091c for NMR, GRADIENT and HESSIAN derivatives : 1092c 1093 lnijr=lnij 1094 lnklr=lnkl 1095c------------------------------------ 1096c reserves memory for amshift 1097c 1098 call memo4b(bl, nbls,igmcnt1) 1099c 1100c------------------------------------ 1101c for NMR derivatives : 1102c 1103 if(where.eq.'shif') then 1104 call txs_second(tderb) 1105c 1106c return to the original values of nsij,nskl and mmax : 1107c 1108 ijderiv=1 1109 klderiv=1 1110 call iobarb(ijderiv,klderiv) 1111c 1112c construct nmr/giao derivatives : (i+j,s|k+l,s)(x,y,z) 1113c 1114 call shift_der(nbls,lnijr,lnklr,nij_uniqe,nkl_uniqe,ngcd, 1115 * idx1,idx2, ixab,ixcd,ixyab,ixycd) 1116c 1117 call txs_second(tdere) 1118 tderiv=tderiv+tdere-tderb 1119 endif 1120c------------------------------------ 1121c for GRADIENT derivatives : 1122c 1123 if(where.eq.'forc') then 1124 call txs_second(tderb) 1125c 1126c return to the original values of nsij,nskl and mmax : 1127c 1128 ijderiv=1 1129 klderiv=1 1130 call iobarb(ijderiv,klderiv) 1131c 1132c construct gradient derivatives : (i+j,s|k+l,s)(x,y,z) 1133c 1134 call force_der(bl,nbls,lnijr,lnklr,nij_uniqe,ngcd,idx1,ixab) 1135c 1136 call txs_second(tdere) 1137 tderiv=tderiv+tdere-tderb 1138 endif 1139c------------------------------------ 1140c for HESSIAN derivatives : 1141c 1142 if(where.eq.'hess') then 1143 call txs_second(tderb) 1144c 1145c return to the original values of nsij,nskl and mmax : 1146c 1147 ijderiv=2 1148 klderiv=2 1149 call iobarb(ijderiv,klderiv) 1150c 1151c construct hessian derivatives : (i+j,s|k+l,s)(x,y,z) 1152c 1153 call hessian_der(bl,nbls,lnijr,lnklr,nij_uniqe,ngcd,idx1,ixab) 1154c 1155 call txs_second(tdere) 1156 tderiv=tderiv+tdere-tderb 1157 endif 1158c------------------------------------ 1159c* shift angular momentum (a->b, c->d) : 1160ctime 1161 call txs_second(tamsb) 1162cnew call amshift(bl,nbls,lnij,lnkl,npij,npklx,ngcd) 1163 call amshift(bl,nbls,lnij,lnkl,nij_uniqe,nkl_uniqe,ngcd) 1164 call txs_second(tamse) 1165 tamshf=tamshf+(tamse-tamsb) 1166c----- 1167 call retmem(igmcnt+igmcnt1-3) 1168c----- 1169c 1170 end 1171c=================================================================== 1172 subroutine erintsp_1(isupb,bl,first,nbls,acc,ikbl,npkl) 1173 implicit real*8 (a-h,o-z) 1174 logical first 1175c 1176 common /primij/ iabprim, ijdim ,ijpar1 1177 common /primkl/ kabprim, kldim ,klpar1 1178 common /howmany/ ntotal,noijkl,nopres,nospec,nrimtot,nrimret 1179 common /types/itype,jtype,ktype,ltype,itype1,jtype1,ktype1,ltype1 1180 common /contr/ ngci,ngcj,ngck,ngcl,lci,lcj,lck,lcl,lc12,lc34 1181 common /gcont/ ngci1,ngcj1,ngck1,ngcl1,ngcd 1182#include "texas_lpar.fh" 1183 COMMON/SHELL/LSHELLT,LSHELIJ,LSHELKL,LHELP,LCAS2(4),LCAS3(4) 1184 common/obarai/ 1185 * lni,lnj,lnk,lnl,lnij,lnkl,lnijkl,mmax, 1186 * nqi,nqj,nqk,nql,nsij,nskl, 1187 * nqij,nqij1,nsij1,nqkl,nqkl1,nskl1,ijbeg,klbeg 1188 common /memor3/ nblok1d 1189 common /memor4/ iwt0,iwt1,iwt2,ibuf,ibuf2, 1190 * ibfij1,ibfij2,ibfkl1,ibfkl2, 1191 * ibf2l1,ibf2l2,ibf2l3,ibf2l4,ibfij3,ibfkl3, 1192 * ibf3l,issss, 1193 * ix2l1,ix2l2,ix2l3,ix2l4,ix3l1,ix3l2,ix3l3,ix3l4, 1194 * ixij,iyij,izij, iwij,ivij,iuij,isij 1195 common /memor5x/ ieab,iecd 1196 common /memor5a/ iaa,ibb,icc,idd,icis,icjs,icks,icls, 1197 * ixab,ixp,ixpn,ixpp,iabnia,iapb,i1apb,ifij,icij,isab, 1198 * ixcd,ixq,ixqn,ixqq,icdnia,icpd,i1cpd,ifkl,ickl,iscd 1199 common /memor5b/ irppq, 1200 * irho,irr1,irys,irhoapb,irhocpd,iconst,ixwp,ixwq,ip1234, 1201 * idx1,idx2,indx 1202c 1203 common /memor5c/ itxab,itxcd,iabcd,ihabcd 1204 common /memor5d/ iabnix,icdnix,ixpnx,ixqnx,ihabcdx 1205 common /memor5e/ igci,igcj,igck,igcl,indgc,igcoef, 1206 * icfg,jcfg,kcfg,lcfg, igcij,igckl 1207 common /memor5f/ indxp 1208 common /memors/ nsym,ijshp,isymm 1209c 1210 common /pnl001/ ispec,ijpres2,klpres2, ijblock,klblock,iqorder 1211 common /pnl005/ isblsize,isblqrts,isblpoint 1212 common /map4uniq/ nij_uniqe, ij_uniqe_p, map_ij, 1213 * nkl_uniqe, kl_uniqe_p, map_kl 1214c 1215 dimension bl(*) 1216c----------------------------------------------------------------- 1217c memory handling for a block 1218c (independent of contractions) 1219c reserves memory for trobsa and assemble 1220c 1221 call memo4a(bl, nbls, l11,l12,mem2,igmcnt) 1222c----------------------------------------------------------------- 1223 first=.true. 1224c----------------------------------------------------------------- 1225 igcoet=1 1226 if(ngcd.gt.1) then 1227 ngcij=ngci1*ngcj1 1228 ngckl=ngck1*ngcl1 1229 call getmem(ngcij*nbls,igcij) 1230 call getmem(ngckl*nbls,igckl) 1231 call getmem(ngcd*nbls ,igcoet) 1232 endif 1233c----------------------------------------------------------------- 1234c get pnl-size of the super-block and the pointer to its quartets: 1235c 1236 call get_nbls_pnl(isupb,bl(isblsize),bl(isblpoint), 1237 * isbl_size,isbl_point) 1238c 1239 call getmem(isbl_size,map_n0) 1240 call get_nbls_n0(isbl_size,isbl_point,bl(isblqrts), 1241 * bl(map_n0), nbls_n0) 1242c----------------------------------------------------------------- 1243c 1244c* loop over contraction : 1245c 1246 lcij=0 1247 do 12 l1=1,lci 1248 do 12 l2=1,lcj 1249 if(ngcd.gt.1) then 1250 call gcparij(nbls, bl(idx1),nij_uniqe, 1251 * l1,l2,ngci1,ngcj1,ngcij, 1252 * bl(igci),bl(igcj), bl(igcij)) 1253 endif 1254 lcij=lcij+1 1255 lckl=0 1256 do 34 l3=1,lck 1257 do 34 l4=1,lcl 1258 if(ngcd.gt.1) then 1259 call gcparkl(nbls, bl(idx2),nkl_uniqe, 1260 * l3,l4,ngck1,ngcl1,ngckl, 1261 * bl(igck),bl(igcl), bl(igckl)) 1262 endif 1263 lckl=lckl+1 1264c 1265c* calculate : rppq,rho, rysx,rhoapb,rhocpd : 1266c 1267 call precspec_1(nbls_n0,bl(map_n0), lcij,lckl,lc12,lc34, 1268 * bl(idx1),bl(idx2), 1269 * bl(ieab),bl(iecd), bl(iapb),bl(icpd),bl(icij),bl(ickl), 1270 * bl(ixp),bl(ixq),bl(i1apb),bl(i1cpd),bl(itxab),bl(itxcd), 1271 * bl(map_ij),bl(map_kl),nij_uniqe,nkl_uniqe, 1272 * bl(irys),bl(iconst),bl(ixwp),bl(ixwq), nbls1,bl(indx) ) 1273c 1274 nrimret=nrimret+nbls1 1275 nrimtot=nrimtot+nbls 1276 if(nbls1.eq.0) go to 34 1277c 1278c note : 1279c* from here a given block (nbls) was reduced to (nbls1) *** 1280c* for given l1,l2,l3,l4 !!! 1281c* and 1282c* rho, rppq, rhoapb, rhocpd, rys, and const 1283c have dimension nbls1 (not (nbls) ) 1284c****************************************************** 1285c* substitute zero for all needed bufors for quartets 1286c* which do not appear first time after neglect 1287c 1288 if(first) then 1289 nblsnot=nbls-nbls1 1290 if(nblsnot.ne.0) then 1291ckw call zeroint(bl,nbls,nbls1,nblsnot,lnij,lnkl,indx,ngcd) 1292 call zeroint(bl,nbls,nbls1,nblsnot,lnij,lnkl,indx,ngcd,1) 1293 endif 1294 endif 1295c 1296c******** 1297c 1298c* special code for (ss ss) and (xs ss), (sx ss), (ss xs), (ss sx) 1299c* itegrals ( x= p or l ) 1300c 1301 if(ngcd.eq.1) then 1302 call specase_1(bl,first,nbls,nbls1,bl(indx),bl(idx1),bl(idx2), 1303 * nij_uniqe,nkl_uniqe,l1,l2,l3,l4, 1304 * bl(icis),bl(icjs),bl(icks),bl(icls), 1305 * bl(ibuf),bl(ibuf2), 1306 * bl(iconst),bl(irys),bl(ixwp),bl(ixwq),bl(irr1) ) 1307 else 1308 call gcqijkl(nbls,nbls1, bl(indx), 1309 * ngci1,ngcj1,ngck1,ngcl1,ngcd, 1310 * bl(indgc),bl(igcoef), 1311 * bl(igcij),ngcij, bl(igckl),ngckl) 1312c 1313c transpose gcoef(ngcd,nbls) into gcoet(nbls,ngcd) 1314 call trspmo(bl(igcoef),ngcd,bl(igcoet),nbls) 1315c 1316 call specasg(bl,first,nbls,nbls1, bl(indx),bl(idx1),bl(idx2), 1317 * bl(ibuf),bl(ibuf2), 1318 * bl(iconst),bl(irys),bl(ixwp),bl(ixwq), 1319 * ngcd,bl(indgc),bl(igcoet),lnijkl) 1320 endif 1321c 1322 34 continue 1323 12 continue 1324c 1325c----------------------------------------------------------------- 1326c release memory 1327c 1328 call retmem(1) ! allocated for map_n0 1329 if(ngcd.gt.1) call retmem(3) 1330c---- 1331 call retmem(igmcnt) 1332c 1333 end 1334c=================================================================== 1335 subroutine erinteg_2(isupb,bl,first,nbls,acc,ikbl,npkl, where) 1336c 1337c character variable WHERE is needed for derivatives only 1338c 1339 implicit real*8 (a-h,o-z) 1340 logical first 1341 character*8 where 1342c 1343 common /primij/ iabprim, ijdim ,ijpar1 1344 common /primkl/ kabprim, kldim ,klpar1 1345c 1346 common /howmany/ ntotal,noijkl,nopres,nospec,nrimtot,nrimret 1347ctime 1348 common /timex/ tconv1,tconv2,ttrobs 1349 common /time1/ tpre4n,txwpq ,tassem,tamshf,tdesti 1350 common /time3/ tzeroi,tspeci 1351 common /time4/ tderiv 1352ctime 1353c 1354 common /types/itype,jtype,ktype,ltype,itype1,jtype1,ktype1,ltype1 1355 common /contr/ ngci,ngcj,ngck,ngcl,lci,lcj,lck,lcl,lc12,lc34 1356 common /gcont/ ngci1,ngcj1,ngck1,ngcl1,ngcd 1357c 1358#include "texas_lpar.fh" 1359c 1360 common/obarai/ 1361 * lni,lnj,lnk,lnl,lnij,lnkl,lnijkl,mmax, 1362 * nqi,nqj,nqk,nql,nsij,nskl, 1363 * nqij,nqij1,nsij1,nqkl,nqkl1,nskl1,ijbeg,klbeg 1364c 1365 common /memor3/ nblok1d 1366c 1367 common /memor4/ iwt0,iwt1,iwt2,ibuf,ibuf2, 1368 * ibfij1,ibfij2,ibfkl1,ibfkl2, 1369 * ibf2l1,ibf2l2,ibf2l3,ibf2l4,ibfij3,ibfkl3, 1370 * ibf3l,issss, 1371 * ix2l1,ix2l2,ix2l3,ix2l4,ix3l1,ix3l2,ix3l3,ix3l4, 1372 * ixij,iyij,izij, iwij,ivij,iuij,isij 1373c 1374 common /memor5x/ ieab,iecd 1375 common /memor5a/ iaa,ibb,icc,idd,icis,icjs,icks,icls, 1376 * ixab,ixp,ixpn,ixpp,iabnia,iapb,i1apb,ifij,icij,isab, 1377 * ixcd,ixq,ixqn,ixqq,icdnia,icpd,i1cpd,ifkl,ickl,iscd 1378 common /memor5b/ irppq, 1379 * irho,irr1,irys,irhoapb,irhocpd,iconst,ixwp,ixwq,ip1234, 1380 * idx1,idx2,indx 1381c 1382 common /memor5c/ itxab,itxcd,iabcd,ihabcd 1383 common /memor5d/ iabnix,icdnix,ixpnx,ixqnx,ihabcdx 1384c new for grad. derivatives: 1385 common /memor5dd/ iaax,ibbx,iccx 1386 common /memor5e/ igci,igcj,igck,igcl,indgc,igcoef, 1387 * icfg,jcfg,kcfg,lcfg, igcij,igckl 1388 common /memor5f/ indxp 1389c nmr der. 1390 common /memor6/ ixyab,ixycd 1391c 1392 common /memors/ nsym,ijshp,isymm 1393c 1394 common /pnl001/ ispec,ijpres2,klpres2, ijblock,klblock,iqorder 1395 common /pnl005/ isblsize,isblqrts,isblpoint 1396 common /map4uniq/ nij_uniqe, ij_uniqe_p, map_ij, 1397 * nkl_uniqe, kl_uniqe_p, map_kl 1398c 1399 dimension bl(*) 1400c----------------------------------------------------------------- 1401c get pnl-size of the super-block and the pointer to its quartets: 1402c 1403 call get_nbls_pnl(isupb,bl(isblsize),bl(isblpoint), 1404 * isbl_size,isbl_point) 1405c----------------------------------------------------------------- 1406c memory handling for a block 1407c (independent of contractions) 1408c reserves memory for trobsa and assemble 1409c 1410 call memo4a(bl, nbls, l11,l12,mem2,igmcnt) 1411c----------------------------------------------------------------- 1412 nfumax=nfu(mmax) 1413 mmax1=mmax-1 1414 nfumax1=nfumax 1415 narray=1 1416 nqijr=nqij 1417 nqklr=nqkl 1418 if(where.eq.'shif') nfumax1=nfu(mmax1) 1419 if(where.eq.'forc') then 1420 nfumax1=nfu(mmax1) 1421 narray=4 1422 nqij=nqij-1 1423 nqkl=nqkl-1 1424 if(nqij.le.0) nqij=1 1425 if(nqkl.le.0) nqkl=1 1426 endif 1427 if(where.eq.'hess') then 1428 nfumax1=nfu(mmax1-1) 1429 narray=10 1430 nqij=nqij-2 1431 nqkl=nqkl-2 1432 if(nqij.le.0) nqij=1 1433 if(nqkl.le.0) nqkl=1 1434 endif 1435c----------------------------------------------------------------- 1436 first=.true. 1437c----------------------------------------------------------------- 1438 igcoet=1 1439 if(ngcd.gt.1) then 1440 ngcij=ngci1*ngcj1 1441 ngckl=ngck1*ngcl1 1442 call getmem(ngcij*nbls,igcijx) 1443 call getmem(ngckl*nbls,igcklx) 1444 call getmem(ngcd*nbls ,igcoet) 1445 endif 1446 if(where.eq.'forc' .or. where.eq.'hess') then 1447 call getmem(nbls,iaax) 1448 call getmem(nbls,ibbx) 1449 call getmem(nbls,iccx) 1450 endif 1451c----------------------------------------------------------------- 1452c get pnl-size of the super-block and the pointer to its quartets: 1453c 1454 call get_nbls_pnl(isupb,bl(isblsize),bl(isblpoint), 1455 * isbl_size,isbl_point) 1456c 1457 call getmem(isbl_size,map_n0) 1458 call get_nbls_n0(isbl_size,isbl_point,bl(isblqrts), 1459 * bl(map_n0), nbls_n0) 1460c----------------------------------------------------------------- 1461c 1462c loop over contraction : 1463c 1464 lcij=0 1465 do 12 l1=1,lci 1466 do 12 l2=1,lcj 1467 lcij=lcij+1 1468 if(ngcd.gt.1) then 1469 call gcpairs(nbls,lcij,ngcij,bl(igcij),bl(igcijx) ) 1470 endif 1471 lckl=0 1472 nblsr=0 1473 do 34 l3=1,lck 1474 do 34 l4=1,lcl 1475 lckl=lckl+1 1476 if(ngcd.gt.1) then 1477 call gcpairs(nbls,lckl,ngckl,bl(igckl),bl(igcklx)) 1478 endif 1479c 1480c* calculate : rppq,rho, rysx,rhoapb,rhocpd : 1481c 1482ctime 1483 call txs_second(t4b) 1484c 1485c2002 call prec4neg_2(isbl_size,isbl_point,bl(isblqrts),npkl, 1486c 1487 call prec4neg_2(nbls_n0,bl(map_n0), 1488 * lcij,lckl, lc12,lc34, 1489 * bl(ieab),bl(iecd),bl(iapb),bl(icpd),bl(icij),bl(ickl), 1490 * bl(ixp),bl(ixq), 1491 * bl(map_ij),bl(map_kl),nij_uniqe,nkl_uniqe, 1492 * bl(irppq),bl(irhoapb),bl(irhocpd),bl(irys),bl(iconst), 1493 * nbls1,bl(indx) ) 1494c 1495ctime 1496 call txs_second(t4e) 1497 tpre4n=tpre4n+(t4e-t4b) 1498c 1499c 1500 nrimret=nrimret+nbls1 1501 nrimtot=nrimtot+nbls 1502 if(nbls1.eq.0) go to 34 1503c 1504c note : 1505c* from here a given block (nbls) was reduced to (nbls1) *** 1506c* for given l1,l2,l3,l4 !!! 1507c* and 1508c* rho, rppq, rhoapb, rhocpd, rys, and const 1509c have dimension nbls1 (not (nbls) ) 1510c****************************************************** 1511c* substitute zero for all needed bufors for quartets 1512c* which do not appear first time after neglect 1513c 1514 if(first) then 1515ctime 1516 call txs_second(tzer) 1517 nblsnot=nbls-nbls1 1518 if(nblsnot.ne.0) then 1519 call zeroint(bl,nbls,nbls1,nblsnot,lnij,lnkl,indx,ngcd,narray) 1520 endif 1521 call txs_second(tzero) 1522 tzeroi=tzeroi+(tzero-tzer) 1523 endif 1524c---------------------------------------------------------------- 1525c calculate products of contraction coefficients 1526c for general contracted basis set 1527c 1528 if(ngcd.GT.1) then 1529 call gcquart(nbls,nbls1, bl(indx), 1530 * ngci1,ngcj1,ngck1,ngcl1,ngcd, 1531 * bl(igcijx),ngcij, bl(igcklx),ngckl, 1532 * bl(indgc),bl(igcoef) ) 1533 endif 1534c---------------------------------------------------------------- 1535c needed for obara-saika recursives : 1536c calculate XWP, XWQ and p1234 where 1537c p1234= - [ b_exp*(A-B) + d_exp*(C-D) ]/(c_exp+d_exp) 1538c 1539 call txs_second(txwpqb) 1540c 1541 call xwpq_2(nbls1,bl(ixwp),bl(ixwq),bl(ip1234), 1542 * lc12,lc34,lcij,lckl, 1543 * bl(idx1),bl(idx2),bl(indx), 1544 * nij_uniqe,nkl_uniqe, 1545 * bl(irppq),bl(ixp),bl(ixq),bl(ixpp),bl(ixqq), 1546c98 * bl(itxab),bl(itxcd),bl(iabcd), 1547 * bl(itxab),bl(itxcd), 1548 * bl(iapb),bl(i1cpd),bl(icpd),bl(i1apb) ) 1549c---------------------------------------------------------------- 1550c needed for tracy's recursives : 1551c depending on numerical stability calculate 1552c abcd=(a+b)/(c+d) or (c+d)/(a+b) 1553c 1554 call abcd_2(nbls1,lcij,lckl,bl(indx), 1555 * bl(iabcd), bl(iapb),bl(i1cpd),bl(icpd),bl(i1apb), 1556 * bl(iconst),bl(igcoef),ngcd) 1557c---------------------------------------------------------------- 1558c convert: abnia,cdnia, xpn,xqn, habcd 1559c from pair to quartet (nbls1) quantities 1560c for trobsa 1561c 1562c ij 1563ctime 1564 call txs_second(tc1xb) 1565 txwpq=txwpq+(tc1xb-txwpqb) 1566c 1567 if(nbls1.ne.nbls .or. nblsr.ne.nbls) then 1568 call conv1x_2(nbls1,mmax1,nij_uniqe,lcij, bl(idx1),bl(indx), 1569 * bl(ixpn), bl(ixpnx) ) 1570 iabnix=iabnia+(lcij-1)*mmax1 1571 if(where.eq.'forc' .or. where.eq.'hess') then 1572 call conv1der_2(nbls1,nij_uniqe,l1,bl(idx1),bl(indx), 1573 * bl(iaa),bl(iaax) ) 1574 call conv1der_2(nbls1,nij_uniqe,l2,bl(idx1),bl(indx), 1575 * bl(ibb),bl(ibbx) ) 1576 endif 1577 endif 1578cxxxxkl 1579 call conv1x_2(nbls1,mmax1,nkl_uniqe,lckl, bl(idx2),bl(indx), 1580 * bl(ixqn), bl(ixqnx) ) 1581 icdnix=icdnia+(lckl-1)*mmax1 1582 if(where.eq.'forc' .or. where.eq.'hess') then 1583 call conv1der_2(nbls1,nkl_uniqe,l3,bl(idx2),bl(indx), 1584 * bl(icc),bl(iccx) ) 1585 endif 1586ctime 1587 call txs_second(tc1xe) 1588 tconv1=tconv1+(tc1xe-tc1xb) 1589cxxxx2 1590c98... if(nsij.ge.nskl) then 1591ccccc call conv2x(nbls1,nfumax1,npklx,lckl, bl(idx2),bl(indx), 1592ccccc* bl(ihabcd),nfumax, bl(ihabcdx) ) 1593 ihabcdx=ihabcd+(lckl-1)*nfumax*3 1594c98... else 1595cccc if(nbls1.ne.nbls .or. nblsr.ne.nbls) then 1596cccc call conv2x(nbls1,nfumax1,npij ,lcij, bl(idx1),bl(indx), 1597cccc * bl(ihabcd),nfumax, bl(ihabcdx) ) 1598c98 ihabcdx=ihabcd+(lcij-1)*nfumax*3 1599cccc endif 1600c98... endif 1601ctime 1602 call txs_second(tc2xe) 1603 tconv2=tconv2+(tc2xe-tc1xe) 1604c 1605c98 call trobsa(bl,nbls1,l11,l12,mem2,immax,kmmax,lobsa) 1606 call trobsa_2(bl,nbls1,l11,l12,mem2) 1607c 1608ctime 1609 call txs_second(ttrob) 1610 ttrobs=ttrobs+(ttrob-tc2xe) 1611cxxxxxx 1612 if(ngcd.eq.1) then 1613 call assemblx(bl,first,nbls,nbls1,lnij,lnkl, 1614 * l1,l2,l3,l4,lcij,lckl,nij_uniqe,nkl_uniqe) 1615cnew * l1,l2,l3,l4,lcij,lckl,npij,npklx) 1616 else 1617c98 call gcquart(nbls,nbls1, bl(indx), 1618c98 * ngci1,ngcj1,ngck1,ngcl1,ngcd, 1619c98 * bl(igcijx),ngcij, bl(igcklx),ngckl, 1620c98 * bl(indgc),bl(igcoef) ) 1621c 1622c transpose gcoef(ngcd,nbls) into gcoet(nbls,ngcd) 1623 call trspmo(bl(igcoef),ngcd,bl(igcoet),nbls) 1624c 1625 call assemblg(bl,first,nbls,nbls1,lnij,lnkl,ngcd,igcoet) 1626c ,npkl) .... spurious argument? 1627 endif 1628ctime 1629 call txs_second(tasse) 1630 tassem=tassem+(tasse-ttrob) 1631cxxxxx 1632c 1633 nblsr=nbls1 1634c 1635 34 continue 1636 12 continue 1637c 1638c grad derivatives: 1639c 1640 nqij=nqijr 1641 nqkl=nqklr 1642c 1643c--------------------------------------------- 1644c release memory : 1645c 1646 call retmem(1) ! allocation for map_n0 1647c 1648 if(where.eq.'forc' .or. where.eq.'hess') call retmem(3) 1649c 1650c transpose but2 into buf2 1651 if(ngcd.gt.1) call retmem(3) 1652c 1653 if( first ) then 1654 call retmem(igmcnt) 1655 return 1656 endif 1657c 1658c release memory from trobsa 1659c 1660 call retmem(3) 1661c 1662c--------------------------------------------- 1663c for NMR , GRADIENT and HESSIAN derivatives : 1664c 1665 lnijr=lnij 1666 lnklr=lnkl 1667c--------------------------------------------- 1668c reserves memory for amshift 1669c 1670 call memo4b(bl, nbls,igmcnt1) 1671c--------------------------------------------- 1672c for NMR derivatives : 1673c 1674c 1675 if(where.eq.'shif') then 1676 call txs_second(tderb) 1677c 1678c return to the original values of nsij,nskl and mmax : 1679c 1680 ijderiv=1 1681 klderiv=1 1682 call iobarb(ijderiv,klderiv) 1683c 1684c construct nmr/giao derivatives : (i+j,s|k+l,s)(x,y,z) 1685c 1686 call shift_der(nbls,lnijr,lnklr,nij_uniqe,nkl_uniqe,ngcd, 1687 * idx1,idx2, ixab,ixcd,ixyab,ixycd) 1688c 1689 call txs_second(tdere) 1690 tderiv=tderiv+tdere-tderb 1691 endif 1692c------------------------------------ 1693c for GRADIENT derivatives : 1694c 1695 if(where.eq.'forc') then 1696 call txs_second(tderb) 1697c 1698c return to the original values of nsij,nskl and mmax : 1699c 1700 ijderiv=1 1701 klderiv=1 1702 call iobarb(ijderiv,klderiv) 1703c 1704c construct gradient derivatives : (i+j,s|k+l,s)(x,y,z) 1705c 1706 call force_der(bl,nbls,lnijr,lnklr,nij_uniqe,ngcd,idx1,ixab) 1707c 1708 call txs_second(tdere) 1709 tderiv=tderiv+tdere-tderb 1710 endif 1711c------------------------------------ 1712c for HESSIAN derivatives : 1713c 1714 if(where.eq.'hess') then 1715 call txs_second(tderb) 1716c 1717c return to the original values of nsij,nskl and mmax : 1718c 1719 ijderiv=2 1720 klderiv=2 1721 call iobarb(ijderiv,klderiv) 1722c 1723c construct hessian derivatives : (i+j,s|k+l,s)(x,y,z) 1724c 1725 call hessian_der(bl,nbls,lnijr,lnklr,nij_uniqe,ngcd,idx1,ixab) 1726c 1727 call txs_second(tdere) 1728 tderiv=tderiv+tdere-tderb 1729 endif 1730c------------------------------------ 1731c 1732c* shift angular momentum (a->b, c->d) : 1733ctime 1734 call txs_second(tamsb) 1735cnew call amshift(bl,nbls,lnij,lnkl,npij,npklx,ngcd) 1736 call amshift(bl,nbls,lnij,lnkl,nij_uniqe,nkl_uniqe,ngcd) 1737 call txs_second(tamse) 1738 tamshf=tamshf+(tamse-tamsb) 1739c----- 1740 call retmem(igmcnt+igmcnt1-3) 1741c----- 1742c 1743 end 1744c=================================================================== 1745 subroutine erintsp_2(isupb,bl,first,nbls,acc,ikbl,npkl) 1746 implicit real*8 (a-h,o-z) 1747 logical first 1748c 1749 common /primij/ iabprim, ijdim ,ijpar1 1750 common /primkl/ kabprim, kldim ,klpar1 1751 common /howmany/ ntotal,noijkl,nopres,nospec,nrimtot,nrimret 1752 common /types/itype,jtype,ktype,ltype,itype1,jtype1,ktype1,ltype1 1753 common /contr/ ngci,ngcj,ngck,ngcl,lci,lcj,lck,lcl,lc12,lc34 1754 common /gcont/ ngci1,ngcj1,ngck1,ngcl1,ngcd 1755#include "texas_lpar.fh" 1756 COMMON/SHELL/LSHELLT,LSHELIJ,LSHELKL,LHELP,LCAS2(4),LCAS3(4) 1757 common/obarai/ 1758 * lni,lnj,lnk,lnl,lnij,lnkl,lnijkl,mmax, 1759 * nqi,nqj,nqk,nql,nsij,nskl, 1760 * nqij,nqij1,nsij1,nqkl,nqkl1,nskl1,ijbeg,klbeg 1761 common /memor3/ nblok1d 1762 common /memor4/ iwt0,iwt1,iwt2,ibuf,ibuf2, 1763 * ibfij1,ibfij2,ibfkl1,ibfkl2, 1764 * ibf2l1,ibf2l2,ibf2l3,ibf2l4,ibfij3,ibfkl3, 1765 * ibf3l,issss, 1766 * ix2l1,ix2l2,ix2l3,ix2l4,ix3l1,ix3l2,ix3l3,ix3l4, 1767 * ixij,iyij,izij, iwij,ivij,iuij,isij 1768 common /memor5x/ ieab,iecd 1769 common /memor5a/ iaa,ibb,icc,idd,icis,icjs,icks,icls, 1770 * ixab,ixp,ixpn,ixpp,iabnia,iapb,i1apb,ifij,icij,isab, 1771 * ixcd,ixq,ixqn,ixqq,icdnia,icpd,i1cpd,ifkl,ickl,iscd 1772 common /memor5b/ irppq, 1773 * irho,irr1,irys,irhoapb,irhocpd,iconst,ixwp,ixwq,ip1234, 1774 * idx1,idx2,indx 1775c 1776 common /memor5c/ itxab,itxcd,iabcd,ihabcd 1777 common /memor5d/ iabnix,icdnix,ixpnx,ixqnx,ihabcdx 1778 common /memor5e/ igci,igcj,igck,igcl,indgc,igcoef, 1779 * icfg,jcfg,kcfg,lcfg, igcij,igckl 1780 common /memor5f/ indxp 1781 common /memors/ nsym,ijshp,isymm 1782c 1783 common /pnl001/ ispec,ijpres2,klpres2, ijblock,klblock,iqorder 1784 common /pnl005/ isblsize,isblqrts,isblpoint 1785 common /map4uniq/ nij_uniqe, ij_uniqe_p, map_ij, 1786 * nkl_uniqe, kl_uniqe_p, map_kl 1787c 1788 dimension bl(*) 1789c----------------------------------------------------------------- 1790c get pnl-size of the super-block and the pointer to its quartets: 1791c 1792 call get_nbls_pnl(isupb,bl(isblsize),bl(isblpoint), 1793 * isbl_size,isbl_point) 1794c----------------------------------------------------------------- 1795c memory handling for a block 1796c (independent of contractions) 1797c reserves memory for trobsa and assemble 1798c 1799 call memo4a(bl, nbls, l11,l12,mem2,igmcnt) 1800c----------------------------------------------------------------- 1801 first=.true. 1802c----------------------------------------------------------------- 1803 igcoet=1 1804 if(ngcd.gt.1) then 1805 ngcij=ngci1*ngcj1 1806 ngckl=ngck1*ngcl1 1807 call getmem(ngcij*nbls,igcijx) 1808 call getmem(ngckl*nbls,igcklx) 1809 call getmem(ngcd*nbls ,igcoet) 1810 endif 1811c----------------------------------------------------------------- 1812c get pnl-size of the super-block and the pointer to its quartets: 1813c 1814 call get_nbls_pnl(isupb,bl(isblsize),bl(isblpoint), 1815 * isbl_size,isbl_point) 1816c 1817 call getmem(isbl_size,map_n0) 1818 call get_nbls_n0(isbl_size,isbl_point,bl(isblqrts), 1819 * bl(map_n0), nbls_n0) 1820c----------------------------------------------------------------- 1821c 1822c* loop over contraction : 1823c 1824 lcij=0 1825 do 12 l1=1,lci 1826 do 12 l2=1,lcj 1827 lcij=lcij+1 1828 if(ngcd.gt.1) then 1829 call gcpairs(nbls,lcij,ngcij,bl(igcij),bl(igcijx) ) 1830 endif 1831 lckl=0 1832 do 34 l3=1,lck 1833 do 34 l4=1,lcl 1834 lckl=lckl+1 1835 if(ngcd.gt.1) then 1836 call gcpairs(nbls,lckl,ngckl,bl(igckl),bl(igcklx)) 1837 endif 1838c 1839c* calculate : rppq,rho, rysx,rhoapb,rhocpd : 1840c 1841c2002 call precspec_2(isbl_size,isbl_point,bl(isblqrts),npkl, 1842c 1843 call precspec_2(nbls_n0,bl(map_n0), 1844 * lcij,lckl, lc12,lc34,bl(idx1),bl(idx2), 1845 * bl(ieab),bl(iecd), 1846 * bl(iapb),bl(icpd),bl(icij),bl(ickl), 1847 * bl(ixp),bl(ixq),bl(i1apb),bl(i1cpd),bl(itxab),bl(itxcd), 1848 * bl(map_ij),bl(map_kl),nij_uniqe,nkl_uniqe, 1849 * bl(irys),bl(iconst),bl(ixwp),bl(ixwq), nbls1,bl(indx) ) 1850c 1851 nrimret=nrimret+nbls1 1852 nrimtot=nrimtot+nbls 1853 if(nbls1.eq.0) go to 34 1854c 1855c note : 1856c* from here a given block (nbls) was reduced to (nbls1) *** 1857c* for given l1,l2,l3,l4 !!! 1858c* and 1859c* rho, rppq, rhoapb, rhocpd, rys, and const 1860c have dimension nbls1 (not (nbls) ) 1861c****************************************************** 1862c* substitute zero for all needed bufors for quartets 1863c* which do not appear first time after neglect 1864c 1865 if(first) then 1866 nblsnot=nbls-nbls1 1867 if(nblsnot.ne.0) then 1868ckw call zeroint(bl,nbls,nbls1,nblsnot,lnij,lnkl,indx,ngcd) 1869 call zeroint(bl,nbls,nbls1,nblsnot,lnij,lnkl,indx,ngcd,1) 1870 endif 1871 endif 1872c 1873c******** 1874c 1875c* special code for (ss ss) and (xs ss), (sx ss), (ss xs), (ss sx) 1876c* itegrals ( x= p or l ) 1877c 1878 if(ngcd.eq.1) then 1879 call specase_2(bl,first,nbls,nbls1, bl(indx), 1880 * nij_uniqe,nkl_uniqe,l1,l2,l3,l4, 1881 * bl(icis),bl(icjs),bl(icks),bl(icls), 1882 * bl(ibuf),bl(ibuf2), 1883 * bl(iconst),bl(irys),bl(ixwp),bl(ixwq),bl(irr1) ) 1884 else 1885 call gcquart(nbls,nbls1, bl(indx), 1886 * ngci1,ngcj1,ngck1,ngcl1,ngcd, 1887 * bl(igcijx),ngcij, bl(igcklx),ngckl, 1888 * bl(indgc),bl(igcoef) ) 1889c 1890c transpose gcoef(ngcd,nbls) into gcoet(nbls,ngcd) 1891 call trspmo(bl(igcoef),ngcd,bl(igcoet),nbls) 1892c 1893 call specasg(bl,first,nbls,nbls1, bl(indx),bl(idx1),bl(idx2), 1894 * bl(ibuf),bl(ibuf2), 1895 * bl(iconst),bl(irys),bl(ixwp),bl(ixwq), 1896 * ngcd,bl(indgc),bl(igcoet),lnijkl) 1897 endif 1898c 1899 34 continue 1900 12 continue 1901c 1902c release memory : 1903c 1904 call retmem(1) ! allocated for map_n0 1905 if(ngcd.gt.1) call retmem(3) 1906c---- 1907 call retmem(igmcnt) 1908c 1909 end 1910c=================================================================== 1911c instead of the dgetmo routine of blas : 1912c transpse A(la,lb) into B(lb,la) 1913c 1914 subroutine trspmo(amat ,lda, bmat ,ldb) 1915 implicit real*8 (a-h,o-z) 1916 dimension amat(lda,ldb), bmat(ldb,lda) 1917c 1918 do 10 i=1,lda 1919 do 10 j=1,ldb 1920 bmat(j,i)=amat(i,j) 1921 10 continue 1922c 1923 end 1924c=============================================================== 1925 subroutine block1_pnl(isupb, isbl_size,isbl_q, ipoint,maxsize, 1926 * ijpres2,klpres2,ijklpres,nblok1,nqrt) 1927 dimension ijpres2(*),klpres2(*),isbl_q(*) ! dim=nquart 1928 dimension ijklpres(*) 1929 dimension nblok1(2,*), nqrt(*) 1930c------------------------------------------------ 1931 ijkl=0 1932 do 100 iqp=1,isbl_size 1933 ijkl=ijkl+1 1934 ijklpres(iqp)=0 ! needed if sup-block was split 1935 iq=isbl_q(ipoint+iqp) 1936 if(iq.eq.0) go to 100 1937 ijcsq=ijpres2(iq) 1938 klcsq=klpres2(iq) 1939c 1940 nblok1(1,ijkl)=ijcsq 1941 nblok1(2,ijkl)=klcsq 1942c 1943 100 continue 1944c 1945 nqrt(1)=min(isbl_size,maxsize) 1946c 1947 end 1948c=============================================================== 1949 subroutine doquarts(isupb,isbl_size,isbl_q,ipoint, 1950 * ijklpres,nblspec) 1951 implicit real*8 (a-h,o-z) 1952 dimension isbl_q(*) ! dim=nquart 1953 dimension ijklpres(*) ! dim.=nbls 1954c------------------------------------------------ 1955c For PNL ijklpres(ijkl) is set up to zeros at the very begining 1956c------------------------------------------------ 1957 ijklspec=0 1958!DIR$ NEXTSCALAR 1959 do iqp=1,isbl_size 1960 iq=isbl_q(ipoint+iqp) 1961 if(iq.ne.0) then 1962 ijklspec=ijklspec+1 1963 ijkl=iqp 1964check it in: 1965 ijklpres(ijkl)=iq 1966 endif 1967 enddo 1968c 1969 nblspec=ijklspec 1970c 1971 end 1972c=============================================================== 1973c23456789.123456789.123456789.123456789.123456789.123456789.123456789.12 1974 subroutine indexp_pnl(isupb,isbl_size,isbl_q,ipoint,map_ij,map_kl, 1975 * indxij,indxkl,ipres,indxp,nblsp) 1976 dimension isbl_q(*) ! dim=nquart 1977 dimension map_ij(isbl_size), map_kl(isbl_size) 1978 dimension indxij(*),indxkl(*),ipres(*),indxp(*) 1979c----------------------------------------------------------------- 1980c setup relation between PRESENT c.s.quartets and pairs 1981c ijklp and ijpar=1,..., and klpar=1,... 1982c----------------------------------------------------------------- 1983 ijkl=0 1984 ijklp=0 1985 do 100 iqp=1,isbl_size 1986 ijkl=ijkl+1 1987 iq=isbl_q(ipoint+iqp) 1988 if(iq.eq.0) go to 100 1989c 1990c----------old----------------- 1991c ijcsq=ijpres2(iq) 1992c klcsq=klpres2(iq) 1993c ijpar=map_ij(ijcsq) 1994c klpar=map_kl(klcsq) 1995c----------old----------------- 1996cnew : map to uniqe 1997c 1998 ijpar=map_ij(iqp) 1999 klpar=map_kl(iqp) 2000c 2001c-----? ijkl=ijkl+1 2002c 2003 if(ipres(ijkl).ne.0) then 2004 ijklp=ijklp+1 2005 indxij(ijklp)=ijpar 2006 indxkl(ijklp)=klpar 2007 indxp(ijklp) =ijkl 2008 endif 2009c 2010 100 continue 2011c 2012c 2013 nblsp=ijklp 2014c 2015c write(82,*)' from indexp_pnl : indxIJ pairs :' 2016c write(82,88) (indxij(ii),ii=1,nblsp) 2017c write(82,*)' from indexp_pnl : indxKL pairs :' 2018c write(82,88) (indxkl(ii),ii=1,nblsp) 2019c 88 format(15(i3,1x)) 2020c 2021 end 2022c==================================================================== 2023 subroutine get_nbls_pnl(isupb,isbl_s,isbl_p, isbl_size,isbl_poin) 2024 dimension isbl_s(*),isbl_p(*) 2025c 2026c returns size of the pnl-block - isbl_size 2027c and pointer to its quartets - iq=isbl_q(IPOINT+iqp) 2028c 2029 isbl_size=isbl_s(isupb) 2030 isbl_poin=isbl_p(isupb) 2031c 2032 end 2033c==================================================================== 2034 subroutine doblock2(ibl,ijblock,idoit) 2035 dimension ijblock(*) 2036c 2037 idoit=ijblock(ibl) 2038c 2039 end 2040c==================================================================== 2041 subroutine dosupblk(isupb,isbl_s,isbl_q,ipoint,idoit,isbl_split) 2042 common /pnl006/ nsplit,isplit,isblsplit, isblpart 2043 dimension isbl_s(*),isbl_q(*),isbl_split(0:*) 2044ctest 2045c write(6,*)'FROM DOSUPBLK :isplit=',isplit,' rep.no=',nsplit 2046c write(6,*)'superblock=',isupb,' its size=',isbl_s(isupb) 2047ctest 2048c 2049 isbl_size=isbl_s(isupb) 2050 idoit=isbl_size 2051 if(idoit.eq.0) return 2052c 2053 icount=0 2054 do 10 iqp=1,isbl_size 2055 iq=isbl_q(ipoint+iqp) 2056 if(iq.gt.0) icount=icount+1 2057 10 continue 2058c 2059 if(icount.eq.0) then 2060 idoit=0 2061 return 2062 endif 2063c 2064 if(isplit.gt.0) then 2065 ibeg=isbl_split(nsplit-1) 2066 iend=isbl_split(nsplit)-1 2067ccc write(6,*)'superblock=',isupb,' in range=',ibeg,iend 2068 if(ibeg.le.isupb .and. isupb.le.iend) then 2069 else 2070 idoit=0 2071 endif 2072 endif 2073c 2074 end 2075c=============================================================== 2076 subroutine copy_block(isbl_size,isbl_q,isbl_point,isbl_copy ) 2077 dimension isbl_q(*),isbl_copy(*) 2078c 2079 do 10 iqp=1,isbl_size 2080 iq=isbl_q(isbl_point+iqp) 2081 isbl_copy(iqp)=iq 2082 10 continue 2083c 2084 end 2085c=============================================================== 2086 subroutine check_in (isupb,isbl_size,isbl_q,isbl_point,isbl_copy, 2087 * maxsize) 2088 implicit real*8 (a-h,o-z) 2089c------------------------------------------------------------ 2090 common /time5/ tcheck,tuniq2,tmap4u 2091c------------------------------------------------------------ 2092 dimension isbl_q(*),isbl_copy(*) 2093c----- 2094 call txs_second(time1) 2095c----- 2096c 2097c copy back first MAXSIZE non-zero quartets from isbl_copy to isbl_q : 2098c 2099 nonzero=0 2100 do 10 iqp=1,isbl_size 2101 isbl_q(isbl_point+iqp)=0 2102c 2103 iq=isbl_copy(iqp) 2104 if(iq.gt.0) then 2105 nonzero=nonzero+1 2106 if(nonzero.le.maxsize) then 2107 isbl_q(isbl_point+iqp)=iq 2108 endif 2109 endif 2110 10 continue 2111c 2112c----- 2113 call txs_second(time2) 2114 tcheck=tcheck+time2-time1 2115c----- 2116c 2117 end 2118c=============================================================== 2119 subroutine check_out(isupb,isbl_s,isbl_q,isbl_point,isbl_split, 2120 * isbl_copy,nparts,maxsize, 2121 * ibl,ijblock,igo_back ) 2122 implicit real*8 (a-h,o-z) 2123c------------------------------------------------------------ 2124 common /time5/ tcheck,tuniq2,tmap4u 2125c------------------------------------------------------------ 2126 common /pnl006/ nsplit,isplit,isblsplit, isblpart 2127 dimension isbl_s(*),isbl_q(*), isbl_copy(*) 2128 dimension isbl_split(0:*) 2129 dimension ijblock(*) 2130c-------------------------------------------------------- 2131 call txs_second(time1) 2132c-------------------------------------------------------- 2133 igo_back=0 2134c 2135 isbl_size=isbl_s(isupb) 2136c 2137 IF(nparts.eq.1) THEN 2138 do 10 iqp=1,isbl_size 2139 isbl_q(isbl_point+iqp)=0 2140 10 continue 2141c 2142 if(isplit.gt.0) then 2143 ibeg=isbl_split(nsplit-1) 2144 iend=isbl_split(nsplit)-1 2145 if(ibeg.le.isupb .and. isupb.le.iend) then 2146 if(isupb.gt.1) ijblock(ibl-1)=0 2147 if(isupb.eq.iend) igo_back=1 2148 endif 2149 endif 2150 ENDIF 2151c 2152 IF(nparts.gt.1) THEN 2153c zero out first MAXSIZE non-zero quartets in isbl_copy : 2154 nonzero=0 2155 do 20 iqp=1,isbl_size 2156 iq=isbl_copy(iqp) 2157 if(iq.gt.0) then 2158 nonzero=nonzero+1 2159 if(nonzero.le.maxsize) then 2160 isbl_copy(iqp)=0 2161 endif 2162 endif 2163 20 continue 2164c 2165c copy from isbl_copy into isbl_q : 2166c 2167 non0_left=0 2168 do 30 iqp=1,isbl_size 2169 iq=isbl_copy(iqp) 2170 if(iq.gt.0) non0_left=non0_left+1 2171 isbl_q(isbl_point+iqp)=iq 2172 30 continue 2173c 2174 if(isplit.eq.0) write(6,*)' NPATRS=',nparts,' ISPLIT=0 !!' 2175c 2176 ibeg=isbl_split(nsplit-1) 2177 iend=isbl_split(nsplit)-1 2178 if(isplit.gt.0) then 2179 if(ibeg.le.isupb .and. isupb.le.iend) then 2180ckwol ? if(isupb.gt.1) ijblock(ibl-1)=0 2181 igo_back=1 2182ckwol ? if(isupb.eq.iend) igo_back=1 2183 endif 2184 endif 2185c 2186 if(non0_left.gt.0) then 2187 nsplit=nsplit-1 2188 else 2189ctry if(isupb.lt.iend) nsplit=nsplit-1 2190 endif 2191c 2192 ENDIF 2193c 2194c-------------------------------------------------------- 2195 call txs_second(time2) 2196 tcheck=tcheck+time2-time1 2197c-------------------------------------------------------- 2198 end 2199c=============================================================== 2200 subroutine get_nparts(isupb,isbl_part,nparts) 2201 dimension isbl_part(*) 2202c 2203 nparts=isbl_part(isupb) 2204c 2205 end 2206c======================================================================= 2207 subroutine make_map2uniq(bl,ibl,kbl,isbl_size,isbl_q,isbl_point, 2208 * ijpres2,klpres2) 2209 implicit real*8 (a-h,o-z) 2210#include "errquit.fh" 2211c------------------------------------------------------------ 2212 common /time5/ tcheck,tuniq2,tmap4u 2213c------------------------------------------------------------ 2214c output : 2215c 2216 common /map4uniq/ nij_uniqe, ij_uniqe_p, map_ij, 2217 * nkl_uniqe, kl_uniqe_p, map_kl 2218c------------------------------------------------------------ 2219 dimension bl(*) 2220 dimension isbl_q(*), ijpres2(*),klpres2(*) 2221c-------------------------------------------------------- 2222 call txs_second(time1) 2223c-------------------------------------------------------- 2224c 2225c number of uniqe pairs in IBL & Kbl (nij_uniqe,nkl_uniqe) and pointers 2226c to the ijcs shell-pairs are obtained in prec2ij & prec2kl 2227c 2228 call memo1_int(isbl_size,map_ij) 2229 call memo1_int(isbl_size,map_kl) 2230c 2231 call map_2_uniq(isbl_size,isbl_q,isbl_point,ijpres2,klpres2, 2232 * nij_uniqe, bl(ij_uniqe_p), 2233 * nkl_uniqe, bl(kl_uniqe_p), 2234 * bl(map_ij),bl(map_kl) ) 2235c-------------------------------------------------------- 2236 call txs_second(time2) 2237 tmap4u=tmap4u+time2-time1 2238c-------------------------------------------------------- 2239c 2240 end 2241c=============================================================== 2242 subroutine map_2_uniq(isbl_size,isbl_q,isbl_point, 2243 * ijpres2,klpres2, 2244 * nij_uniqe, ij_uniqe, 2245 * nkl_uniqe, kl_uniqe, 2246 * map_ij,map_kl) 2247 dimension isbl_q(*), ijpres2(*),klpres2(*) 2248 dimension ij_uniqe(nij_uniqe),kl_uniqe(nkl_uniqe) 2249 dimension map_ij(isbl_size),map_kl(isbl_size) 2250c 2251 do 10 iqp=1,isbl_size 2252 iq=isbl_q(isbl_point+iqp) 2253 if(iq.eq.0) go to 10 2254 ijcsq=ijpres2(iq) 2255 klcsq=klpres2(iq) 2256c 2257 call seartch_4_uniqe(ijcsq,ij_uniqe,nij_uniqe,ij_u) 2258 call seartch_4_uniqe(klcsq,kl_uniqe,nkl_uniqe,kl_u) 2259c 2260 map_ij(iqp)=ij_u 2261 map_kl(iqp)=kl_u 2262 10 continue 2263c 2264 end 2265c=============================================================== 2266 subroutine seartch_4_uniqe(ijcsq,ij_uniqe,nij_uniqe,ij_u) 2267 dimension ij_uniqe(nij_uniqe) 2268c 2269 if(nij_uniqe.eq.0) then 2270 write(6,*)' number of uniqe pair is ZERO ; wrong' 2271 call errquit('s_4_u:execution stopped in seartch_4_uniqe',0, 2272 & INT_ERR) 2273 endif 2274 if(nij_uniqe.eq.1) then 2275 ij_u=1 2276 return 2277 endif 2278c 2279 ijhalf=nij_uniqe/2 2280 ijdelta=ijhalf/2 2281c 2282 100 continue 2283 ijcsh =ij_uniqe(ijhalf) 2284c write(6,*)' ijhalf=',ijhalf,' ijcsh=',ijcsh 2285 if(ijcsq.eq.ijcsh) then 2286 ij_u=ijhalf 2287 return 2288 endif 2289c 2290 if(ijcsq.gt.ijcsh) then 2291 ijhalf=ijhalf + ijdelta 2292 if(ijhalf.gt.nij_uniqe) ijhalf=nij_uniqe 2293 ijdelta=ijdelta/2 2294 if(ijdelta.eq.0) ijdelta=1 2295 go to 100 2296 endif 2297 if(ijcsq.lt.ijcsh) then 2298 ijhalf=ijhalf - ijdelta 2299 if(ijhalf.le. 0 ) ijhalf=1 2300 ijdelta=ijdelta/2 2301 if(ijdelta.eq.0) ijdelta=1 2302 go to 100 2303 endif 2304c 2305 end 2306c======================================================================= 2307 subroutine get_nbls_n0(isbl_size,isbl_point,isbl_q,map_n0,nbls_n0) 2308 dimension isbl_q(*) 2309 dimension map_n0(*) 2310c 2311 ijklp=0 2312 do 100 iqp=1,isbl_size 2313 iq=isbl_q(isbl_point+iqp) 2314 if(iq.eq.0) go to 100 2315 ijklp=ijklp+1 2316 map_n0(ijklp)=iqp 2317 100 continue 2318 nbls_n0=ijklp 2319c 2320 end 2321c======================================================================= 2322