1! 2! Dalton, a molecular electronic structure program 3! Copyright (C) by the authors of Dalton. 4! 5! This program is free software; you can redistribute it and/or 6! modify it under the terms of the GNU Lesser General Public 7! License version 2.1 as published by the Free Software Foundation. 8! 9! This program is distributed in the hope that it will be useful, 10! but WITHOUT ANY WARRANTY; without even the implied warranty of 11! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 12! Lesser General Public License for more details. 13! 14! If a copy of the GNU LGPL v2.1 was not distributed with this 15! code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html. 16! 17! 18C 19#ifdef OLD_REV_OLG 20!=========================================================================== 21! 22!--- 23!971201-hjaaj: TRACTL: make check for transformation already there 24! also work when AOSRTINT.DA has been deleted !!! 25!951219-hjaaj: TRACTL: defined DM1=-1.0D0 (was undefined); never 26! modify input KTRLVL (was ITRLVL, which is abs(KTRLVL)) 27!941115-hjaaj: SIR_INTOPEN: included FLAG(34) in def. of OLDTRA 28!940927-hjaaj: SIR_INTOPEN: s/FLAG(9)/DOMP2/ and s/FLAG(8)/DORSP/ 29!940909-hjaaj: added PARAGON define label (same code as AIX) 30!940517-hjaaj: Ide: gem integraler paa LUINTM uden labels (afhaengigt af ITRLVL) 31! Brug ITRLVL til at bestemme LBINTM (ikke NMORBT men f.eks. NMOCCT + ...) 32! Beregn ogsaa LBINTD v.hj.af LBMXSQ ?? 33!MAERKE 930305: usikkert,indfoer USEDRC el. LVLDRC i parameterliste til TRACTL 34!-------- Revision n05 ends here --------- 35!940517-hjaaj 36!SIR_INTOPEN: deleted write DUMMY for CRAY; attempt open old LUINTM file if FLAG(14) 37! true, irrespective of FLAG(34); fixed ICASE print error on Cray 38! by extending IF (OPINTM) THEN to after ICASE print 39!-------- Revision n04 ends here --------- 40!931201-hjaaj 41!TRACTL: corrected error for TRAAB; LRBUF for 600 but LBINTM could be up to 1024 42! new parameter LBMXSQ specifying max buffer length for seq. files, used for 43! calculating LBINTM and as buffer length for LUHALF 44!930720-hjaaj 45!TRACTL: new parameter PROCC false in CALL PRORB 46!921203-hjaaj revision n04 (based on revision m02) 47!- added define label DEC (same selection as AIX); removed define label FPS 48!SORTA : LDAMIN = LDAMAX / 4 (was LDAMAX / 2) 49!TRACTL: print KWORK,MWORK when IPRTRA .gt. 5 (was .ge. 20); 50! SAVE LRBUF; LRBUF defined for 600 integrals/buffer in AO file 51!910319-hjaaj 52!SORTA,SORTB: reduced NBLOCK to max NMBASX, thus only calling DZERO with 53! what is used later; at the same time this gives the most even 54! distribution of distributions per chain. 55!TRACD,TRAAB: insert symmetry check in write-to-file loop, gprof showed 56! for Neon that a lot of time was used in skipping zero integrals. 57!901025-hjaaj 58! TRACTL: JTER=1 only if IPRTRA.gt.1, print TIME IN SORTA if iprtra.gt.0 59!900801-hjaaj 60! Moved writing of pre-MOLTWOEL information on LUINTM from TRAAB to TRACTL; 61! added NSYM,NORB,NBAS and CMO to information. 62! TRACTL: check ITRLVL,CMO against LUINTM information, abandon transformation 63! if the needed mo integrals already available. 64! SIR_INTOPEN: revised ITRLVL check against LUINTM for levels 2 and 3; 65! rationalized file openings for different computer types; 66! write warning if flag(14) and integrals not found. 67!900307-hjaaj 68! TRACTL: empty LUINTM in beginning of transformation to give more disk 69! space for temporary files during the transformation. 70!900305-hjaaj 71! TRACTL: only call drcctl if ITRLVL.gt.0 (not ci) and usedrc 72! MAERKE: usikkert,indfoer USEDRC el. LVLDRC i parameterliste til TRACTL 73!900302-hjaaj 74! TRACTL: reset /CBGETD/ to tell any H2AC on disk now obsolete. 75!900226-hjaaj 76! TRACTL: if (usedrc) then call DRCCTL 77!900115-hjaaj 78! NXTH2M,NXTH2D: needtp(ityp) .lt. 0 now means at least one distribution 79! of that type needed, but not all (in that case needtp(ityp) .gt. 0) 80!900111-hjaaj 81! SIR_INTOPEN: local logical variable OPINTM, added FLAG(8:9) [RESPONSE and 82! MP2] to cases when LUINTM is opened. 83!900108-hjaaj 84! NXTH2D: read LBINTD,LVLDRC after label DRCINFO on LUINTD 85!891208-hjaaj 86! SORTA,TRACD,TRAAB: print threshold for discarding integrals 87! TRACTL: if (thrp.le.0) thrp=1.d-12 (was 1.d-9) 88!891116-hjaaj 89! NXTH2M,NXTH2D: more tests that buffers are allocated/have not been 90! corrupted. Release buffer allocation when no more buffers. 91!890710-hjaaj 92! Write label 'END INTM' at end of LUINTM in TRAAB 93! Write total time in tracd,sortb,traab in tractl if iprtra .gt.0 94! and only write TRANSFORMATION level each time if iprtra .gt. 0 95! and corrected 890615 tractl output error in TRACTL 96!=========================================================================== 97#endif 98C /* Deck trauth */ 99 SUBROUTINE TRAUTH(IWUNIT) 100C 101C 891117-hjaaj 102C Author of two-electron integral transformation routines. 103C 104 WRITE(IWUNIT,100) 105 100 FORMAT(/T2,'Two-electron integral transformation:'/ 106 * T5,'Peter Taylor, University of Lund, Sweden'/, 107 * T5,'Hans Joergen Aa. Jensen, University of Aarhus, Denmark') 108 RETURN 109 END 110C /* Deck sir_intopen */ 111 SUBROUTINE SIR_INTOPEN 112C 113C 3-Nov-1989 Hans Joergen Aa. Jensen 114C 115C Open mo integral units required by this Sirius calculation, 116C as specified in command input. 117C 118#include "implicit.h" 119#include "dummy.h" 120C 121#include "iratdef.h" 122#include "lbmxsq.h" 123C 124C Used from common blocks: 125C INFINP : FLAG(*),ITRLVL,ITRFIN,ICI0 126C INFORB : NNORBT 127C INFDIM : LBUFMA 128C INFTRA : LSRTAO, THRP 129C CCOM : THRS ! from hermit 130C INFTAP : LU* 131C 132#include "maxorb.h" 133#include "priunit.h" 134#include "infinp.h" 135#include "inforb.h" 136#include "infdim.h" 137#include "inftra.h" 138#include "maxaqn.h" 139#include "ccom.h" 140#include "inftap.h" 141C 142 LOGICAL OPINTM, FNDLAB, OLDTRA, FEXIST 143C 144C 145C Find out if we need LUINTM (MOTWOINT): 146C we need LUINTM in MC/CI optimization or for WRGEOM 147C or for MP2 or for RESPONSE 148C 149 OPINTM = .NOT.FLAG(34) .OR. (FLAG(25) .AND. ABS(ITRFIN).LE.10) 150 & .OR. DOMP2 .OR. DORSP 151 152 IF (.NOT. OPINTM) RETURN 153 154 155 IF (THRP .LT. THRS) THEN 156 WRITE(LUPRI,'(/A,1P,D10.2,A,D10.2)') 157 & ' INFO: Changed THRP threshold used in'// 158 & ' integral transformation from',THRP, 159 & ' to integral evaluation threshold',THRS 160 THRP = THRS 161 END IF 162 IF (NBAST .GT. 255 .AND. .NOT. NEWTRA) THEN 163 NEWTRA = .TRUE. 164 WRITE(LUPRI,'(/4X,A/4X,A)') 'INFORMATION: Switched to "new" '/ 165 & /'integral transformation (from 1988 ;-) )', 166 & 'because more than 255 basis functions.' 167 END IF 168 IF (FCKTRA_TYPE .gt. 0) THEN 169 CALL SIR_INTOPEN_FCKTRA 170 ! RETURN -- no, continue below to create MOTWOINT for now 171 ELSE IF (NEWTRA) THEN 172 CALL SIR_INTOPEN_NEWTRA 173 RETURN 174 END IF 175C 176C *** Specify LSRTAO = 0 to tell TRACTL to sort ao integrals 177C on first call; if old LUINTM LSRTAO is reset to -1, 178C signalling TRACTL to check if LUDA exist. 179C 180 LSRTAO = 0 181C 182C 183 IF (OPINTM) THEN 184 OLDTRA = FLAG(14) .AND. .NOT.FLAG(34) 185C ... FLAG(14) means .OLD TRANSFORMATION has been set in input 186C or, if also FLAG(34), that transformation is not needed for 187C next module (which will then be RHF). 188 NBINTM = (NNORBT-1)/LBMXSQ + 1 189 LBINTM = (NNORBT-1)/NBINTM + 1 190C 191C Always attempt to open LUINTM to see if AOSRTINT.DA information 192C is available 193C 194 ICASE = 1 195 CALL GPINQ(FNINTM,'EXIST',FEXIST) 196 IF (FEXIST) THEN 197 IF (LUINTM .LE. 0) CALL GPOPEN(LUINTM,FNINTM,'OLD',' ', 198 & 'UNFORMATTED',IDUMMY,.FALSE.) 199 REWIND (LUINTM) 200 IF (FCKTRA_TYPE .GT. 0) THEN 201 ICASE = 0 202 ELSE IF (.NOT.FNDLAB('MOLTWOEL',LUINTM)) THEN ! not for FCKTRA, MOLTWOEL only for "old" transformation 203 ICASE = 2 204 FLAG(14) = .FALSE. 205 ELSE 206 ICASE = 0 207 LSRTAO = -1 208C ... use of old AOSRTINT.DA may be possible 209 END IF 210 IF (FLAG(14)) THEN 211 KTRLVL = ABS(ITRLVL) 212 REWIND (LUINTM) 213 READ (LUINTM) 214 READ (LUINTM) LBINTM, JTRLVL 215 ! Special tests because JTRLVL .gt. KTRLVL does not always mean all required integrals available 216 IF (JTRLVL .EQ. 3 .AND. KTRLVL .EQ. 2) THEN 217 FLAG(14) = .FALSE. ! (aa|ii) and (ai|ai) missing 218 ELSE IF (JTRLVL .EQ. 2 .AND. KTRLVL .EQ. 3) THEN 219 ! OK 220 ELSE IF (JTRLVL .EQ. 6 .AND. KTRLVL .EQ. 4) THEN 221 FLAG(14) = .FALSE. ! (gg|ii) and (gi|gi) missing 222 ELSE IF (JTRLVL .EQ. 6 .AND. KTRLVL .EQ. 5) THEN 223 FLAG(14) = .FALSE. ! (gg|gi) missing 224 ELSE IF (JTRLVL .LT. KTRLVL) THEN 225 FLAG(14) = .FALSE. ! some integrals missing 226 END IF 227 IF (.NOT. FLAG(14)) ICASE = 3 228 END IF 229 ELSE 230 FLAG(14) = .FALSE. 231C no old luintm opened, then old integrals not available 232 IF (LUINTM .LE. 0) CALL GPOPEN(LUINTM,FNINTM,'UNKNOWN',' ', 233 & 'UNFORMATTED',IDUMMY,.FALSE.) 234 END IF 235 IF (OLDTRA .AND. ICASE .GT. 0) THEN 236 NINFO = NINFO + 1 237 WRITE (LUPRI,'(/A,I2,3(/A))') 238 * ' INFO SIR_INTOPEN, old integral transformation not' 239 * //' found as expected, ICASE =',ICASE, 240 * ' - ICASE = 1: MO integral file MOTWOINT does not exist', 241 * ' - ICASE = 2: no MO integrals on MOTWOINT', 242 * ' - ICASE = 3: transformation level on MOTWOINT' 243 * //' not sufficient' 244 IF (ICASE .EQ. 3) WRITE(LUPRI,'(A,2I5)') 245 * ' ITRLVL requested & ITRLVL on MOTWOINT :',ITRLVL,JTRLVL 246 END IF 247 LBUFMA = MAX(LBUFMA,NNORBT,LBINTM) 248 END IF 249C 250C END OF SIR_INTOPEN (Open MO integral files) 251C 252 RETURN 253 END 254C /* Deck sorta */ 255 SUBROUTINE SORTA(LUDA,RINT,IINTPK,SCR,ISCR,MEMSX,MEMTX,IT) 256C 257C.... THIS SUBROUTINE PERFORMS A SORT OF THE AO INTEGRAL FILE PRIOR 258C.... TO THE FIRST HALF OF THE TWO-ELECTRON TRANSFORMATION. THE 259C.... ALGORITHM INVOLVES MULTIPLE PASSES OVER THE AO INTEGRAL FILE 260C.... AND THE NUMBER OF PASSES AND THE BUFFER SIZE ON THE BACK- 261C.... CHAINED DIRECT ACCESS SORT FILE IS DETERMINED BY BOTH THE 262C.... AVAILABLE CORE DURING THE TRANSFORMATION AND BY THE REQUIRE- 263C.... MENT THAT THE TOTAL NUMBER OF I/O REQUESTS SHOULD BE A 264C.... MINIMUM. THE ROUTINE ASSUMES THAT THE ENTIRE SORT FILE 265C.... CAN BE HELD ON SYSTEM DISK. 266C 267C.... PETER R. TAYLOR MAR *79 268C 269C.... FORMAL PARAMETERS 270C 271C.... RINT -- START ADDRESS OF INTEGRAL BUFFER FOR AO FILE 272C... IINTPK-- (AS RINT - ALLOWS INTEGER*4 ADDRESSING) 273C.... SCR -- START OF SCRATCH WORK AREA 274C.... ISCR -- (AS SCR - ALLOWS INTEGER ADDRESSING) 275C.... MEMS -- LENGTH OF SCRATCH SORT AREA 276C.... MEMT -- ** * ** TRANSFORMATION AREA 277C 278#include "implicit.h" 279 DIMENSION ISCR(*),SCR(*),RINT(*) 280 INTEGER*4 IINTPK(*) 281 INTEGER IT(*) 282 INTEGER, ALLOCATABLE :: IINDX4(:,:) 283#include "iratdef.h" 284#include "dummy.h" 285C 286 PARAMETER (IBUFL = 500) 287 INTEGER IBUF(IBUFL), NAOS(8) 288C 289 LOGICAL FEXIST, OLDDX 290C 291#include "priunit.h" 292#include "inttra.h" 293#include "inftra.h" 294#include "inftap.h" 295#include "ibtdef.h" 296C 297 CALL QENTER('SORTA ') 298C 299C.... PARAMETERS USED HERE 300C 301C.... THRP: THRESHOLD TO DISCARD INTEGRALS DURING SORT 302C -- FROM /INFTRA/ , DEFINED IN INPUT -- 303C 304C.... LDAMAX: MAXIMUM DA BUFFER LENGTH (-- now set in TRACTL --) 305C 306C 307C.... MEMS is the memory available for the sort in this routine. 308C MEMT is the memory available for the integral transformation 309C to follow after the sort. 310C 311C Determine number of chains from MEMT: 312C NBLOCK is max number of simult. distrib. in transf. step 313C NCHAIN is number of chains for this blocking in the 314C transformation step 315C 316C 317 MEMS = MEMSX 318 MEMT = MEMTX 319 NBLOCK = MEMT/NMBASX 320 NCHAIN = 1 + (NMBASX-1)/NBLOCK 321 NBLOCK = 1 + (NMBASX-1)/NCHAIN 322C 323C.... FIND NSTEP, THE NUMBER OF PASSES OVER SEQUENTIAL AO FILE 324C 325C 326C here LDABUF = DA buffer length for NSTEP = 1 327C 328 LDABUF = (IRAT*MEMS) / ((1+IRAT)*NCHAIN) - 1 329 IF (LDABUF .GE. LDAMAX) THEN 330 NSTEP = 1 331 ELSE 332 LDAMIN = LDAMAX / 4 333 NSTEP = LDAMIN / LDABUF + 1 334 END IF 335C 336C.... THIS DEFINES THE NUMBER OF PASSES AS NSTEP. 337C.... NOW DETERMINE THE REMAINING SORT STATISTICS 338C 339C.... NBATCH: THE NUMBER OF BATCH BUFFERS USED IN EACH PASS 340C HERE IN SORTA 341C 342 NBATCH = 1 + (NCHAIN-1)/NSTEP 343C 344C.... CHECK STATISTICS FOR OVERFLOW OF MAXIMUM INSTALLED VALUES 345C 346 IF (NCHAIN .GT. MAXCHA) THEN 347 WRITE(LUPRI,200) NSTEP,NBLOCK,NCHAIN,NBATCH 348 WRITE(LUPRI,153) MAXCHA 349 CALL QUIT('SORTA: TOO MANY CHAINS') 350 END IF 351 153 FORMAT(/' TRACTL.SORTA: TOO MANY CHAINS, MORE THAN',I5) 352C 353#if defined (VAR_OLDCODE) 354 IF (NBATCH .GT. IBUFL) THEN 355 WRITE(LUPRI,200) NSTEP,NBLOCK,NCHAIN,NBATCH 356 WRITE(LUPRI,151) IBUFL 357 CALL QUIT('SORTA: TOO MANY BUFFERS') 358 END IF 359 151 FORMAT(/' TRACTL.SORTA: MORE THAN',I5,' BUFFERS PER PASS', 360 1 /' Increase IBUFL in SORTA (and probably also in SORTB).') 361#else 362 IF (NBATCH .GT. IBUFL) THEN 363 NSTEP = (NCHAIN-1)/IBUFL + 1 364 NBATCH = 1 + (NCHAIN-1)/NSTEP 365 END IF 366#endif 367C 368 IF (JTER.EQ.1) WRITE(LUPRI,200) NSTEP,NBLOCK,NCHAIN,NBATCH 369C 370C.... NUMBER OF PAIR INDICES PER CHAIN (EXCEPT THE LAST) 371C 372 NPQ = NBLOCK 373 NDIV = NBATCH*NBLOCK 374C 375C.... BUFFER LENGTH 376C 377 LDABUF = (MEMS/(NBATCH+NBATCH/IRAT+1)) - 1 378 LDABUF = MIN(LDABUF,LDAMAX) 379 LDABUF = (LDABUF/IRAT) * IRAT 380 IF (JTER.EQ.1) WRITE (LUPRI,201) LDABUF 381C 382 IDABUF = IRAT*LDABUF 383 LBUF = IDABUF + LDABUF + 2 384 NDAREC = (2 * LPQRS - NMBASX)/ LDABUF + 1 385C 386C... OPEN AO INTEGRAL FILE AO2INTFILE_LABEL, typically = "AOTWOINT" 387C 388 CALL GPOPEN(LUINTA,AO2INTFILE_LABEL,'OLD',' ','UNFORMATTED', 389 & IDUMMY,.FALSE.) 390 CALL MOLLAB('BASINFO ',LUINTA,LUPRI) 391 READ (LUINTA) NSYMAO, NAOS,LBUFAO,NIBUFAO,NBITSAO, LENINT4 392C 393 allocate(IINDX4(4,LBUFAO)) 394C 395C... OPEN DIRECT ACCESS FILE FOR SORTED AO INTEGRALS "AOSRTINT.DA" 396C 397 CALL GPINQ('AOSRTINT.DA','EXIST',FEXIST) 398 IF (FEXIST) THEN 399 IF (LUDA .LE. 0) CALL GPOPEN(LUDA,'AOSRTINT.DA','OLD','DIRECT', 400 & ' ',LBUF,OLDDX) 401 CALL GPCLOSE(LUDA,'DELETE') 402 END IF 403 CALL GPOPEN(LUDA,'AOSRTINT.DA','NEW','DIRECT',' ',LBUF,OLDDX) 404C 405C.... PREPARE LOOK-UP TABLE OF BUFFER STARTING ADDRESSES 406C 407 KAP = 0 408 DO 25 K = 1,NBATCH 409 IBUF(K) = KAP 410 KAP = KAP + LBUF 41125 CONTINUE 412C 413C.... AO INTEGRAL FILE SPEC; MOLECULE FORMAT 414C 415 LRINT = 2*LBUFAO ! factor 2 because IINTPK is always integer*4 416 LENINT = (2+NIBUFAO)*LBUFAO + 1 417 IF (LENINT .NE. LENINT4) 418 & CALL QUIT('LENINT .ne. LENINT4 from AOTWOINT') 419C 420C.... BEGIN THE LOOP OVER PASSES THROUGH AO INTEGRAL FILE 421C 422 ISUM = 0 423 IDISK = 1 424 NFIN = 0 425 DO 100 ISQ = 1,NSTEP 426 NST = NFIN + 1 427 NFIN = NFIN + NDIV 428 NOB = NBATCH 429C 430C.... RESET THE NUMBER OF BUFFERS(IF NECESSARY) FOR THE LAST PASS 431C 432 IF (NFIN.GT.NMBASX) THEN 433 NOB = (NMBASX-NST+NBLOCK)/NBLOCK 434 NFIN = NMBASX 435 END IF 436C 437C.... NOB NOW GIVES THE NUMBER OF BUFFERS IN USE ON THIS PASS 438C 439C.... INITIALIZE BUFFERS 440C 441 KAP = 0 442 DO K = 1,NOB 443 KAP = KAP + LBUF 444 ISCR(KAP-1) = -1 445 ISCR(KAP) = 0 446 END DO 447C 448C.... BEGIN LOOP OVER THE INTEGRAL FILE 449C 450C.... REWIND IT FIRST 451C 452 REWIND(LUINTA) 453 CALL MOLLAB('BASTWOEL',LUINTA,LUPRI) 454 1 CONTINUE 455 CALL READI4(LUINTA,LENINT,IINTPK) 456 CALL AOLAB4(IINTPK(LRINT+1),LBUFAO,NIBUFAO,NBITSAO, 457 & IINDX4,LENGTH) 458 IF(LENGTH.EQ. 0)GOTO 1 459 IF(LENGTH.EQ.-1)GOTO 2 460C 461C.... LOOP OVER INTEGRALS 462C 463 DO 3 M = 1,LENGTH 464C 465C.... PICK UP INTEGRAL, LABEL AND UNPACK THE LATTER 466C.... NOTE THAT STANDARD ORDER IS NOT CERTAIN IN SYMINT, SO 467C.... P.LT.Q ETC IS TESTED FOR. 468C 469 VALUE = RINT(M) 470C 471C.... CHECK INTEGRAL AGAINST THRESHOLD 472C 473 IF (ABS(VALUE).LT.THRP) GOTO 3 474 475 JS = IINDX4(4,M) 476 JR = IINDX4(3,M) 477 JQ = IINDX4(2,M) 478 JP = IINDX4(1,M) 479 IF (JP.GE.JQ) THEN 480 IPQ = IT(JP) + JQ 481 ELSE 482 IPQ = IT(JQ) + JP 483 END IF 484 IF (JR.GE.JS) THEN 485 IRS = IT(JR) + JS 486 ELSE 487 IRS = IT(JS) + JR 488 END IF 489C 490C.... CHECK PQ FIRST 491C 492 IF (IPQ.GE.NST .AND. IPQ.LE.NFIN) THEN 493C 494C.... PQ IS IN RANGE. DETERMINE APPROPRIATE BUFFER 495C 496 IPQS = IPQ - NST + 1 497 IPQB = (IPQS+NPQ-1)/NPQ 498C 499C.... ALLOCATE INTEGRAL AND LABEL TO BUFFER 500C 501 IBAD = IBUF(IPQB) 502 IRBAD = IBAD/IRAT 503 LENGTH = ISCR(IBAD+LBUF) + 1 504 SCR(IRBAD+LENGTH) = VALUE 505C ISCR(IBAD+IDABUF+LENGTH) = IPQ*(2**16) + IRS 506 ISCR(IBAD+IDABUF+LENGTH) = 507 * IOR( ISHFT( IPQ , 16 ) , IRS ) 508 ISCR(IBAD+LBUF) = LENGTH 509 IF (LENGTH .EQ. LDABUF) THEN 510C 511C.... THIS BUFFER IS NOW FULL AND MUST BE EMPTIED 512C 513 ISUM = ISUM + LDABUF 514 IDO = IDISK 515 CALL WRITDX(LUDA,IDISK,LBUF,ISCR(IBAD+1)) 516 ISCR(IBAD+LBUF-1) = IDO 517 ISCR(IBAD+LBUF) = 0 518 IDISK = IDISK + 1 519 END IF 520 END IF 521C 522C.... NOW CHECK RS 523C 524 IF (IPQ.EQ.IRS) GOTO 3 525 IF (IRS.GE.NST .AND. IRS.LE.NFIN) THEN 526 IRSS = IRS - NST + 1 527 IRSB = (IRSS+NPQ-1)/NPQ 528C 529C.... ALLOCATE TO APPROPRIATE BUFFER 530C 531 IBAD = IBUF(IRSB) 532 IRBAD = IBAD/IRAT 533 LENGTH = ISCR(IBAD+LBUF) + 1 534 SCR(IRBAD+LENGTH) = VALUE 535C ISCR(IBAD+IDABUF+LENGTH) = IRS*(2**16) + IPQ 536 ISCR(IBAD+IDABUF+LENGTH) = 537 * IOR( ISHFT( IRS , 16 ) , IPQ ) 538 ISCR(IBAD+LBUF) = LENGTH 539 IF (LENGTH .EQ. LDABUF) THEN 540C 541C.... BUFFER IS FULL AND MUST BE EMPTIED 542C 543 ISUM = ISUM + LDABUF 544 IDO = IDISK 545 CALL WRITDX(LUDA,IDISK,LBUF,ISCR(IBAD+1)) 546 ISCR(IBAD+LBUF-1) = IDO 547 ISCR(IBAD+LBUF) = 0 548 IDISK = IDISK + 1 549 END IF 550 END IF 551 3 CONTINUE 552C 553C.... THIS COMPLETES THE LOOP OVER THIS BUFFER OF AO INTEGRALS. 554C.... LOOP BACK TO READ ANOTHER 555C 556 GOTO 1 557 2 CONTINUE 558C 559C.... AT THIS POINT END-OF-FILE HAS BEEN ENCOUNTERED ON THE AO 560C... INTEGRAL FILE, THUS SIGNIFYING THE END OF THIS PASS. 561C.... IT IS NECESSARY TO WRITE LAST DA BACK-CHAIN ADDRESSES TO 562C.... THE LAST ADDRESS TABLE, EMPTYING THE BUFFERS IF THEY STILL 563C.... CONTAIN INFORMATION. 564C 565 KAP = -LBUF 566 IOFF = (ISQ-1)*NBATCH 567 DO 6 K = 1,NOB 568 KAP = KAP + LBUF 569 IDO = IDISK 570 IF (ISCR(KAP+LBUF) .NE. 0) THEN 571C 572C.... THIS BUFFER CONTAINS INFORMATION SO IT IS EMPTIED 573C 574 ISUM = ISUM + ISCR(KAP+LBUF) 575 CALL WRITDX(LUDA,IDISK,LBUF,ISCR(KAP+1)) 576 LASTAD(IOFF+K) = IDO 577 IDISK = IDISK + 1 578 ELSE 579 IDO = ISCR(KAP+LBUF-1) 580 LASTAD(IOFF+K) = IDO 581 END IF 582 6 CONTINUE 583C 584C.... END OF THIS PASS. CHECK TIMING 585C 586100 CONTINUE 587 CALL GPCLOSE(LUINTA,'KEEP') 588 IF (JTER.EQ.1) THEN 589 PNZ = (LPQRS*2-NMBASX) * 100 590 PNZ = ISUM / PNZ 591 WRITE (LUPRI,202) NDAREC,(IDISK-1),PNZ 592 WRITE (LUPRI,203) THRP 593 END IF 594 deallocate(IINDX4) 595 CALL QEXIT('SORTA ') 596 RETURN 597C 598C.... FORMAT STATEMENTS 599C 600 200 FORMAT(/' FIRST HALF-SORT STATISTICS' 601 1 /' NUMBER OF PASSES OVER AO FILE',T33,I4 602 2 /' NUMBER OF BLOCKS PER CHAIN',T30,I7 603 3 /' TOTAL NUMBER OF CHAINS',T30,I7 604 4 /' NUMBER OF BUFFERS PER PASS',T30,I7/) 605 201 FORMAT(' NUMBER OF INTEGRALS PER DIRECT ACCESS BUFFER',I6) 606 202 FORMAT(/' NUMBER OF DA RECORDS FOR ALL INTEGRALS',I6 607 1 /' NUMBER OF DA RECORDS ACTUALLY USED ',I6 608 2 /' PERCENT NON-ZERO AO INTEGRALS',F15.2) 609 203 FORMAT(/' THRESHOLD FOR DISCARDING INTEGRALS',1P,D10.2) 610 END 611C /* Deck sortb */ 612 SUBROUTINE SORTB(ITRLVL,LUHALF,LUDA2,RINT,IINT,SCR,ISCR,MEMSX, 613 & MEMTX,IT) 614C 615C.... THIS ROUTINE SORTS HALF-TRANSFORMED INTEGRALS FOR THE SECOND 616C.... HALF-TRANSFORMATION. AS IN *SORTA*, THE SORT IS OF MULTI-PASS 617C.... TYPE, AND THE OPTIMIZING ALGORITHM IS AGAIN DESIGNED TO 618C.... MINIMIZE THE NUMBER OF DISK ACCESSES. 619C 620C.... NOTE THAT THIS VERSION IS DESIGNED TO COPE ONLY WITH MO 621C.... INDICES OF THE FORM IJ, OR IA, WHERE I,J ARE OCCUPIED AND 622C.... A IS A VIRTUAL INDEX. 623C 624C.... FORMAL PARAMETERS AND EXTERNALS REQUIRED ARE THE SAME AS 625C.... FOR *SORTA* (Q.V.) 626C 627C.... PETER R. TAYLOR LUND MAR *79 628C 629#include "implicit.h" 630 PARAMETER (IBUFL = 500) 631#include "iratdef.h" 632 DIMENSION IINT(*),ISCR(*),SCR(*),RINT(*),IBUF(IBUFL) 633 INTEGER IT(*) 634 LOGICAL FEXIST 635#include "lbmxsq.h" 636C 637#include "priunit.h" 638#include "inttra.h" 639#include "inftra.h" 640#include "inftap.h" 641#include "ibtdef.h" 642C 643 CALL QENTER('SORTB ') 644C 645C.... SCRATCH CORE AVAILABILITY 646C 647 MEMS = MEMSX 648 MEMT = MEMTX 649C 650C.... NUMBER OF MO PAIRS 651C MAXD = number of orbitals for full integral transformation 652C else number of occupied orbitals. 653C 654C 655 IF (ITRLVL.EQ.0) THEN 656 NCDA = NMASHX 657 MAXD = MOCCT 658 ELSE IF (ITRLVL .EQ. 10) THEN 659 NCDA = NMORBX 660 MAXD = MORBT 661 ELSE 662 NCDA = NMOCCX + MOCCT*MSSHT 663 MAXD = MOCCT 664 END IF 665C 666C NBLOCK is max number of simult. distrib. in transf. step 667C NCHAIN is number of chains for this blocking in the 668C transformation step 669C 670 NBLOCK = MEMT/NMBASX 671 NCHAIN = 1 + (NCDA-1)/NBLOCK 672 NBLOCK = 1 + (NCDA-1)/NCHAIN 673C 674C.... DETERMINE NUMBER OF PASSES 675C 676C 677C here LDABUG = DA buffer length for NSTEP = 1 678C 679 LDABUG = (IRAT*MEMS) / ((1+IRAT)*NCHAIN) - 1 680 IF (LDABUG .GE. LDAMAX) THEN 681 NSTEP = 1 682 ELSE 683 LDAMIN = LDAMAX / 4 684 NSTEP = LDAMIN / LDABUG + 1 685 END IF 686C 687C.... NOW DETERMINE THE REMAINING DATA 688C 689 NBATCH = (NCHAIN+NSTEP-1)/NSTEP 690C 691C.... CHECK FOR OVERFLOW 692C 693 IF (NCHAIN .GT. MAXCHA) THEN 694 WRITE(LUPRI,200) NSTEP,NBLOCK,NCHAIN,NBATCH 695 WRITE(LUPRI,153) MAXCHA 696 CALL QUIT('TRACTL.SORTB: TOO MANY CHAINS') 697 END IF 698 153 FORMAT(/' TRACTL.SORTB: TOO MANY CHAINS; MORE THAN',I5) 699C 700#if defined (VAR_OLDCODE) 701 IF (NBATCH .GT. IBUFL) THEN 702 WRITE(LUPRI,200) NSTEP,NBLOCK,NCHAIN,NBATCH 703 WRITE(LUPRI,151) IBUFL 704 CALL QUIT('TRACTL.SORTB: TOO MANY BUFFERS') 705 END IF 706 151 FORMAT(/' TRACTL.SORTB: MORE THAN',I5,' BUFFERS PER PASS', 707 1 /' Increase IBUFL in SORTB.') 708#else 709 IF (NBATCH .GT. IBUFL) THEN 710 NSTEP = (NCHAIN-1)/IBUFL + 1 711 NBATCH = 1 + (NCHAIN-1)/NSTEP 712 END IF 713#endif 714C 715 IF(JTER.EQ.1)WRITE(LUPRI,200)NSTEP,NBLOCK,NCHAIN,NBATCH 716C 717C.... NO. OF INDICES PER CHAIN 718C 719 NCD = NBLOCK 720 NDIV = NBATCH*NBLOCK 721C 722C.... BUFFER LENGTH 723C 724 LDABUG = (MEMS/(1+NBATCH+NBATCH/IRAT)) - 1 725 LDABUG = MIN(LDABUG,LDAMAX) 726 LDABUG = (LDABUG/IRAT) * IRAT 727 IF (JTER.EQ.1) WRITE (LUPRI,201) LDABUG 728C 729 IDABUG = IRAT*LDABUG 730 LBUF = IDABUG + LDABUG + 2 731 NDAREC = LPQCD / LDABUG + 1 732C 733C delete any existing LUDA2 734C 735 CALL GPINQ('AOMOTRA.DA','EXIST',FEXIST) 736 IF (FEXIST) THEN 737 IF (LUDA2 .LE. 0) CALL GPOPEN(LUDA2,'AOMOTRA.DA','OLD', 738 & 'DIRECT','UNFORMATTED',LBUF,OLDDX) 739 CALL GPCLOSE(LUDA2,'DELETE') 740 END IF 741 CALL GPOPEN(LUDA2,'AOMOTRA.DA','NEW','DIRECT','UNFORMATTED', 742 & LBUF,OLDDX) 743C 744C.... PREPARE LOOK-UP TABLE OF BUFFER ADDRESSES 745C 746 KAP = 0 747 DO 25 K = 1,NBATCH 748 IBUF(K) = KAP 749 KAP = KAP + LBUF 75025 CONTINUE 751C 752C.... HALF-TRANSFORMED INTEGRAL FILE SPEC 753C 754 LRINT = IRAT*LBMXSQ 755 LENINT = LRINT + LBMXSQ + 2 756C 757C.... BEGIN PASSES 758C 759 ISUM = 0 760 IDISK = 1 761 NFIN = 0 762 DO 100 ISQ = 1,NSTEP 763 NST = NFIN + 1 764 NFIN = NFIN + NDIV 765 NOB = NBATCH 766C 767C.... CHECK TO SEE WHETHER IT IS NECESSARY TO RESET THESE VALUES 768C.... ON THE LAST PASS 769C 770 IF (NFIN.GT.NCDA) THEN 771 NFIN = NCDA 772 NOB = (NFIN-NST+NBLOCK)/NBLOCK 773 END IF 774C 775C.... NOB GIVES THE ACTUAL NUMBER OF BUFFERS USED IN THIS PASS 776C 777C.... INITIALIZE BUFFERS 778C 779 KAP = 0 780 DO 30 K = 1,NOB 781 KAP = KAP + LBUF 782 ISCR(KAP-1) = -1 783 ISCR(KAP) = 0 784 30 CONTINUE 785C 786C.... BEGIN READING THE INTEGRAL FILE 787 REWIND(LUHALF) 788 1 CONTINUE 789 CALL READI(LUHALF,LENINT,IINT) 790 LENGTH = IINT(LENINT-1) 791 IF(LENGTH.EQ.0)GOTO 1 792 IF(LENGTH.EQ.-1)GOTO 2 793C 794C.... LOOP OVER THE INTEGRALS IN THIS BUFFER 795C 796 DO 3 M = 1,LENGTH 797C 798C.... PICK UP INTEGRAL AND LABEL, UNPACK IPQ, IC, AND ID 799C 800 LABEL = IINT(M+LRINT) 801 VALUE = RINT(M) 802 IPQ=IAND(ISHFT(LABEL,-16),IBT16) 803 IC =IAND(ISHFT(LABEL, -8),IBT08) 804 ID =IAND( LABEL, IBT08) 805C 806C.... PACK IC AND ID TO THE INDEX FORM USED IN BUFFER ALLOCATION 807C 808 IF (IC.GT.MAXD) THEN 809 ICDX = MOCCT*IC - NMOCCX + ID 810 ELSE 811 ICDX = IT(IC) + ID 812 END IF 813 IF (ICDX.LT.NST.OR.ICDX.GT.NFIN) GOTO 3 814C 815C.... IN RANGE. PICK UP BUFFER ADDRESS 816C 817 ICDB = (ICDX - NST)/NCD + 1 818C 819C.... ALLOCATE INTEGRAL AND PACK LABEL 820C 821 IBAD = IBUF(ICDB) 822 IRBAD = IBAD/IRAT 823 LENGTH = ISCR(IBAD+LBUF) + 1 824 SCR(IRBAD+LENGTH) = VALUE 825C ISCR(IBAD+IDABUG+LENGTH) = IPQ*2**16 + ICDX 826 ISCR(IBAD+IDABUG+LENGTH) = 827 * IOR( ISHFT( IPQ , 16 ) , ICDX ) 828 ISCR(IBAD+LBUF) = LENGTH 829 IF (LENGTH.GE.LDABUG) THEN 830C 831C.... THIS BUFFER IS FULL AND MUST BE EMPTIED 832C 833 ISUM = ISUM + LDABUG 834 IDO = IDISK 835 CALL WRITDX(LUDA2,IDISK,LBUF,ISCR(IBAD+1)) 836 ISCR(IBAD+LBUF-1) = IDO 837 ISCR(IBAD+LBUF) = 0 838 IDISK= IDISK + 1 839 END IF 840 3 CONTINUE 841C 842C.... LOOP OVER THIS BUFFER COMPLETE. READ ANOTHER 843C 844 GOTO 1 845 2 CONTINUE 846C 847C.... END OF HALF-TRANSFORMED INTEGRAL FILE. EMPTY BUFFERS OF 848C.... LAST RECORDS 849C 850 KAP = -LBUF 851 IOFF = (ISQ-1)*NBATCH 852 DO 6 K = 1,NOB 853 KAP = KAP + LBUF 854 IF (ISCR(KAP+LBUF).NE.0) THEN 855C 856C.... NEEDS EMPTYING 857C 858 ISUM = ISUM + ISCR(KAP+LBUF) 859 IDO = IDISK 860 CALL WRITDX(LUDA2,IDISK,LBUF,ISCR(KAP+1)) 861 IDISK = IDISK + 1 862 ELSE 863 IDO = ISCR(KAP+LBUF-1) 864 END IF 865 LASTAE(IOFF+K) = IDO 866 6 CONTINUE 867C 868C.... CHECK TIMING AND THEN END THIS PASS 869C 870 100 CONTINUE 871 IF (JTER.EQ.1) THEN 872 PNZ = LPQCD * 100 873 PNZ = ISUM / PNZ 874 WRITE (LUPRI,202) NDAREC,(IDISK-1),PNZ 875 END IF 876 CALL QEXIT('SORTB ') 877 RETURN 878C 879C.... FORMAT STATEMENTS 880C 881 200 FORMAT(/,' SECOND HALF-SORT STATISTICS'/ 882 1 ' NUMBER OF PASSES OVER FILE',T33,I4/ 883 2 ' NUMBER OF BLOCKS PER CHAIN',T30,I7/ 884 3 ' TOTAL NUMBER OF CHAINS',T30,I7/ 885 4 ' NUMBER OF BUFFERS PER PASS',T30,I7/) 886 201 FORMAT(' NUMBER OF INTEGRALS PER DIRECT ACCESS BUFFER',I6) 887 202 FORMAT(/' NUMBER OF DA RECORDS FOR ALL INTEGRALS',I6 888 1 /' NUMBER OF DA RECORDS ACTUALLY USED ',I6 889 2 /' PERCENT NON-ZERO INTEGRALS',F18.2) 890 END 891C /* Deck tracd */ 892#include "vdir.h" 893 SUBROUTINE TRACD(ITRLVL,LUHALF,LUDA,CMOT,BUF,IBUF,OUTB,IOUTB,PQCD, 894 & PQRD,TRI,ITSMO,ITSAO,ITSX,IT) 895C 896C Revisions 897C 26-Jul-1984 Hans Joergen Aa. Jensen 898C 25-Nov-1984 hjaaj (corrected code for ITRLVL.eq.2) 899C 30-Jul-1986 hjaaj (use CMOT instead of CMO) 900C 10-Feb-1989 hjaaj (ITRLVL=10) 901C 902C CASSCF:TRANSFORMATION SECTION.TRANSFORMATION OF LAST TWO INDICES 903C 904C THIS SUBROUTINE PERFORMS THE TRANSFORMATION OF THE SECOND INDEX 905C PAIR IN THE TWO-ELECTRON INTEGRALS (PQ/RS). THE RESULT IS A LIST 906C OF HALF TRANSFORMED INTEGRALS (PQ/CD),STORED ON UNIT LUHALF. 907C 908C CALLED FROM TRACTL 909C 910C CMOT = CMO(transposed) 911C 912C ********** RELEASE 79 03 28 ********** 913C 914#include "implicit.h" 915 DIMENSION CMOT(*),BUF(*),IBUF(*),OUTB(*),IOUTB(*), 916 * PQCD(MORBT,*),PQRD(*),TRI(*) 917 INTEGER ITSMO(*), ITSAO(*), ITSX(*), IT(*) 918#include "iratdef.h" 919#include "lbmxsq.h" 920C 921#include "priunit.h" 922#include "inttra.h" 923#include "inftra.h" 924#include "inftap.h" 925C 926 LOGICAL SKIPR 927C 928#include "ibtdef.h" 929C 930 CALL QENTER('TRACD ') 931 IF (ITRLVL.EQ.2 .OR. ITRLVL.EQ.3 .OR. ITRLVL.EQ.6) THEN 932CHJ-840726: skip (some) inactive orbitals 933 NXISHT = MORBT-MSSHT-MASHT 934 ELSE 935 NXISHT = 0 936 END IF 937 IF (IPRTRA .GE. 100) THEN 938 WRITE (LUPRI,*) 'Test output from TRACD' 939 WRITE (LUPRI,*) 'ITRLVL =',ITRLVL 940 WRITE (LUPRI,*) 'NXISHT =',NXISHT 941 END IF 942C ***** AUXILIARY PARAMETERS ***** 943C ***** LOUT22: LENGTH OF BUFFER FOR PQCD 944C ***** LDA22 : LENGTH OF DIRECT ACCESS BUFFER 945C ***** LPQCD : NUMBER OF INTEGRALS (PQ/CD) 946 IOUT = 0 947 LOUT = LBMXSQ 948 LOUTI = IRAT*LOUT 949 LOUT2 = LOUTI + LOUT 950 LOUT21 = LOUT2 + 1 951 LOUT22 = LOUT2 + 2 952 IDABUF = IRAT*LDABUF 953 LDA2 = IDABUF + LDABUF 954 LDA21 = LDA2 + 1 955 LDA22 = LDA21 + 1 956 LPQCD = 0 957C 958C ***** REWIND FILE FOR HALF TRANSFORMED INTEGRALS LUHALF ***** 959C 960 REWIND(LUHALF) 961C 962C ***** START TRANSFORMATION ***** 963C ***** LOOP OVER ALL PAIRS PQ ***** 964C 965 NPQTOT = (MBAST*MBAST+MBAST)/2 966 LPQ = MIN(NPQ,NPQTOT) 967 LTRIPQ = LPQ*NMBASX 968 IF (LTRIPQ .GT. LTRI) CALL ERRWRK('TRACD for TRI(pq)',LTRIPQ,LTRI) 969C This will only happen if TRACTL is called with (much) less work 970C space than was available when SORTA was called (first call). 971 ICHAIN = 0 972 IPQ = 0 973 LPQ = NPQ 974 ISROLD = 0 975 SKIPR = .FALSE. 976 DO 100 NP = 1,MBAST 977 ISP = ITSAO(NP) 978 DO 100 NQ = 1,NP 979 ISQ = ITSAO(NQ) 980 ISPQ = MUL(ISP,ISQ) 981 IPQ = IPQ + 1 982 IPQL16 = ISHFT(IPQ,16) 983 LPQ = LPQ + 1 984 IF (LPQ .GT. NPQ) THEN 985 ICHAIN = ICHAIN + 1 986 LPQ = MIN(NPQ,NPQTOT+1-IPQ) 987 LTRIPQ = LPQ*NMBASX 988 CALL DZERO(TRI,LTRIPQ) 989 LPQ = 1 990C ***** FIND LAST ADDRESS OF NEW CHAIN AND START READ ***** 991 IADR = LASTAD(ICHAIN) 992 IF (IADR.LE.0) GO TO 100 993 IPQRS0 = NMBASX*(-1 - NPQ*(ICHAIN-1)) 994 103 CALL READDX(LUDA,IADR,LDA22,IBUF) 995 IADR = IBUF(LDA21) 996 LENGTH = IBUF(LDA22) 997 DO 104 I = 1,LENGTH 998 LDAI = IDABUF + I 999 IBL = IBUF(LDAI) 1000 MPQ = IAND(ISHFT(IBL,-16),IBT16) 1001 MRS = IAND( IBL, IBT16) 1002C ***** FIND ADDRESS AND DISTRIBUTE INTEGRAL ***** 1003 IPQRS = IPQRS0 + NMBASX*MPQ + MRS 1004 TRI(IPQRS) = BUF(I) 1005 104 CONTINUE 1006 IF (IADR.NE.-1) GO TO 103 1007 END IF 1008C 1009C ***** A NEW CHAIN IS NOW TRANSFERRED INTO CORE. START ***** 1010C ***** TRANSFORMATION OF LAST INDICES FOR THIS PQ ***** 1011C 1012C ***** SET THE ARRAY PQCD TO ZERO ***** 1013C ***** SET START ADDRESS IN TRI FOR INTEGRALS WITH ***** 1014C ***** FIRST PAIR INDEX EQUAL TO IPQ. ***** 1015C ***** LOOP OVER THIRD INDEX R ***** 1016C 1017 CALL DZERO(PQCD,MORBT*MORBT) 1018 NTPQ = NMBASX*(LPQ-1) 1019 DO 150 NR = 1,MBAST 1020 ISR = ITSAO(NR) 1021 IF (ISR .NE. ISROLD) THEN 1022 ISROLD = ISR 1023 SKIPR = .FALSE. 1024C ***** DETERMINE SYMMETRY FOR FOURTH INDEX ***** 1025 ISS = MUL(ISPQ,ISR) 1026 IF (ISS .GT. ISR) THEN 1027 IF (ITRLVL .EQ. 0 .OR. ITRLVL .EQ. 10) GO TO 149 1028 END IF 1029C 1030C ***** NUMBER OF BASIS FUNCTIONS AND ORBITALS ***** 1031 IF (ITRLVL .EQ. 10) THEN 1032 NENDD = MORB(ISS) 1033 ELSE 1034 NENDD = MOCC(ISS) 1035 END IF 1036 IF (ITRLVL.EQ.0) THEN 1037 MASHC = MASH(ISR) 1038 MASHD = MASH(ISS) 1039 IF (MASHC.EQ.0 .OR. MASHD.EQ.0) GO TO 149 1040C V------------------------------------------------ 1041 NENDC=MOCC(ISR) 1042 ELSE 1043 NENDC=MORB(ISR) 1044 IF (NENDC.EQ.0 .OR. NENDD.EQ.0) GO TO 149 1045C V------------------------------------------------ 1046 END IF 1047 MORBR=MORB(ISR) 1048 MORBS=MORB(ISS) 1049 MBASS=MBAS(ISS) 1050C ***** NUMBER OF ORBITALS IN PRECEDING SYMMETRIES *** 1051C ***** NUMBER OF BASIS FUNCTIONS OF EARLIER SYMMETRIES 1052 JORBC=JORB(ISR) 1053 JORBD=JORB(ISS) 1054 JBASR=JBAS(ISR) 1055 JBASS=JBAS(ISS) 1056C ***** SET LOOP LIMITS FOR ORBITAL C ***** 1057 NCRMAX=NENDC 1058 IF (ITRLVL.EQ.0) THEN 1059 NCRMIN = MISH(ISR)+1 1060 NDRMIN = MISH(ISS)+1 1061 ELSE IF (ITRLVL.EQ.10) THEN 1062 NCRMIN = 1 1063 NDRMIN = 1 1064 ELSE IF (ITRLVL.EQ.2 .OR. ITRLVL.EQ.3) THEN 1065 IF (ISS.GT.ISR) THEN 1066 NCRMIN = MOCC(ISR)+1 1067 NDRMIN = MISH(ISS)+1 1068C (we only need (ai/ integrals 1069C where ITSMO(a) = ITSMO(i) ) 1070 ELSE 1071 NCRMIN = 1 1072 NDRMIN = 1 1073C (we only need (ii/ integrals when 1074C ISP=ISQ=ISR=ISS, but we always 1075C need (iu/ and (ui/, so no gain here ) 1076 END IF 1077 ELSE 1078C -- ITRLVL.EQ.1 or 4 or 5 or 6: 1079 IF (ISS.GT.ISR) THEN 1080 NCRMIN = MOCC(ISR)+1 1081 ELSE 1082 NCRMIN = 1 1083 END IF 1084 NDRMIN = 1 1085 END IF 1086 NNDR = 1 + NENDD - NDRMIN 1087 IF (NCRMIN.GT.NCRMAX) GO TO 149 1088 IF (NNDR .LE. 0) GO TO 149 1089 ELSE IF (SKIPR) THEN 1090 GO TO 150 1091 END IF 1092C 1093C ***** START ADDRESSES FOR MOs ***** 1094C ***** SET THE ARRAY PQRD TO ZERO ***** 1095 IMOR = JCMO(ISR) + NCRMIN + (NR-JBASR-1)*MORBR 1096 IMOS = JCMO(ISS) + NDRMIN 1097 CALL DZERO(PQRD,NENDD) 1098C ***** LOOP OVER RELATIVE INDEX S IN SYMMETRY ISS ***** 1099 IPQRS1 = NTPQ + NR 1100 IPQRS2 = NTPQ + IT(NR) 1101 DO 130 NS = (JBASS+1),(JBASS+MBASS) 1102C ***** COMPUTE POSITION OF APPROPRIATE INTEGRAL ***** 1103 IF (NS.GT.NR) THEN 1104 IPQRS = IPQRS1 + IT(NS) 1105 ELSE 1106 IPQRS = IPQRS2 + NS 1107 END IF 1108 PQRS = TRI(IPQRS) 1109C ***** LOOP OVER MOs OF SYMMETRY ISS ***** 1110 IF (ABS(PQRS).GE.THRP) 1111 * CALL DAXPY(NNDR,PQRS,CMOT(IMOS),1,PQRD(NDRMIN),1) 1112 IMOS = IMOS + MORBS 1113 130 CONTINUE 1114C 1115C ***** PQRD NOW CONTAINS INTEGRALS (PQ/RD) FOR A GIVEN ***** 1116C ***** TRIPLE PQR AND ALL SYMMETRY ALLOWED VALUES OF D ***** 1117C ***** USE THESE TO COMPUTE CONTRIBUTIONS TO (PQ/CD) ***** 1118C 1119C 1120C ***** LOOP OVER ORBITALS C ***** 1121 DO 140 NCR = NCRMIN,NCRMAX 1122 CMOR = CMOT(IMOR) 1123 IMOR = IMOR + 1 1124 IF (ABS(CMOR) .LE. THRP) GO TO 140 1125 NC = NCR + JORBC 1126C ***** SET LOOP LIMITS FOR ORBITALS D ***** 1127 IF (ISS.EQ.ISR) THEN 1128 NDRMAX = MIN(NCR,NENDD) 1129 ELSE 1130 NDRMAX = NENDD 1131 END IF 1132 DO 135 NDR = NDRMIN,NDRMAX 1133 PQCD(JORBD+NDR,NC) = PQCD(JORBD+NDR,NC) 1134 * + CMOR*PQRD(NDR) 1135 135 CONTINUE 1136 140 CONTINUE 1137C END DO NCR 1138 GO TO 150 1139C 1140 149 CONTINUE 1141C ***** NOTHING NEEDED FROM THIS SYMMETRY OF R (C) 1142 SKIPR = .TRUE. 1143C 1144 150 CONTINUE 1145C END DO NR 1146C 1147C ***** ALL INTEGRALS (PQ/CD) FOR A GIVEN PAIR PQ HAVE ***** 1148C ***** BEEN CREATED. WRITE THEM ON UNIT LUHALF WITH ***** 1149C ***** INDICES IPQ AND ICD (REORDERED PAIR INDEX) ***** 1150C 1151 IF (ITRLVL.EQ.0) THEN 1152 ICMAX=MASHT 1153 ELSE 1154 ICMAX=MORBT 1155 END IF 1156 DO 170 IC = 1,ICMAX 1157 NC = ITSX( IC ) 1158 ISC = ITSMO( NC ) 1159 ISD = MUL(ISC,ISPQ) 1160 IF (ITRLVL.EQ.10) THEN 1161C ... full integral transformation 1162 IDMIN = 1 1163 IDMAX = IC 1164 ELSE IF (IC.GT.MOCCT) THEN 1165 IF (ITRLVL.EQ.2 .AND. ISC.EQ.ISP .AND. ISP.EQ.ISQ) THEN 1166 IDMIN = 1 1167C 1168C ( (ai/ distributions only needed for (ai/ai) integrals, 1169C ITRLVL.EQ.2, if ITRLVL.eq.3 then (ai/ai) are neglected 1170C if ITRLVL.eq.4 all (ai/bj) are included) 1171C 1172 ELSE 1173 IDMIN = NXISHT + 1 1174 END IF 1175 IDMAX = MIN(IC,MOCCT) 1176 ELSE 1177C -- IC .LE. MOCCT 1178 IF (IC.LE.NXISHT) THEN 1179 IF (ITRLVL.NE.2) GO TO 170 1180C ( if ITRLVL.eq.2 include (ii/aa) integrals ) 1181 IDMIN = IC 1182C ( only (ii/ inac-inac distributions needed ) 1183 ELSE 1184 IDMIN = 1 1185 END IF 1186 IDMAX = MIN(IC,MOCCT) 1187 END IF 1188 DO 160 ID = IDMIN,IDMAX 1189 ND = ITSX( ID ) 1190 IF (ITSMO( ND ) .NE. ISD) GO TO 160 1191 IF (NC .EQ. ND) THEN 1192 P = PQCD(NC,NC) 1193 ELSE 1194 P = PQCD(NC,ND) + PQCD(ND,NC) 1195 END IF 1196 IF (ABS(P).LT.THRP) GO TO 160 1197 LPQCD = LPQCD + 1 1198 IOUT = IOUT + 1 1199 IF (IOUT.GT.LOUT) THEN 1200C ***** WRITE THIS BUFFER ***** 1201 IOUT = 1 1202 IOUTB(LOUT21) = LOUT 1203 CALL WRITI(LUHALF,LOUT22,IOUTB) 1204 IF (IPRTRA .GE. 100) THEN 1205 WRITE(LUPRI,*) 1206 & 'Writing buffer for NP =',NP,' NQ = ', NQ 1207 WRITE(LUPRI,*) 'IC =',IC,' ID = ', ID 1208 WRITE(LUPRI,*) 'First 10 integrals' 1209 WRITE(LUPRI,'(5F14.7)') (OUTB(I),I=1,10) 1210 END IF 1211 END IF 1212 OUTB(IOUT) = P 1213C IOUTB(LOUTI+IOUT) = IPQ*2**16 + IC*2**8 + ID 1214 ICD = IOR(ISHFT(IC, 8),ID) 1215 IOUTB(LOUTI+IOUT) = IOR(IPQL16,ICD) 1216C ... IPQL16 = ISHFT(IPQ,16) 1217 160 CONTINUE 1218 170 CONTINUE 1219 100 CONTINUE 1220C END DO NQ 1221C END DO NP 1222C 1223C ***** WRITE LAST BUFFER ***** 1224C 1225 IOUTB(LOUT21)=IOUT 1226 CALL WRITI(LUHALF,LOUT22,IOUTB) 1227 IF (IPRTRA .GE. 100) THEN 1228 WRITE(LUPRI,*) 'Writing last buffer in TRACD' 1229 WRITE(LUPRI,*) 'First 10 integrals, length =',IOUT 1230 LPRI = MIN(10,IOUT) 1231 WRITE(LUPRI,'(5F14.7)') (OUTB(I),I=1,LPRI) 1232 END IF 1233C 1234C ***** WRITE END BUFFER ***** 1235C 1236 IOUTB(LOUT21)=-1 1237 CALL WRITI(LUHALF,LOUT22,IOUTB) 1238 REWIND(LUHALF) 1239 IF (JTER.EQ.1) THEN 1240 WRITE(LUPRI,1000) LUHALF,LPQCD 1241 WRITE(LUPRI,1100) THRP 1242 END IF 1243 1000 FORMAT(/,' NUMBER OF INTEGRALS (PQ/CD) WRITTEN ON UNIT', 1244 1 I3,' IS',I10) 1245 1100 FORMAT(/' THRESHOLD FOR DISCARDING INTEGRALS',1P,D10.2) 1246 CALL QEXIT('TRACD ') 1247 RETURN 1248 END 1249C /* Deck traab */ 1250#include "vdir.h" 1251 SUBROUTINE TRAAB(ITRLVL,LUDA2,CMOT,BUF,IBUF,OUTB,IOUTB,ABCD,PBCD, 1252 & TRI,ITSMO,ITSAO,ITSW,ITSX,IT) 1253C 1254C Revised for SIRIUS fall 1983, see below. 1255C Revisions: 1256C 25-Nov-1984 hjaaj (corrected code for ITRLVL.eq.2) 1257C 30-Jul-1986 hjaaj (use CMOT instead of CMO) 1258C CMOT = CMO(transposed) 1259C 28-Jun-1987 hjaaj (ITRLVL=0, removed AB .ge. CD restriction) 1260C 1261C TRAAB: TRANSFORMATION SECTION-TRANSFORM TWO FIRST INDICES 1262C 1263C THIS SUBROUTINE PERFORMS THE TRANSFORMATION OF THE FIRST 1264C INDEX PAIR IN THE TWO-ELECTRON INTEGRALS (PQ/CD).THE 1265C RESULT IS A LIST OF TRANSFORMED INTEGRALS (AB/CD),STORED 1266C ON UNIT LUINTM IN BUFFERS OF LENGTH 1344 1267C (671 INTEGRALS FOLLOWED BY 671 INDICES PLUS A POSITION 1268C GIVING THE NUMBER OF INTEGRALS IN THE BUFFER.LAST POSITION 1269C IS EMPTY). 1270C INDEX PACKING: NA*2**24+NB*2**16+NC*2**8+ND 1271C 1272C NOTE:CANONICAL ORDERING IS NOT ASSURED IN 1273C PARTIAL TRANSFORMATION 1274C 1275C ********** RELEASE 79 04 03 ********** 1276C 1277C Revision 14-Oct-1983 by Hans Jorgen Aa. Jensen, Uppsala 1278C This version writes also (AB/CD) where (AB) .lt. (CD), 1279C which makes the one-index transformation in NEO easier 1280C (could of course also be achieved with a SORTC similar to SORTA) 1281C Also start new record when (CD) increments, this means only the 1282C (CD) index of the first integral of each record needs to be 1283C checked and maybe decoded. 1284C 1285C Last revisions 5-Sep-1984 hjaaj / 1286C 27-Jul-1984 hjaaj / 1287C 23-Mar-1984 hjaaj / 1288C 11-Dec-1983 ha/hjaaj 1289C 1290#include "implicit.h" 1291 DIMENSION CMOT(*),BUF(*),IBUF(*),OUTB(*),IOUTB(*), 1292 * ABCD(MORBT,*),PBCD(*),TRI(*) 1293 INTEGER ITSMO(*), ITSAO(*), ITSW(*), ITSX(*), IT(*) 1294#include "iratdef.h" 1295C 1296#include "priunit.h" 1297#include "inforb.h" 1298#include "inttra.h" 1299#include "inftra.h" 1300#include "inftap.h" 1301C 1302 CHARACTER*8 TABLE(2),LAB123(3) 1303C 1304#include "ibtdef.h" 1305C 1306 DATA TABLE /'MOLTWOEL','END INTM'/ 1307 DATA LAB123/'********','********','********'/ 1308C 1309 CALL QENTER('TRAAB ') 1310 DST=SECOND() 1311 IF (IPRTRA .GE. 100) THEN 1312 WRITE (LUPRI,*) 'Test output from TRAAB' 1313 WRITE (LUPRI,*) 'ITRLVL =',ITRLVL 1314 END IF 1315C 1316 CALL GETDAT(LAB123(2),LAB123(3)) 1317C ... place date in LAB123(2) and time in LAB123(3) 1318C 1319 IF (IPRTRA .GE. 200) THEN 1320 WRITE (LUPRI,1008) 1321 1008 FORMAT(/,' *** TRAAB-INFORMATION, two-electron molecular ', 1322 * 'integrals'/2X,'Int. no.',4X,'C',4X,'D',4X,'A',4X,'B', 1323 * 10X,'Value') 1324 END IF 1325 IF (ITRLVL.EQ.2 .OR. ITRLVL.EQ.3 .OR. ITRLVL.EQ.6) THEN 1326 NXISHT = MORBT-MASHT-MSSHT ! skip inactive 1327 ELSE 1328 NXISHT = 0 1329 END IF 1330C ************** LENGTH OF BUFFER FOR ABCD on LUINTM 1331 LOUT =LBINTM 1332 LOUTI =IRAT*LOUT 1333 LOUT2 =LOUTI+LOUT 1334 LOUT21=LOUT2+1 1335 LOUT22=LOUT2+2 1336C ************* LENGTH OF DIRECT ACCESS BUFFER on LUDA2 1337 IDABUG=IRAT*LDABUG 1338 LDA2 =IDABUG+LDABUG 1339 LDA21 =LDA2+1 1340 LDA22 =LDA21+1 1341C ***** initialize counter for NUMBER OF INTEGRALS (AB/CD) 1342 LABCD =0 1343 IOUT =0 1344C 1345C LUINTM has been positioned before call of TRAAB (900801-hjaaj) 1346C 1347 WRITE(LUINTM)LAB123,TABLE(1) 1348C 1349 NRINTM = 0 1350 IF (LPQCD .EQ. 0) THEN 1351 LABCD = 0 1352 GO TO 8000 1353 END IF 1354C 1355C ***** START TRANSFORMATION ***** 1356C ***** LOOP OVER ALL PAIRS CD ***** 1357 ICD =0 1358 LTRICD=NCD*NMBASX 1359 IF (LTRICD .GT. LTRI) CALL ERRWRK('TRAAB for TRI(cd)',LTRICD,LTRI) 1360C should never happen. 1361 ICHAIN=0 1362 LCD =NCD 1363 IF (ITRLVL.EQ.0) THEN 1364 ICMAX=MASHT 1365 ELSE 1366 ICMAX=MORBT 1367 END IF 1368 IF (ITRLVL.LE.6) THEN 1369 MAXID = MOCCT 1370 ELSE 1371 MAXID = MORBT 1372C ... full integral transformation 1373 END IF 1374 DO 100 IC=1,ICMAX 1375 NC = ITSX(IC) 1376 ISC = ITSMO(NC) 1377 ITISC = IT(ISC) 1378 NCR=NC-JORB(ISC) 1379CHJ IDMAX=IC 1380CHJ IF (IC.GT.MOCCT) IDMAX=MOCCT 1381 IF (IC.LE.NXISHT) THEN 1382 IDMIN = IC 1383 IDMAX = IC 1384 ELSE 1385 IDMIN = 1 1386 IDMAX = MIN(MAXID,IC) 1387 END IF 1388 DO 100 ID=1,IDMAX 1389 ICD=ICD+1 1390 LCD=LCD+1 1391 IF (LCD.LE.NCD) GO TO 110 1392C ***** ONE CHAIN COMPLETED. READ AND DISTRIBUTE A NEW CHAIN 1393 LCD=1 1394 ICHAIN=ICHAIN+1 1395C ***** FIND LAST ADDRESS OF NEW CHAIN AND START READ ***** 1396 IADR=LASTAE(ICHAIN) 1397 IF(IADR.LE.0) GO TO 100 1398 CALL DZERO(TRI,LTRICD) 1399 MCDOFF = 1 + NCD*ICHAIN - NCD 1400 103 CALL READDX(LUDA2,IADR,LDA22,IBUF) 1401 IADR=IBUF(LDA21) 1402 LENGTH=IBUF(LDA22) 1403 IF(LENGTH.EQ.0) GO TO 106 1404 DO 104 I=1,LENGTH 1405 IBL=IBUF(IDABUG+I) 1406 MPQ=IAND(ISHFT(IBL,-16),IBT16) 1407 MCD=IAND( IBL, IBT16) 1408C ***** FIND ADDRESS AND DISTRIBUTE INTEGRAL ***** 1409 IPQCD=NMBASX*(MCD-MCDOFF)+MPQ 1410 TRI(IPQCD)=BUF(I) 1411 104 CONTINUE 1412 106 IF(IADR.NE.-1) GO TO 103 1413CHJ 1 ! IADR = 1 means non-zero integrals in this chain 1414 IADR = 1 1415C ***** A NEW CHAIN IS NOW TRANSFERRED INTO CORE. START ***** 1416C ***** TRANSFORMATION OF FIRST INDICES FOR THIS CD ***** 1417 110 CONTINUE 1418CHJ 2 1419 IF (IADR.LE.0) GO TO 100 1420CHJ-S-841018 1421C This patch (instead of DO 100 ID = IDMIN,IDMAX) was necessary 1422C to get right CHAIN in core (keep ICD and LCD in agreement 1423C with the counting in SORTB). 1424 IF (ID.LT.IDMIN) GO TO 100 1425CHJ-E-841018 1426 ND =ITSX(ID) 1427 ISD=ITSMO(ND) 1428CHJ-S-841125 1429 IF (IC.GT.MOCCT .AND. ID.LE.NXISHT) THEN 1430 IF (ISC.NE.ISD) GO TO 100 1431C ( we only need (ai/ distributions for (ai/ai), 1432C sym(a)=sym(i) ) 1433 END IF 1434CHJ-E-841125 1435 NDR = ND - JORB(ISD) 1436 ISCD=MUL(ISC,ISD) 1437 IF (ISD.LE.ISC) THEN 1438 NSCD=ITISC+ISD 1439 ELSE 1440 NSCD=IT(ISD)+ISC 1441 END IF 1442 IF (NC.GE.ND) THEN 1443C INDCD = NC*2**24 + ND*2**16 1444 INDCD = NC*2**8 + ND 1445 ELSE 1446C INDCD = ND*2**24 + NC*2**16 1447 INDCD = ND*2**8 + NC 1448 END IF 1449 INDCD = ISHFT( INDCD , 16 ) 1450C ***** SET THE ARRAY ABCD TO ZERO ***** 1451 CALL DZERO(ABCD,MORBT*MORBT) 1452C ***** SET START ADDRESS IN TRI FOR INTEGRALS WITH ***** 1453C ***** FIRST PAIR INDEX EQUAL TO ICD. ***** 1454 NTCD=NMBASX*(LCD-1) 1455C ***** LOOP OVER FIRST INDEX P ***** 1456 DO 150 NP=1,MBAST 1457 ISPA=ITSAO(NP) 1458C ***** DETERMINE SYMMETRY FOR SECOND INDEX ***** 1459 ISQB=MUL(ISPA,ISCD) 1460C ***** NUMBER OF ORBITALS AND BASIS FUNCTIONS ***** 1461 IF(ITRLVL.EQ.0.AND.MASH(ISQB).EQ.0)GO TO 150 1462 MORBA=MORB(ISPA) 1463 IF(MORBA.EQ.0) GO TO 150 1464 MORBB=MORB(ISQB) 1465 IF(MORBB.EQ.0) GO TO 150 1466 MOCCA=MOCC(ISPA) 1467 MOCCB=MOCC(ISQB) 1468 MORBP=MORB(ISPA) 1469 MORBQ=MORB(ISQB) 1470 MBASQ=MBAS(ISQB) 1471C ***** DETERMINE LOOP LIMITS FOR FIRST INDEX A ***** 1472 IF (ITRLVL .GE. 5) THEN 1473C ... general index on A and B. 1474 IF (ISPA.LT.ISQB) GO TO 150 1475 NARMIN = 1 1476 NARMAX = MORBA 1477 ELSE IF (IC.LE.MOCCT) THEN 1478C ***** OCCUPIED-OCCUPIED PAIR CD ***** 1479 IF(ISPA.LT.ISQB.AND.ITRLVL.NE.0) GO TO 150 1480C ... A .ge. B condition, 1481C for ITRLVL=0 A may be .lt. B when A secondary. 1482C (870628-hjaaj) 1483 IF (IC.LE.NXISHT) THEN 1484C ( we now know IC .eq. ID, 1485C this distribution is for (ii/aa) integrals ) 1486 IF (ITRLVL.NE.2 .OR. ISPA.NE.ISC) GO TO 150 1487C ( we only need (ii/aa) with sym(a) .eq. sym(i) ) 1488C ( if ITRLVL.eq.3 or 6, neglect (ii/aa) integrals ) 1489 NARMIN = MOCCA+1 1490 NARMAX = MORBA 1491 NBRMIN = NARMIN 1492 NBRMAX = NARMAX 1493 GO TO 145 1494 ELSE 1495 NARMIN=1 1496 NARMAX=MORBA 1497 END IF 1498 ELSE 1499C ***** VIRTUAL-OCCUPIED PAIR CD ***** 1500CHJ-S-841125 1501 IF (ID.LE.NXISHT) THEN 1502C ( symmetry tested above, this is (ai/, 1503C used for (ai/ai) ) 1504C ( if ITRLVL.eq.2 include (ai/ai) integrals ) 1505 IF (ITRLVL.NE.2 .OR. ISPA.NE.ISC) GO TO 150 1506 NARMIN = NCR 1507 NARMAX = NCR 1508 NBRMIN = NDR 1509 NBRMAX = NDR 1510 GO TO 145 1511 ELSE IF (ISPA.LT.ISQB) THEN 1512C ( (au/ case ) 1513 NARMIN = MOCCA+1 1514 ELSE 1515 NARMIN = 1 1516 END IF 1517CHJ-E-841125 1518 NARMAX=MORBA 1519#if defined (VAR_OLDCAS) 1520 IF (ITRLVL.LT.2) THEN 1521 IF (ISPA.LT.ISC) GO TO 150 1522 IF (ISPA.EQ.ISC) NARMIN=NCR 1523 END IF 1524#endif 1525 IF(NARMIN.GT.NARMAX) GO TO 150 1526 END IF 1527C 1528C *****SET ABSOLUTE LOOP LIMITS FOR B ***** 1529 IF (ITRLVL.GE.5) THEN 1530 NBRMIN = 1 1531 NBRMAX = MORBB 1532 ELSE IF (ITRLVL.EQ.0) THEN 1533 NBRMIN = MISH(ISQB) + 1 1534 NBRMAX = MOCCB 1535 ELSE 1536 NBRMIN = 1 1537 NBRMAX = MORBB 1538 IF (IC.GT.MOCCT) NBRMAX = MOCCB 1539 END IF 1540 IF (NBRMAX.LT.NBRMIN) GO TO 150 1541CHJ-S-841125 1542 145 CONTINUE 1543CHJ-E-841125 1544 NNBR = 1 + NBRMAX - NBRMIN 1545C ***** NUMBER OF ORBITALS OF PRECEDING SYMMETRIES ***** 1546 JORBA=JORB(ISPA) 1547 JORBB=JORB(ISQB) 1548 JBASP=JBAS(ISPA) 1549 JBASQ=JBAS(ISQB) 1550 NSAB = IT(MAX(ISPA,ISQB)) + MIN(ISPA,ISQB) 1551C ***** START ADDRESSES FOR MO*S ***** 1552C ***** SET THE ARRAY PBCD TO ZERO ***** 1553 IMOP = JCMO(ISPA) + NARMIN-1 + (NP-JBASP-1)*MORBP 1554 IMOQ = JCMO(ISQB) + NBRMIN 1555 CALL DZERO(PBCD(NBRMIN),NNBR) 1556 NCDITP = NTCD + IT(NP) 1557 NCDP = NTCD + NP 1558C ***** LOOP OVER RELATIVE INDICES Q IN SYMMETRY ISQB ***** 1559 DO 130 NQ = (JBASQ+1),(JBASQ+MBASQ) 1560C ***** FIND POSITION OF INTEGRAL (PQ/CD) ***** 1561 IF (NP.GE.NQ) THEN 1562 IPQCD = NCDITP + NQ 1563 ELSE 1564 IPQCD = NCDP + IT(NQ) 1565 END IF 1566 PQCD=TRI(IPQCD) 1567C ***** LOOP OVER MOs OF SYMMETRY ISQB ***** 1568C ***** AND TRANSFORM SECOND INDEX ***** 1569 IF (ABS(PQCD) .GE. THRP) 1570 * CALL DAXPY(NNBR,PQCD,CMOT(IMOQ),1,PBCD(NBRMIN),1) 1571 IMOQ = IMOQ + MORBQ 1572 130 CONTINUE 1573C ***** PBCD NOW CONTAINS INTEGRALS (PB/CD) FOR A ***** 1574C ***** GIVEN TRIPLE PCD AND ALL APPROPRIATE B ***** 1575C ***** LOOP OVER ORBITALS A AND ADD CONTRIBUTIONS***** 1576C ***** TO INTEGRALS (AB/CD) FOR THIS TRIPLE ***** 1577 DO 140 NAR=NARMIN,NARMAX 1578 IMOP = IMOP + 1 1579 CMOP = CMOT(IMOP) 1580 IF (ABS(CMOP) .LE. THRP) GO TO 140 1581C ***** SET ABSOLUTE INDEX ***** 1582 NA = NAR + JORBA 1583C ***** SET LOOP LIMITS FOR ORBITAL B ***** 1584 NBRS = NBRMIN 1585 NBRE = NBRMAX 1586 IF (ITRLVL.GE.2) THEN 1587 IF (ISPA.EQ.ISQB) NBRE = MIN(NBRE,NAR) 1588 ELSE IF (ITRLVL.EQ.0) THEN 1589 IF (ITSW(NA).LE.MASHT) THEN 1590C ... A is active 1591 IF(ISPA.LT.ISQB) GO TO 140 1592#if defined (VAR_OLDCAS) 1593 IF(NSAB.LT.NSCD) GO TO 140 1594 IF(ISPA.EQ.ISC.AND.NAR.LT.NCR)GO TO 140 1595 IF(NA.EQ.NC)NBRS=NDR 1596#endif 1597 IF(ISPA.EQ.ISQB)NBRE=NAR 1598 END IF 1599 ELSE 1600 IF(NSAB.LT.NSCD.AND.NAR.LE.MOCCA) NBRS=MOCCB+1 1601 IF(NA.EQ.NC) NBRS=NDR 1602 IF(ISPA.EQ.ISC.AND.NAR.LT.NCR) NBRS=MOCCB+1 1603 IF(ISPA.EQ.ISQB.AND.IC.LE.MOCCT) NBRE=NAR 1604C (HJ: why IC.le.MOCCT? I have removed it for ITRLVL.eq.2 or 3; 1605C IA.le.MOCCT would make more sense) 1606 END IF 1607CHJ-S-841125 1608 IF (IC.LE.NXISHT) THEN 1609C ( (ii/aa) case, symmetry etc. has already been checked) 1610 NBRS = NAR 1611 NBRE = NAR 1612 END IF 1613CHJ-E-841125 1614C ***** LOOP OVER ORBITALS B AND ADD CONTRIBUTION ***** 1615C ***** TO INTEGRALS (AB/CD) ***** 1616 DO 135 NBR = NBRS,NBRE 1617! that could be done with daxpy, could not it? 1618 ABCD(JORBB+NBR,NA) = ABCD(JORBB+NBR,NA) + CMOP*PBCD(NBR) 1619 135 CONTINUE 1620 140 CONTINUE 1621 150 CONTINUE 1622C ***** ALL INTEGRALS (AB/CD) FOR A GIVEN PAIR CD HAVE ***** 1623C ***** BEEN CREATED. WRITE THEM ON UNIT LUINTM WITH ***** 1624C ***** INDICES A,B,C AND D. ***** 1625CHJ: new index order: C,D; A,B 1626 DO 170 NA=1,MORBT 1627 INDCDA = IOR(INDCD , ISHFT( NA, 8 ) ) 1628 ISA = ITSMO( NA ) 1629 ISB = MUL( ISA, ISCD ) 1630 NBST = JORB(ISB) + 1 1631 NBEND = JORB(ISB) + MORB(ISB) 1632C ... not NBEND = MIN(NA,NBEND) 1633C because e.g. (ph/ph) integrals do not follow NB .le. NA 1634 DO 160 NB = NBST,NBEND 1635 P = ABCD(NB,NA) 1636 IF (ABS(P).LT.THRP) GO TO 160 1637 LABCD=LABCD+1 1638 IOUT =IOUT+1 1639 IF (IPRTRA .GE. 200) THEN 1640 WRITE (LUPRI,1010)LABCD,NC,ND,NA,NB,P 1641 1010 FORMAT(I10,4I5,F20.15) 1642 END IF 1643 IF (IOUT.GT.LOUT) THEN 1644 IOUT=1 1645 IOUTB(LOUT21)=LOUT 1646C ***** WRITE THIS BUFFER ***** 1647 CALL WRITI(LUINTM,LOUT22,IOUTB) 1648 NRINTM = NRINTM + 1 1649 IF (IPRTRA .GE. 100 .AND. IPRTRA .LT. 200) THEN 1650 WRITE(LUPRI,*) 1651 & ' Writing buffer for NC = ',NC,' ND = ',ND 1652 WRITE(LUPRI,*) ' NA = ',NA,' NB = ',NB 1653 WRITE(LUPRI,*) 'Buffer no.',NRINTM 1654 WRITE(LUPRI,*) ' First 10 integrals' 1655 WRITE(LUPRI,'(5F14.7)') (OUTB(I),I=1,10) 1656 END IF 1657 END IF 1658 OUTB(IOUT) = P 1659 IOUTB(LOUTI+IOUT) = IOR( INDCDA , NB ) 1660 160 CONTINUE 1661 170 CONTINUE 1662CHJ 3 *** Going to next (CD), empty this buffer *** 1663 IOUTB(LOUT21) = IOUT 1664 CALL WRITI(LUINTM,LOUT22,IOUTB) 1665 NRINTM = NRINTM + 1 1666 IF (IPRTRA .GE. 100 .AND. IPRTRA .LT. 200) THEN 1667 WRITE(LUPRI,*) 1668 & ' Writing last buffer for NC = ',NC,' ND = ',ND 1669 WRITE(LUPRI,*) 'Buffer no.',NRINTM 1670 WRITE(LUPRI,*) ' First 10 integrals, length =',IOUT 1671 LPRI = MIN(10,IOUT) 1672 WRITE(LUPRI,'(5F14.7)') (OUTB(I),I=1,LPRI) 1673 END IF 1674 IOUT = 0 1675 100 CONTINUE 1676C ***** WRITE LAST BUFFER ***** 1677 8000 CONTINUE 1678 IOUTB(LOUT21)=-1 1679 CALL WRITI(LUINTM,LOUT22,IOUTB) 1680 WRITE(LUINTM)LAB123,TABLE(2) 1681 REWIND(LUINTM) 1682 DFIN=SECOND() 1683 DTIM=DFIN-DST 1684 IF (JTER.EQ.1) THEN 1685 WRITE(LUPRI,49)DTIM 1686 WRITE(LUPRI,1000) LUINTM,LABCD,NRINTM,LBINTM 1687 WRITE(LUPRI,1100) THRP 1688 END IF 1689 49 FORMAT(' TIME IN TRAAB',F10.2) 1690 1000 FORMAT(/,' NUMBER OF TRANSFORMED INTEGRALS WRITTEN ON UNIT' 1691 1 ,I3,' IS',I10, 1692 2 /,' USING',I6,' + 1 RECORDS WITH BUFFER LENGTH',I6) 1693 1100 FORMAT(/' THRESHOLD FOR DISCARDING INTEGRALS',1P,D10.2) 1694 CALL QEXIT('TRAAB ') 1695 RETURN 1696 END 1697C /* Deck tractl */ 1698 SUBROUTINE TRACTL(KTRLVL,CMO,WORK,KWORK) 1699C 1700C Revisions 1701C 050125-hjaaj: new TRACTL_1 for dynamic allocation 1702C 1703C Input: 1704C KTRLVL, abs(KTRLVL) is transformation level, DELDA = (KTRLVL.lt.0) 1705C CMO, MO coefficients 1706C WORK, work array 1707C KWORK, actual length of work array 1708! 1709! module dependencies 1710 use lucita_mcscf_ci_cfg 1711C 1712#include "implicit.h" 1713#include "dummy.h" 1714 DIMENSION CMO(*), WORK(KWORK) 1715#include "iratdef.h" 1716#include "thrzer.h" 1717#include "lbmxsq.h" 1718C 1719C Used from common blocks: 1720C exeinf.h : FTRCTL, ITRLVL_LAST, LVLDRC_LAST 1721C 1722#include "priunit.h" 1723#include "inforb.h" 1724#include "exeinf.h" 1725#include "inttra.h" 1726#include "inftra.h" 1727#include "inftap.h" 1728#include "gnrinf.h" 1729C 1730 LOGICAL VFIRST, FOPEN, DELDA, FEXIST 1731 SAVE VFIRST, IPRTRA_SAVE, LDBUF, LRBUF 1732 DATA VFIRST /.TRUE./ 1733C 1734C VFIRST is is used to control setting of print level, see below. 1735C 1736C MAKE SURE BLOCK DATA MODULE FOR INFTRA IS INCLUDED : 1737C 1738 EXTERNAL BDTRA 1739C 1740 IF (VFIRST) THEN 1741C If SIRIUS has been called, the print level IPRTRA_SAVE will have the value 1742C from SIRIUS in all TRACTL calls. 1743C If SIRIUS has not been called before ABACUS is called, then 1744C the value specified in RHSINP for DERTRA is used. 1745C This makes sure that the SIRIUS value is 1746C always used if it has been defined, and the abacus 1747C value is different. 1748C 1749 IPRTRA_SAVE = IPRTRA 1750 VFIRST = .FALSE. 1751 ELSE 1752 IPRTRA = IPRTRA_SAVE 1753 END IF 1754C 1755C FTRCTL in exeinf.h forces integral sort and transformation, 1756C FTRCTL is true if Hermit has calculated a new AOTWOINT file 1757C since last call of TRACTL, or if integral type has changed 1758C (example: switch between AOSR2INT and AOTWOINT for srDFT) 1759C 1760 IF (FTRCTL) THEN ! new geometry, force new integral transformation 1761 ITRLVL_LAST = -999 1762 LVLDRC_LAST = -999 1763 LSRTAO = 0 ! force new integral sort and transformation in TRACTL_1 1764 FTRCTL = .FALSE. 1765 END IF 1766C 1767 ITRLVL = ABS( KTRLVL ) 1768C 1769 IF (IPRTRA .GE. 0) THEN 1770 CALL GETTIM(TCPU0,TWALL0) 1771 END IF 1772 1773C ***** SET THRESHOLDS ***** 1774C *********** THRP: THRESHOLD USED IN TRACD AND TRAAB 1775C *********** AND THRESHOLD USED IN SORTA FOR INTEGRALS (PQ/RS) 1776 IF (THRP.LE.1.D-16) THEN 1777 WRITE(LUPRI,'(/A,1P,D10.2,A,D10.2)') 1778 & ' INFO: Changed THRP threshold used in 2-el.'// 1779 & ' integral transformation from',THRP,' to',1.D-14 1780 THRP=1.D-14 1781 END IF 1782C 1783 IF (FCKTRA_TYPE .gt. 0) THEN 1784 NEWTRA_USEDRC = .FALSE. ! FCKTRA_DRCCTL not implemented yet 1785 CALL SIR_FCKTRA_CTL('COUL',KTRLVL,THRP,CMO,WORK,KWORK) 1786 ELSE IF (NEWTRA) THEN 1787 IF (ITRLVL .EQ. 0 .OR. ITRLVL .GE. 5) THEN ! zeroth, third or fourth order transformation 1788 NEWTRA_USEDRC = .FALSE. ! use DRCCTL and not Dirac integrals from N_TRACTL 1789 ! because a lot of MO integrals are calculated twice in N_TRACTL if ITRLVL .ge. 5 1790 ELSE 1791 NEWTRA_USEDRC = USEDRC 1792 END IF 1793 CALL N_TRACTL(KTRLVL,CMO,WORK,KWORK) 1794 IF ( NEWTRA_USEDRC ) USEDRC = .TRUE. 1795 ITRLVL_LAST = ITRLVL 1796 IF (NEWTRA_USEDRC) THEN 1797 IF (ITRLVL .EQ. 1 .OR. ITRLVL .EQ. 3) THEN ! see sirntra.F(TRALVL) 1798 LVLDRC_SAVE = 0 ! act-act Dirac integrals 1799 ELSE 1800 LVLDRC_SAVE = 1 ! occ-occ Dirac integrals 1801 END IF 1802 END IF 1803 ELSE 1804 NEWTRA_USEDRC = .FALSE. 1805 KFRSAV = 1 1806 KFREE = 1 1807 LFREE = KWORK 1808 CALL MEMGET('INTE',KITSMO,NBAST,WORK,KFREE,LFREE) 1809 CALL MEMGET('INTE',KITSAO,NBAST,WORK,KFREE,LFREE) 1810 CALL MEMGET('INTE',KITSW ,NBAST,WORK,KFREE,LFREE) 1811 CALL MEMGET('INTE',KITSX ,NBAST,WORK,KFREE,LFREE) 1812 CALL MEMGET('INTE',KIT ,NBAST+1,WORK,KFREE,LFREE) 1813 CALL TRACTL_1(KTRLVL,CMO,WORK(KITSMO),WORK(KITSAO), 1814 & WORK(KITSW),WORK(KITSX),WORK(KIT),WORK(KFREE),LFREE) 1815 CALL MEMREL('TRACTL_1',WORK,KFRSAV,KFRSAV,KFREE,LFREE) 1816 ! ITRLVL_LAST is set inside TRACTL_1, if integral transformation is not abandoned 1817 END IF 1818 1819! set flag for new integrals (required for parallel MCSCF/CI with lucita) 1820 mcscf_ci_update_ijkl = .true. 1821 1822 IF (IPRTRA .GE. 0) THEN 1823 CALL GETTIM(TCPU,TWALL) 1824 TCPU = TCPU -TCPU0 1825 TWALL = TWALL-TWALL0 1826 WRITE(LUPRI,2100) ITRLVL, TCPU, TWALL 1827 END IF 1828 2100 FORMAT(/' 2-el. integral transformation level ',I0,':', 1829 & ' Total CPU and WALL times (sec)',2F12.3) 1830C 1831C Sort mo integrals to Dirac form, if requested 1832C 1833 IF (IPRTRA .GE. 5) THEN 1834 WRITE(LUPRI,*) 'USEDRC, NEWTRA_USEDRC, ITRLVL:', 1835 & USEDRC, NEWTRA_USEDRC, ITRLVL 1836 END IF 1837 IF (USEDRC .AND. .NOT.NEWTRA_USEDRC .AND. ITRLVL .GT. 0) THEN 1838 CALL GETTIM(TCPU0,TWALL0) 1839 IF (ITRLVL .GE. 4) THEN 1840 LVLDRC = 1 1841C ... all occ-occ distributions 1842 ELSE 1843 LVLDRC = 0 1844C ... all active-active distributions 1845 END IF 1846 CALL DRCCTL(LVLDRC,CMO,WORK,KWORK) 1847 IF (IPRTRA .GE. 0) THEN 1848 CALL GETTIM(TCPU,TWALL) 1849 TCPU = TCPU -TCPU0 1850 TWALL = TWALL-TWALL0 1851 WRITE(LUPRI,2200) TCPU, TWALL 1852 END IF 1853 CALL FLSHFO(LUPRI) 1854 LVLDRC_LAST = LVLDRC 1855 END IF 1856 2200 FORMAT(/' Sorting integrals to Dirac format:', 1857 & ' Total CPU and WALL times (sec)',2F12.3) 1858 RETURN 1859 END 1860C /* Deck tractl_1 */ 1861 SUBROUTINE TRACTL_1(KTRLVL,CMO,ITSMO,ITSAO,ITSW,ITSX,IT, 1862 & WORK,KWORK) 1863C 1864C Revisions 1865C 18-Apr-1984 hjaaj / 4-Dec-1983 hjaaj 1866C 13-Dec-1984 hjaaj (ITRLVL=3 option) 1867C 28-Aug-1986 hjaaj (test on LTRI, ITRLVL=4 option) 1868C 18-May-1987 hjaaj (ITRLVL=5 option) 1869C 19-Dec-1988 hjaaj (ITRLVL,CMO in parameter list) 1870C 10-Feb-1989 hjaaj (ITRLVL=10, full integral transf.) 1871C 890615-hjaaj disabled if (first) for orb.spec., so it can be used 1872C in combination runs. 1873C 900801-hjaaj: added CMO to LUINTM 1874C 050125-hjaaj: TRACTL divided in TRACTL and TRACTL_1 for dynamic allocation 1875C of IT, ITSMO, ITSAO, ITSW, ITSX 1876C 1877C Input: 1878C KTRLVL, abs(KTRLVL) is transformation level, DELDA = (KTRLVL.lt.0) 1879C CMO, MO coefficients 1880C WORK, work array 1881C KWORK, actual length of work array 1882C MWORK, desired length of work array 1883C (dependent on size of physical memory) 1884C 1885C SUBROUTINE CALLS: 1886C SORTA (BIN SORT OF INTEGRALS (PQ/RS)) 1887C TRACD (TRANSFORMATION OF (PQ/RS) INTO (PQ/CD)) 1888C SORTB (BIN SORT OF INTEGRALS (PQ/CD)) 1889C TRAAB (TRANSFORMATION OF (PQ/CD) INTO (AB/CD)) 1890C MOLLAB 1891C 1892C Based on: 1893C CASSCF:TRANSFORMATION SECTION.CONTROL ROUTINE 1894C ********** RELEASE 79 04 20 ********** 1895C 1896#include "implicit.h" 1897#include "dummy.h" 1898 INTEGER ITSMO(*), ITSAO(*), ITSW(*), ITSX(*), IT(*) 1899 DIMENSION CMO(*), WORK(KWORK), MISH_TEST(8), MASH_TEST(8) 1900#include "iratdef.h" 1901#include "thrzer.h" 1902#include "lbmxsq.h" 1903 PARAMETER (DM1 = -1.0D0) 1904C 1905C Used from common blocks: 1906C exeinf.h : ITRLVL_LAST 1907C CBGETDIS : IAD* 1908C INFDIM : MWORK 1909C 1910#include "priunit.h" 1911#include "inforb.h" 1912#include "exeinf.h" 1913#include "infdim.h" 1914#include "inttra.h" 1915#include "inftra.h" 1916#include "inftap.h" 1917#include "gnrinf.h" 1918#include "cbgetdis.h" 1919C 1920 LOGICAL FOPEN, DELDA, FEXIST 1921 SAVE LDBUF, LRBUF 1922C 1923 CALL QENTER('TRACTL_1') 1924 ITRLVL = ABS(KTRLVL) 1925 DELDA = (KTRLVL .LT. 0) 1926 IF (IPRTRA .GT. 5) THEN 1927 WRITE(LUPRI,*) 'TRACTL_1(KTRLVL,CMO,WORK,KWORK,MWORK)' 1928 WRITE(LUPRI,*) '=====================================' 1929 WRITE(LUPRI,*) 'KTRLVL = ', KTRLVL 1930 WRITE(LUPRI,*) 'KWORK = ', KWORK 1931 WRITE(LUPRI,*) 'MWORK = ', MWORK 1932 END IF 1933 IF (IPRTRA .GE. 20) THEN 1934 WRITE(LUPRI,*) 'CMO array' 1935 CALL PRORB(CMO,.FALSE.,LUPRI) 1936 END IF 1937C 1938 IF (ITRLVL .NE. ITRLVL_LAST) THEN 1939 IF (IPRTRA .GT. 1) THEN 1940 JTER = 1 1941 ELSE 1942 JTER = 2 1943 END IF 1944 END IF 1945 LWORK = MIN(KWORK,MWORK) 1946 IF (JTER.EQ.1) WRITE(LUPRI,903) LWORK 1947 903 FORMAT(/' WORKING AREA IN TRANSFORMATION SECTION ',I10) 1948C 1949C ***** Reset /CBGETD/ because any H2AC on disk will be 1950C obsolete after transformation. 1951C 1952 IADINT = -1 1953 IADH2 = -1 1954 IADH2X = -1 1955 IADH2D = -1 1956 LUHALF = -1 1957 LUDA = -1 1958 LUDA2 = -1 1959C 1960C ***** TYPE OF TRANSFORMATION ***** 1961C ***** ITRLVL=0:FIRST INDEX,ALL ORBITALS ***** 1962C ***** OTHER INDICES,ACTIVE ORBITALS ONLY ***** 1963C (CAS-SCI, CI, MC gradient transformation) 1964C ***** ITRLVL=1:TWO INDICES,ALL ORBITALS ***** 1965C ***** TWO INDICES,PRIMARY ORBITALS ONLY ***** 1966C (CAS-NR transformation) 1967C 1968C ITRLVL=0: zeroth order transformation (uv|xy) 1969C (all four indices active orbitals) 1970C ITRLVL=1: first order transformation (gv|xy) 1971C (first index all orbitals, other indices active orbitals) 1972C ITRLVL=2: SIRIUS 2. order transformation, include (ii|aa) and (ia|ia) 1973C included: (uv|ab), (ua|vb), (ii|aa), and (ia|ia) 1974C (two indices active orbitals, other indices secondary orbitals) 1975C ITRLVL=3: SIRIUS 2. order transformation, neglect (ii|aa) and (ia|ia) 1976C included: (uv|ab) and (ua|vb) 1977C ITRLVL=4: full 2. order transformation (o1 o2|g1 g2) and (o1 g1|o2 g2) 1978C (two indices occupied orbitals, other indices general orbitals) 1979C 1980C ITRLVL=5: third order transformation (g1 g2|g3 o) 1981C (Three indices all orbitals, last index occupied orbitals) 1982C 1983C ITRLVL=6: third order transformation (g1 g2|g3 u) 1984C (Three indices all orbitals, last index active orbitals) 1985C 1986C ITRLVL=10: Full integral transformation (g1 g2|g3 g4) 1987C 1988 IF (IPRTRA.GT.0) WRITE (LUPRI,'(/A,I4)') 1989 & ' 2-ELECTRON INTEGRAL TRANSFORMATION SECTION; PARAMETER =', 1990 & ITRLVL 1991 IF (ITRLVL.GT.6 .AND. ITRLVL .NE. 10) THEN 1992 WRITE (LUPRI,'(//A)') 1993 1 ' TRACTL_1, parameter outside valid range' 1994 CALL QTRACE(LUPRI) 1995 CALL QUIT('ERROR, INVALID TRANSFORMATION PARAMETER (TRACTL_1)') 1996 END IF 1997C 1998C ***** Assign INTTRA from INFORB ***** 1999 MSYM = NSYM 2000 DO 5 ISYM = 1,8 2001 MISH(ISYM) = NISH(ISYM) 2002 MASH(ISYM) = NASH(ISYM) 2003 MSSH(ISYM) = NSSH(ISYM) 2004 MBAS(ISYM) = NBAS(ISYM) 2005 5 CONTINUE 2006 IF (IPRTRA .GE. 5) THEN 2007 WRITE(LUPRI,*) 'MSYM',MSYM 2008 WRITE(LUPRI,*) 'MISH',MISH(1:8) 2009 WRITE(LUPRI,*) 'MASH',MASH(1:8) 2010 WRITE(LUPRI,*) 'MSSH',MSSH(1:8) 2011 WRITE(LUPRI,*) 'MBAS',MBAS(1:8) 2012 END IF 2013C 2014C If first call of TRACTL_1 set up orbital information in /INTTRA/ 2015C ... 890615 - hjaaj: in new combination runs 2016C (e.g. RHF-MP2-CI) number of active orbitals will change 2017C therefore we now must recalculate this information 2018C each time. 2019C 2020C ***** SET UP ORBITAL DATA ***** 2021 MOCCT=0 2022 MSSHT=0 2023 MASHT=0 2024 MBAST=0 2025 NMORBT=0 2026 NMBAST=0 2027 DO 10 ISYM=1,NSYM 2028 MOCC(ISYM)=MISH(ISYM)+MASH(ISYM) 2029 MORB(ISYM)=MOCC(ISYM)+MSSH(ISYM) 2030 MOCCT=MOCCT+MOCC(ISYM) 2031 MSSHT=MSSHT+MSSH(ISYM) 2032 MASHT=MASHT+MASH(ISYM) 2033 MBAST=MBAST+MBAS(ISYM) 2034 NMORBT=NMORBT+(MORB(ISYM)**2+MORB(ISYM))/2 2035 NMBAST=NMBAST+(MBAS(ISYM)**2+MBAS(ISYM))/2 2036 10 CONTINUE 2037C 2038 MORBT=MOCCT+MSSHT 2039 IF (MORBT .GT. MBAST) CALL QUIT('MORBT .gt. MBAST :(') 2040CHJ-S-841208 2041 IF (MORBT.GT.255 .OR. MBAST.GT.255) THEN 2042 WRITE (LUPRI,*) 2043 * 'TRACTL_1-ERROR; TRACTL_1 is limited to 255 orbitals' 2044 WRITE (LUPRI,*) 2045 * 'NORBT =',MORBT,', NBAST =',MBAST 2046 CALL QTRACE(LUPRI) 2047 CALL QUIT('TRACTL_1: TOO MANY ORBITALS') 2048 END IF 2049CHJ-E-841208 2050 NMBASX=MBAST*(MBAST+1)/2 2051 NMORBX=MORBT*(MORBT+1)/2 2052 NMASHX=MASHT*(MASHT+1)/2 2053 M2ORBX=MORBT*MORBT 2054C ***** NUMBER OF ORBITALS AND BASIS FUNCTIONS OF EARLIER ***** 2055C ***** SYMMETRIES. START ADDRESSES FOR MO*S OF GIVEN ***** 2056C ***** SYMMETRY. ***** 2057 JBAS(1) = 0 2058 JORB(1) = 0 2059 JCMO(1) = 0 2060 DO 20 ISYM = 2,NSYM 2061 ISYM1 = ISYM - 1 2062 JBAS(ISYM) = JBAS(ISYM1) + MBAS(ISYM1) 2063 JORB(ISYM) = JORB(ISYM1) + MORB(ISYM1) 2064 JCMO(ISYM) = JCMO(ISYM1) + MORB(ISYM1)*MBAS(ISYM1) 2065 20 CONTINUE 2066C ***** COMPUTE SYMMETRY INDEX ITSAO and ITSMO ***** 2067 II = 0 2068 JJ = 0 2069 DO 30 ISYM = 1,NSYM 2070 MBASI = MBAS(ISYM) 2071 MORBI = MORB(ISYM) 2072 DO 25 I = 1,MBASI 2073 II = II + 1 2074 ITSAO(II) = ISYM 2075 IF (I.GT.MORBI) GO TO 25 2076 JJ = JJ + 1 2077 ITSMO(JJ) = ISYM 2078 25 CONTINUE 2079 30 CONTINUE 2080C 2081C MCMOT : SIZE OF THE MO-SPACE 2082C NNOCX : NUMBER OF OCCUPIED PAIRS 2083C MAXCHA: MAXIMUM NUMBER OF DA CHAINS 2084C LDAMAX: MAXIMUM SIZE OF DA RECORDS 2085C LRBUF : Max size of AO/MO integral records 2086C 2087 MCMOT = 0 2088 DO 50 ISYM = 1,NSYM 2089 50 MCMOT = MCMOT + MORB(ISYM)*MBAS(ISYM) 2090 NMOCCX = MOCCT*(MOCCT+1)/2 2091C 2092 LDAMAX = MAX(1706,NMBAST) 2093 LDAMAX = ((LDAMAX+IRAT-1)/IRAT)*IRAT 2094C ... 1706 gives buffer length of 20kB when IRAT=2 2095 LDBUF = (LDAMAX/IRAT)*(IRAT+1) + (1+1/IRAT) 2096 LRBUF = (LBMXSQ/IRAT)*(IRAT+1) + (1+1/IRAT) 2097C 2098C All orbital information now found and written in /INTTRA/. 2099C ... 890615 - hjaaj: in new combination runs 2100C (e.g. RHF-MP2-CI) number of active orbitals will change 2101C therefore we now must recalculate this information 2102C each time. 2103C 2104C ***** REORDERING INDICES FOR MOLECULAR ORBITALS ***** 2105C ***** ITSW REORDERS TO PRIMARY-SECONDARY ORDERING ***** 2106C ***** ITSX REORDERS BACK ***** 2107C 2108 IF (ITRLVL.EQ.0) GO TO 42 2109 IF (ITRLVL.EQ.10) GO TO 43 2110 IF (ITRLVL.EQ.1 .OR. ITRLVL.EQ.4 .OR. ITRLVL.EQ.5) GO TO 199 2111C 2112C --- ITRLVL.EQ.2 .OR. ITRLVL.EQ.3 .OR. ITRLVL.EQ.6: 2113C (inactive, active and secondary orbital lists) 2114C 2115 I = 0 2116 IO = 0 2117 IACT = MOCCT - MASHT 2118 ISEC = MOCCT 2119 DO 180 ISYM = 1,NSYM 2120 MISHI = MISH(ISYM) 2121 DO 120 K = 1,MISHI 2122 I = I + 1 2123 IO = IO + 1 2124 ITSW(I) = IO 2125 ITSX(IO) = I 2126 120 CONTINUE 2127 MASHI = MASH(ISYM) 2128 DO 140 K = 1,MASHI 2129 I = I + 1 2130 IACT = IACT + 1 2131 ITSW(I) = IACT 2132 ITSX(IACT) = I 2133 140 CONTINUE 2134 MSSHI = MSSH(ISYM) 2135 DO 160 K = 1,MSSHI 2136 I = I + 1 2137 ISEC = ISEC + 1 2138 ITSW(I) = ISEC 2139 ITSX(ISEC) = I 2140 160 CONTINUE 2141 180 CONTINUE 2142 GO TO 44 2143C 2144C --- ITRLVL.EQ.1 .OR. ITRLVL.EQ.4 .OR. ITRLVL.EQ.5 2145C (inactive orbitals in occupied list) 2146C 2147 199 CONTINUE 2148 I = 0 2149 ISEC = MOCCT 2150 IO = 0 2151 DO 40 ISYM = 1,NSYM 2152 MOCCI = MOCC(ISYM) 2153 DO 32 KOCC = 1,MOCCI 2154 I = I + 1 2155 IO = IO + 1 2156 ITSW(I) = IO 2157 ITSX(IO) = I 2158 32 CONTINUE 2159 MSSHI = MSSH(ISYM) 2160 DO 37 KSSH = 1,MSSHI 2161 I = I + 1 2162 ISEC = ISEC + 1 2163 ITSW(I) = ISEC 2164 ITSX(ISEC) = I 2165 37 CONTINUE 2166 40 CONTINUE 2167 GO TO 44 2168C 2169C --- ITRLVL.EQ.0: 2170C (inactive orbitals in secondary list) 2171C 2172 42 I = 0 2173 ISEC = MASHT 2174 IO = 0 2175 DO 400 ISYM = 1,NSYM 2176 MISHI = MISH(ISYM) 2177 DO 320 KISH = 1,MISHI 2178 ISEC = ISEC + 1 2179 I = I + 1 2180 ITSW(I) = ISEC 2181 ITSX(ISEC) = I 2182 320 CONTINUE 2183 MASHI = MASH(ISYM) 2184 DO 340 KASH = 1,MASHI 2185 IO=IO+1 2186 I=I+1 2187 ITSW(I)=IO 2188 ITSX(IO)=I 2189 340 CONTINUE 2190 MSSHI=MSSH(ISYM) 2191 DO 360 KSSH=1,MSSHI 2192 ISEC=ISEC+1 2193 I=I+1 2194 ITSW(I)=ISEC 2195 ITSX(ISEC)=I 2196 360 CONTINUE 2197 400 CONTINUE 2198 GO TO 44 2199C 2200 43 CONTINUE 2201C ... ITRLVL .EQ. 10 2202 DO 430 I = 1,MORBT 2203 ITSW(I) = I 2204 ITSX(I) = I 2205 430 CONTINUE 2206C 2207 44 CONTINUE 2208C 2209C ***** SET UP DYNAMIC STORAGE FOR TRACD ***** 2210C 2211 LW1 = 1 2212 LW2 = LW1 + MCMOT 2213 LW3 = LW2 + LDBUF 2214 LW4 = LW3 + LRBUF 2215 LW5 = LW4 + M2ORBX 2216 LW6 = LW5 + MORBT 2217 LTRI = LWORK - LW6 + 1 2218 NPQMIN = 1 + (NMBASX-1)/MAXCHA 2219 LTRIMIN = NPQMIN*NMBASX 2220 IF (LTRI.LT.LTRIMIN) THEN 2221 LWORK = LW6 - 1 + LTRIMIN 2222 IF (LWORK.LE.KWORK) THEN 2223 LTRI = LTRIMIN 2224 MWORK = MAX(LWORK,MWORK) 2225 WRITE(LUPRI,9030) LWORK 2226 ELSE 2227 WRITE(LUPRI,9000) LTRI 2228 WRITE(LUPRI,'(/A,6I9)') 2229 * ' LW1, LW2, ..., LW6 :',LW1,LW2,LW3,LW4,LW5,LW6 2230 CALL QTRACE(LUPRI) 2231 CALL QUIT('TRACTL_1: INSUFFICIENT SPACE FOR TRACD') 2232 END IF 2233 9030 FORMAT(/' TRACTL_1, not enough work space for TRACD,', 2234 * /' work space increased to',I10/) 2235 9000 FORMAT(/' WORK SPACE IN TRACD,',I10,', IS MUCH TOO SMALL.') 2236 END IF 2237C 2238C Allocate sorting area for SORTA/SORTB. 2239C hjaaj nov 2001: Now allow SORTA/SORTB to use all memory 2240C (i.e. KWORK) for sorting (LWORK is the memory available 2241C for TRACD, based on MWORK)) 2242C 2243 KWORK1 = 1 + LRBUF 2244 LWORK1 = KWORK - KWORK1 2245C 2246C calculate LPQRS, max. number of ao integrals 2247C 2248 LPQRS = 0 2249 DO 11 I = 1,NSYM 2250 DO 12 J = 1,I 2251 NSIJ = MUL(I,J) 2252 NBIJ = MBAS(I)*MBAS(J) 2253 IF (I.EQ.J) NBIJ = MBAS(I)*(MBAS(I)+1)/2 2254 DO 13 K = 1,I 2255 LMAX = K 2256 IF (K.EQ.I) LMAX = J 2257 DO 14 L = 1,LMAX 2258 NSKL = MUL(K,L) 2259 IF(NSKL.NE.NSIJ)GO TO 14 2260 NBKL = MBAS(K)*MBAS(L) 2261 IF(K.EQ.L)NBKL = MBAS(K)*(MBAS(K)+1)/2 2262 ITERM = NBIJ*NBKL 2263 IF (I.EQ.K.AND.J.EQ.L) ITERM = NBIJ*(NBIJ+1)/2 2264 LPQRS = LPQRS + ITERM 2265 14 CONTINUE 2266 13 CONTINUE 2267 12 CONTINUE 2268 11 CONTINUE 2269 IF (JTER.EQ.1) WRITE (LUPRI,902) LPQRS 2270 902 FORMAT(' MAXIMUM NUMBER OF AO INTEGRALS = ',I10) 2271C 2272 NBINTM = (NMORBT-1)/LBMXSQ + 1 2273 LBINTM = (NMORBT-1)/NBINTM + 1 2274 LBINTM = MAX(4,LBINTM) 2275C ... max distributions for any symmetry = 2276C number of distributions of symmetry 1 = NMORBT (890215-hjaaj) 2277C 2278C ***** CALLING SEQUENCE FOR FIRST SORT ***** 2279C 2280C.... LOOP TO SET UP TRIANGULAR INDEXING ARRAY 2281C 2282 ITLIM = MBAST + 1 2283 II = 0 2284 DO 45 I = 1,ITLIM 2285 IT(I) = II 2286 II = II + I 2287 45 CONTINUE 2288C 2289 DST = SECOND() 2290C 2291C 900801-hjaaj: check if transformation already done 2292C 971201-hjaaj: code moved to here (was too late) 2293C LSRTAO .ne. 0 should cover all situations where 2294C it is possible that MO integrals already exist 2295C on LUINTM. NOTE: this requires that SIR_INTOPEN has been 2296C called before first call to TRACTL_1 !!!!!! 2297C SIR_INTOPEN will set LSRTAO = -1 if LUINTM exists with 2298C a 'MOLTWOEL' label. Whenever TRACTL_1 has performed 2299C a transformation LSRTAO will be = +1. 2300C 2301C 2302 IF (LSRTAO .NE. 0) THEN 2303 REWIND(LUINTM) 2304 READ (LUINTM) 2305 READ (LUINTM) LBINTM_1,JTRLVL 2306 IF (JTRLVL .EQ. 3 .AND. ITRLVL .EQ. 2) JTRLVL = -1 ! (aa|ii) and (ai|ai) missing 2307 IF (JTRLVL .EQ. 2 .AND. ITRLVL .EQ. 3) JTRLVL = 3 2308 IF (JTRLVL .EQ. 5 .AND. ITRLVL .EQ. 6) JTRLVL = 6 2309 IF (JTRLVL .EQ. 6 .AND. ITRLVL .EQ. 4) JTRLVL = -1 ! (gg|ii) and (gi|gi) missing 2310 IF (JTRLVL .EQ. 6 .AND. ITRLVL .EQ. 5) JTRLVL = -1 ! (gg|gi) missing 2311 IF (JTRLVL .GE. ITRLVL) THEN 2312C read CMO matrix from LUINTM 2313C and subtract from input CMO matrix 2314 READ (LUINTM, IOSTAT=IOSVAL) WORK(1:NCMOT) 2315 IF (IOSVAL .EQ. 0) THEN 2316 CALL DAXPY(NCMOT,DM1,CMO,1,WORK,1) 2317 I = IDAMAX(NCMOT,WORK,1) 2318 READ (LUINTM) MISH_TEST(1:8), MASH_TEST(1:8) 2319 MISH_TEST(1:8) = MISH_TEST(1:8) - MISH(1:8) 2320 MASH_TEST(1:8) = MASH_TEST(1:8) - MASH(1:8) 2321 J = 0 2322 DO ISYM = 1,8 2323 J = J + ABS(MISH_TEST(ISYM)) + ABS(MASH_TEST(ISYM)) 2324 END DO 2325 IF (ABS(WORK(I)) .LE. THRZER .AND. J .EQ. 0) THEN 2326 IF (IPRTRA.GE.0) WRITE (LUPRI,'(/A)') 2327 & ' TRACTL_1: Integral transformation abandoned,'// 2328 & ' the required MO integrals are already available.' 2329 LBINTM = LBINTM_1 2330 GO TO 9999 2331 END IF 2332 END IF ! IOSVAL test 2333 END IF ! JTRLVL test 2334 END IF 2335 ITRLVL_LAST = ITRLVL 2336 2337 IF (LSRTAO .LT. 0) THEN 2338C ... check if old LUDA exist, if yes then check if file 2339C is consistent with LUINTM information /890302-hjaaj 2340 REWIND(LUINTM) 2341 READ (LUINTM) NPQ,LDABUF,(LASTAD(I),I=1,2500) 2342 REWIND(LUINTM) 2343 LBUF = (IRAT + 1)*LDABUF + 2 2344 CALL GPINQ('AOSRTINT.DA','EXIST',FEXIST) 2345 IF (FEXIST) THEN 2346 IF (LUDA .LE. 0) CALL GPOPEN(LUDA,'AOSRTINT.DA','OLD', 2347 & 'DIRECT',' ',LBUF,OLDDX) 2348C LUDA exists, try to read last record acc. to LUINTM ... 2349 IPQ = (NMBASX-1)/NPQ + 1 2350 IADR = LASTAD(1) 2351 DO 110 I = 2,IPQ 2352 IADR = MAX(IADR,LASTAD(IPQ)) 2353 110 CONTINUE 2354 LBUF = LBUF/IRAT 2355 READ (LUDA,REC=IADR,IOSTAT=IOSVAL) WORK(1:LBUF) 2356 IF (IOSVAL .NE. 0) THEN 2357C ... LUDA not consistent with LUINTM information 2358 CALL GPCLOSE(LUDA,'DELETE') 2359 LSRTAO = 0 2360 ELSE 2361 LSRTAO = 1 2362 WRITE(LUPRI,'(/A)') 2363 * ' Old direct access AO file found, SORTA skipped.' 2364 END IF 2365 ELSE 2366C ... LUDA does not exist 2367 LSRTAO = 0 2368 END IF 2369 END IF 2370 IF (LSRTAO.EQ.0) THEN 2371C hjaaj 11-jan-2007: make sure AOTWOINT exists 2372 CALL MAKE_AOTWOINT(WORK,KWORK) 2373 2374C HJ-840710 was: IF (JTER.EQ.1) THEN 2375 CALL GPINQ(FNINTM,'EXIST',FEXIST) 2376 if(FEXIST) CALL GPCLOSE(LUINTM,'DELETE') 2377 CALL GPOPEN(LUINTM,FNINTM,'UNKNOWN',' ','UNFORMATTED',IDUMMY, 2378 & .FALSE.) 2379C ... empty old mo integral file to save disk space during 2380C the transformation 2381C Now done by deleting the file, just in case it has been split 2382C K.Ruud, April 19 2000 2383C 2384 CALL SORTA(LUDA,WORK,WORK,WORK(KWORK1),WORK(KWORK1), 2385 * LWORK1,LTRI,IT) 2386C CALL SORTA(RINT,IINT,SCR,ISCR,MEMSX,MEMTX,IT) 2387 DFIN = SECOND() 2388 DTIM = DFIN - DST 2389 DST = DFIN 2390 IF (IPRTRA.GT.0) WRITE (LUPRI,46) DTIM 2391 46 FORMAT(/' TIME IN SORTA',F10.2) 2392 REWIND(LUINTM) 2393 WRITE (LUINTM) NPQ,LDABUF,(LASTAD(I),I=1,2500) 2394 REWIND(LUINTM) 2395 LSRTAO = 1 2396 ELSE 2397 REWIND(LUINTM) 2398 READ (LUINTM) NPQ,LDABUF,(LASTAD(I),I=1,2500) 2399 REWIND(LUINTM) 2400C 2401 LTRIPQ = NPQ*NMBASX 2402 LBUF = (IRAT + 1)*LDABUF + 2 2403 IF (LTRI .LT. LTRIPQ) THEN 2404C This will only happen if TRACTL_1 is called with (much) less 2405C work space than was available when SORTA was called. 2406 LTRI = KWORK - LW6 + 1 2407 IF (LTRI .LT. LTRIPQ) 2408 * CALL ERRWRK('TRACTL_1.TRACD, too large chains',LTRIPQ,LTRI) 2409 END IF 2410 IF (LUDA .LE. 0) CALL GPOPEN(LUDA,'AOSRTINT.DA','OLD', 2411 & 'DIRECT','UNFORMATTED',LBUF,OLDDX) 2412 END IF 2413C 2414 DSTTRA = DST 2415C 2416 CALL GPOPEN(LUHALF,' ',' ',' ','UNFORMATTED',IDUMMY,.FALSE.) 2417C 2418C 2419C ***** first transpose CMO 2420C 2421 DO 720 ISYM=1,NSYM 2422 MORBI=MORB(ISYM) 2423 IF(MORBI.EQ.0) GO TO 720 2424 MBASI=NBAS(ISYM) 2425 I1 = 1 + ICMO(ISYM) 2426 I2 = LW1 + JCMO(ISYM) 2427 CALL MTRSP(MBASI,MORBI,CMO(I1),MBASI,WORK(I2),MORBI) 2428 720 CONTINUE 2429C 2430C ***** CALLING SEQUENCE FOR FIRST TRANSFORMATION STEP ***** 2431C 2432 CALL TRACD(ITRLVL,LUHALF,LUDA,WORK(LW1),WORK(LW2),WORK(LW2), 2433 & WORK(LW3),WORK(LW3),WORK(LW4),WORK(LW5),WORK(LW6), 2434 & ITSMO,ITSAO,ITSX,IT) 2435 DFIN = SECOND() 2436 DTIM = DFIN - DST 2437 DST = DFIN 2438 IF (JTER.EQ.1) WRITE (LUPRI,47) DTIM 2439 47 FORMAT(' TIME IN TRACD',F10.2) 2440C 2441 IF (DELDA) THEN 2442 CALL GPCLOSE(LUDA,'DELETE') 2443 LSRTAO = 0 2444 ELSE 2445 CALL GPCLOSE(LUDA,'KEEP') 2446 END IF 2447C 2448 IF (LPQCD.EQ.0) THEN 2449 WRITE (LUPRI,'(/A)') ' --- NO NON-ZERO MO INTEGRALS ---' 2450 GO TO 7000 2451C 7000: CALL TRAAB TO OPEN LUINTM AND WRITE INFORMATION 2452C ABOUT NO INTEGRALS. 2453 END IF 2454C 2455C ***** SET UP DYNAMIC STORAGE FOR TRAAB ***** 2456C 870522: is now the same as for TRACD /hjaaj 2457C 2458C 2459C ***** CALLING SEQUENCE FOR SECOND SORT ***** 2460C 2461 CALL SORTB(ITRLVL,LUHALF,LUDA2,WORK,WORK,WORK(KWORK1), 2462 1 WORK(KWORK1),LWORK1,LTRI,IT) 2463 CALL GPCLOSE(LUHALF,'DELETE') 2464 DFIN = SECOND() 2465 DTIM = DFIN - DST 2466 IF (JTER.EQ.1) WRITE (LUPRI,48) DTIM 2467 48 FORMAT(' TIME IN SORTB',F10.2) 2468C 2469C ***** CALL SEQUENCE FOR SECOND TRANSFORMATION STEP ***** 2470C 2471 7000 CONTINUE 2472C 2473C ***** first transpose CMO 2474C 2475 DO 820 ISYM=1,NSYM 2476 MORBI=MORB(ISYM) 2477 IF(MORBI.EQ.0) GO TO 820 2478 MBASI=NBAS(ISYM) 2479 I1 = 1 + ICMO(ISYM) 2480 I2 = LW1 + JCMO(ISYM) 2481 CALL MTRSP(MBASI,MORBI,CMO(I1),MBASI,WORK(I2),MORBI) 2482 820 CONTINUE 2483C 2484C ***** save information about AO DA file (LUDA) on LUINTM 2485C 2486 REWIND(LUINTM) 2487 READ (LUINTM) 2488C 831014-hjaaj: maybe here also save LUINTM buffer length 2489C 860731-hjaaj: done 2490 WRITE (LUINTM) LBINTM,ITRLVL,NSYM,NORB,NBAS 2491 WRITE (LUINTM) CMO(1:NCMOT) 2492 WRITE (LUINTM) MISH(1:8), MASH(1:8) 2493C 2494 CALL TRAAB(ITRLVL,LUDA2,WORK(LW1),WORK(LW2),WORK(LW2), 2495 & WORK(LW3),WORK(LW3),WORK(LW4),WORK(LW5),WORK(LW6), 2496 & ITSMO,ITSAO,ITSW,ITSX,IT) 2497 IF (LPQCD.GT.0) CALL GPCLOSE(LUDA2,'DELETE') 2498C 2499 DFIN = SECOND() 2500 DTIM = DFIN - DSTTRA 2501 IF (IPRTRA .GT. 0) WRITE (LUPRI,49) DTIM 2502 49 FORMAT(' TOTAL TIME IN TRACD,SORTB,TRAAB',F10.2) 2503C 2504C set JTER to 2 to suppress printing on next call of TRACTL_1/hj-840710 2505 IF (IPRTRA .LT. 10) JTER = 2 2506C 2507C 2508 9999 CALL FLSHFO(LUPRI) 2509 CALL QEXIT('TRACTL_1') 2510 RETURN 2511C 2512C End of TRACTL_1 2513C 2514 END 2515C /* Deck bdtra */ 2516 BLOCK DATA BDTRA 2517C 2518C Define MUL in /INTTRA/ 2519C 2520#include "inttra.h" 2521C 2522C 2523C MULTIPLICATION TABLE FOR SYMMETRIES 2524C 2525 DATA MUL/1,2,3,4,5,6,7,8, 2526 * 2,1,4,3,6,5,8,7, 2527 * 3,4,1,2,7,8,5,6, 2528 * 4,3,2,1,8,7,6,5, 2529 * 5,6,7,8,1,2,3,4, 2530 * 6,5,8,7,2,1,4,3, 2531 * 7,8,5,6,3,4,1,2, 2532 * 8,7,6,5,4,3,2,1/ 2533C 2534 END 2535C /* Deck nxth2m */ 2536 SUBROUTINE NXTH2M(IC,ID,H2CD,NEEDTP,WRK,KFREE,LFREE,IDIST) 2537C 2538C Written by Hans Joergen Aa. Jensen October 1989 2539C This version is interface routine for old integral transformation. 2540C 2541C NOTE: The space allocated in WRK must not be touched outside 2542C until all desired distributions have been read. 2543C 2544C Purpose: 2545C Read next Mulliken two-electron integral distribution (**|cd) 2546C where (cd) distribution is needed according to NEEDTP(ITYPCD) 2547C (if needtp(itypcd) .gt. 0 all distributions of that type needed; 2548C if needtp(itypcd) .lt. 0 at least one distribution needed; 2549C if needtp(itypcd) .eq. 0 no distributions of that type needed). 2550C 2551C Usage: 2552C Set IDIST = 0 before first call of NXTH2M. 2553C DO NOT CHANGE IDIST or WRK(KFREE1:KFREE2-1) in calling routine 2554C until last distribution has been read (signalled by IDIST .eq. -1) 2555C Prototype code: 2556C IDIST = 0 2557C define NEEDTP(-4:6) 2558C 100 CALL NXTH2M(IC,ID,H2CD,NEEDTP,WRK,KFREE,LFREE,IDIST) 2559C IF (IDIST .GT. 0) THEN 2560C KW1 = KFREE 2561C LW1 = LFREE 2562C use (**|cd) distribution in H2CD as desired 2563C WRK(KW1:KW1-1+LW1) may be used 2564C GO TO 100 2565C END IF 2566C 2567C 2568#include "implicit.h" 2569#include "iratdef.h" 2570#include "priunit.h" 2571 DIMENSION H2CD(*),NEEDTP(-4:6),WRK(LFREE) 2572C 2573C Used from common blocks: 2574C INFORB : N2ORBX 2575C INFTAP : LUINTM,LBINTM 2576C 2577#include "inforb.h" 2578#include "inftap.h" 2579#include "inftra.h" 2580C 2581C 2582 SAVE KBUF, KIBUF, KNEXT 2583 DATA KNEXT /-1/ 2584C 2585 IF (FCKTRA_TYPE .gt. 0) THEN 2586 IF (LFREE .LT. N2ORBX) THEN 2587 write(lupri,*) 'KFREE, LFREE, IDIST',KFREE,LFREE,IDIST 2588 write(lupri,*) 'N2ORBX',N2ORBX 2589 call QUIT('too little memory for FCKTRA_NXTH2M') 2590 END IF 2591 CALL FCKTRA_NXTH2M(IC,ID,H2CD,NEEDTP,WRK(KFREE),IDIST) 2592 RETURN 2593 ELSE IF (NEWTRA) THEN 2594 CALL N_NXTH2M(IC,ID,H2CD,NEEDTP,WRK,KFREE,LFREE,IDIST) 2595 RETURN 2596 END IF 2597 CALL QENTER('NXTH2M') 2598C 2599C If first read (IDIST .eq. 0) then 2600C get buffer length and which distributions have been saved 2601C (lvlmol=0 for CI; lvlmol.ge.4 for inact-inact and inact-act; 2602C lvlmol=10 for full integral transformation) 2603C and allocate space for buffers. 2604C 2605 IF (IDIST .EQ. 0) THEN 2606 REWIND(LUINTM) 2607 READ (LUINTM,ERR=998,END=999) 2608 READ (LUINTM,ERR=998,END=999) LBINTM, LVLMOL 2609 LVLMIN = 0 2610 IF (NEEDTP(1).NE.0 .OR. NEEDTP(2).NE.0. .OR. 2611 & NEEDTP(4).NE.0 .OR. NEEDTP(5).NE.0) LVLMIN = 2 2612 IF (NEEDTP(1).GT.0 .OR. NEEDTP(4).GT.0) LVLMIN = 4 2613C level 4 is needed if all inac-inac or all inac-sec 2614C distributions must be available (needtp .gt. 0). 2615 IF (NEEDTP(6).NE.0) LVLMIN = 10 2616 IF (LVLMOL .LT. LVLMIN) THEN 2617 WRITE (LUPRI,*) 2618 & 'NXTH2M error, LVLMOL on LUINTM too small' 2619 WRITE (LUPRI,*) 'LVLMOL =',LVLMOL 2620 WRITE (LUPRI,*) 'need :',LVLMIN 2621 WRITE (LUPRI,*) 'NEEDTP :',(NEEDTP(J),J=-4,6) 2622 CALL QTRACE(LUPRI) 2623 CALL QUIT('NXTH2M error, LVLMOL too small') 2624 END IF 2625C 2626C We need to keep BUF and IBUF consecutively in memory. 2627C IBUF(LBINTM + 1) = LENGTH of record 2628C K.Ruud, April 2000 2629C 2630 LBINTT = LBINTM*(IRAT + 1) + 2 2631 CALL MEMGET('INTE',KBUF ,LBINTT,WRK,KFREE,LFREE) 2632 KIBUF = KBUF + LBINTM 2633 KNEXT = KFREE 2634 ELSE 2635C ... check that BUF,IBUF has not been destroyed by calling 2636C routine. 2637 IF (KNEXT.EQ. -1 ) THEN 2638 WRITE (LUPRI,*) 2639 & 'NXTH2M error, IDIST must be zero in first call' 2640 WRITE (LUPRI,*) 'IDIST =',IDIST 2641 CALL QTRACE(LUPRI) 2642 CALL QUIT('NXTH2M error, IDIST must be zero in first call') 2643 END IF 2644 IF (KFREE.LT.KNEXT) THEN 2645 WRITE (LUPRI,*) 2646 & 'NXTH2M error, KFREE lower than buffer allocation' 2647 WRITE (LUPRI,*) 'KFREE ',KFREE 2648 WRITE (LUPRI,*) 'KBUF ',KBUF 2649 WRITE (LUPRI,*) 'KIBUF ',KIBUF 2650 WRITE (LUPRI,*) 'KNEXT ',KNEXT, 2651 & ' ( next avail. address after buffers)' 2652 END IF 2653 CALL MEMCHK('IDIST .ne. 0 MEMCHK in NXTH2M',WRK,KBUF) 2654 IF (KFREE.LT.KNEXT) THEN 2655 CALL QTRACE(LUPRI) 2656 CALL QUIT('NXTH2M error:KFREE lower than buffer allocation') 2657 END IF 2658 END IF 2659C 2660 CALL NX2H2M(IC,ID,H2CD,NEEDTP,WRK(KBUF),WRK(KIBUF),IDIST) 2661C 2662C If no more buffers (IDIST .lt. 0) then release buffer space 2663C 2664 IF (IDIST .LT. 0) THEN 2665 CALL MEMREL('Releasing buffer space in NXTH2M',WRK,KBUF, 2666 & KBUF,KFREE,LFREE) 2667 KNEXT = -1 2668 END IF 2669 CALL QEXIT('NXTH2M') 2670 RETURN 2671 998 CONTINUE 2672 WRITE (LUPRI,*) 'ERROR reading MOTWOINT in NXTH2M' 2673 CALL QUIT('ERROR reading MOTWOINT in NXTH2M') 2674 999 CONTINUE 2675 WRITE (LUPRI,*) 'END OF FILE reading MOTWOINT in NXTH2M' 2676 CALL QUIT('END OF FILE reading MOTWOINT in NXTH2M') 2677 END 2678C /* Deck nx2h2m */ 2679 SUBROUTINE NX2H2M(IC,ID,H2CD,NEEDTP,BUF,IBUF,IDIST) 2680C 2681C Written by Hans Joergen Aa. Jensen October 1989 2682C 2683C Purpose: 2684C 2685C Read next Mulliken two-electron integral distribution (**|cd) 2686C where (cd) distribution is needed according to NEEDTP(ITYPCD) 2687C 2688C Input: 2689C NEEDTP(-4:6); non-zero for needed (cd) distribution types 2690C if .gt. 0 all distributions of that type needed 2691C IDIST; .eq. 0 first read 2692C .gt. 1 intermediate read 2693C .lt. 0 end-of-file has been reached previously 2694C Output: 2695C H2CD(NORBT,NORBT); H2CD(a,b) = (ab|cd) 2696C IC,ID; value of c and d 2697C IDIST; .gt. 0 when next distribution IC,ID available in H2CD 2698C = -1 when no more distributions 2699C Scratch: 2700C BUF, IBUF 2701C 2702C **************************************************************** 2703C 2704#include "implicit.h" 2705#include "iratdef.h" 2706#include "inftap.h" 2707 DIMENSION H2CD(NORBT,NORBT),BUF(LBINTM), IBUF(LBINTM) 2708#include "priunit.h" 2709 DIMENSION NEEDTP(-4:6) 2710C 2711C 2712C Used from common blocks: 2713C INFORB : NORBT 2714C INFIND : JTACT,JTSEC,ISW(*),... 2715C INFTAP : LUINTM,LBINTM 2716C 2717#include "maxash.h" 2718#include "maxorb.h" 2719#include "inforb.h" 2720#include "infind.h" 2721C 2722#include "orbtypdef.h" 2723C 2724 SAVE INDCDN, LENGTH 2725C 2726#include "ibtdef.h" 2727C 2728C **************************************************************** 2729C 2730C If first read (IDIST .EQ. 0) 2731C then setup for reading MO integrals ... 2732C 2733 IF (IDIST .EQ. 0) THEN 2734 REWIND(LUINTM) 2735 CALL MOLLAB('MOLTWOEL',LUINTM,LUPRI) 2736 100 CONTINUE 2737 READ (LUINTM) BUF, IBUF, LENGTH 2738 IF (LENGTH.EQ.0) GOTO 100 2739 IF (LENGTH.EQ.-1) THEN 2740 WRITE (LUPRI,'(/A)') 2741 & ' **** NXTH2M-ERROR *** no MO integrals on LUINTM' 2742 CALL QTRACE(LUPRI) 2743 CALL QUIT('*** ERROR ***(NXTH2M) no mo integrals on LUINTM') 2744 END IF 2745 INDCDN = IAND(ISHFT(IBUF(1),-16),IBT16) 2746 END IF 2747C 2748C *** Initialize H2CD 2749C 2750 CALL DZERO(H2CD,N2ORBX) 2751C 2752C *** Read next distribution which is needed according to NEEDTP(6) 2753C into H2CD 2754C 2755C ITYPCD values: 1=i*i : 2=t*i : 3=t*t : 4=a*i : 5=a*t : 6=a*a 2756C 0 for not wanted type. 2757C 2758C The CD distributions are stored by the present transformation 2759C program with IC.ge.ID 2760C 2761 200 CONTINUE 2762 IF (INDCDN .EQ. -1) THEN 2763C Last distribution has been read 2764 IDIST = -1 2765 GO TO 9999 2766 END IF 2767 INDCD = INDCDN 2768 IDIST = IDIST + 1 2769 IC = IAND(ISHFT(INDCD,-8),IBT08) 2770 ID = IAND( INDCD, IBT08) 2771 ITYPC = IOBTYP(IC) 2772 ITYPD = IOBTYP(ID) 2773 ITYPCD = IDBTYP(ITYPC,ITYPD) 2774 IF ( NEEDTP(ITYPCD) .EQ. 0 ) THEN 2775 ITYPCD = 0 2776 ELSE 2777 GO TO 400 2778C ... first buffer is already in BUF,IBUF 2779 END IF 2780C 2781 300 CONTINUE 2782C 2783 READ (LUINTM) BUF, IBUF, LENGTH 2784 IF (LENGTH.EQ.0) GO TO 300 2785 IF (LENGTH.EQ.-1) THEN 2786 INDCDN = -1 2787 ELSE 2788 INDCDN = IAND(ISHFT(IBUF(1),-16),IBT16) 2789 END IF 2790 IF (INDCDN.NE.INDCD) THEN 2791C ... finished reading the INDCD distribution 2792C if not needed (ITYPCD .eq. 0) go find type of 2793C new INDCDN distribution 2794 IF (ITYPCD.EQ.0) THEN 2795 GO TO 200 2796 ELSE 2797 GO TO 9999 2798 END IF 2799 END IF 2800 IF (ITYPCD.EQ.0) GO TO 300 2801C 2802 400 DO 450 I = 1,LENGTH 2803 IA = IAND(ISHFT(IBUF(I),-8),IBT08) 2804 IB = IAND( IBUF(I), IBT08) 2805 H2CD(IA,IB) = BUF(I) 2806 H2CD(IB,IA) = BUF(I) 2807 450 CONTINUE 2808 GO TO 300 2809C 2810C******************************************************************* 2811C 2812C End of subroutine NXTH2M 2813C 2814 9999 CONTINUE 2815 RETURN 2816 END 2817C /* Deck nxth2d */ 2818 SUBROUTINE NXTH2D(IC,ID,H2CD,NEEDTP,LUINTD,WRK,KFREE,LFREE,IDIST) 2819C 2820C Written by Hans Joergen Aa. Jensen October 1989 2821C This version is interface routine for old integral transformation. 2822C 2823C NOTE: The space allocated in WRK must not be touched outside 2824C until all desired distributions have been read. 2825C 2826C Purpose: 2827C Read next Dirac two-electron integral distribution <**|cd> 2828C where (cd) distribution is needed according to NEEDTP(ITYPCD) 2829C (if needtp(itypcd) .gt. 0 all distributions of that type needed; 2830C if needtp(itypcd) .lt. 0 at least one distribution needed; 2831C if needtp(itypcd) .eq. 0 no distributions of that type needed). 2832C 2833C Usage: 2834C Set IDIST = 0 before first call of NXTH2D. 2835C DO NOT CHANGE IDIST or WRK(KFREE1:KFREE2-1) in calling routine 2836C until last distribution has been read (signalled by IDIST .eq. -1) 2837C Prototype code: 2838C IDIST = 0 2839C define NEEDTP(-4:6) 2840C 100 CALL NXTH2D(IC,ID,H2CD,NEEDTP,WRK,KFREE,LFREE,IDIST) 2841C IF (IDIST .GT. 0) THEN 2842C KW1 = KFREE 2843C LW1 = LFREE 2844C use (**|cd) distribution in H2CD as desired 2845C WRK(KW1:KW1-1+LW1) may be used 2846C GO TO 100 2847C END IF 2848C 2849C 2850#include "implicit.h" 2851#include "priunit.h" 2852#include "dummy.h" 2853 DIMENSION H2CD(*),NEEDTP(-4:6),WRK(LFREE) 2854C 2855C Used from common blocks: 2856C INFORB : NCMOT 2857C INFTAP : LBINTD 2858C INFTRA : NEWTRA, NEWTRA_USEDRC, ? 2859C 2860#include "inforb.h" 2861#include "inftap.h" 2862#include "inftra.h" 2863C 2864 LOGICAL FEXIST 2865 INTEGER KBUF, KIBUF, KNEXT 2866 SAVE KBUF, KIBUF, KNEXT 2867 DATA KNEXT /-1/ 2868C 2869 IF (FCKTRA_TYPE .gt. 0 .AND. NEWTRA_USEDRC) THEN 2870 CALL FCKTRA_NXTH2D(IC,ID,H2CD,NEEDTP,IDIST) 2871 RETURN 2872 ELSE IF (NEWTRA .AND. NEWTRA_USEDRC) THEN 2873 CALL N_NXTH2D(IC,ID,H2CD,NEEDTP,WRK,KFREE,LFREE,IDIST) 2874 RETURN 2875 END IF 2876 CALL QENTER('NXTH2D') 2877C 2878C If first read (IDIST .eq. 0) then 2879C get buffer length and which distributions have been saved 2880C (lvldrc=0: act-act; =1: occ-occ; else orb-orb) 2881C and allocate space for buffers. 2882C 2883 IF (IDIST .EQ. 0) THEN 2884 CALL GPINQ('MODRCINT','EXIST',FEXIST) 2885 IF (FEXIST) THEN 2886 IF (LUINTD .LE. 0) CALL GPOPEN(LUINTD,'MODRCINT', 2887 & 'OLD',' ','UNFORMATTED',IDUMMY,.FALSE.) 2888 ELSE 2889 CALL QUIT('NXTH2D called, but MODRCINT does not exist.') 2890 END IF 2891 REWIND LUINTD 2892 CALL MOLLAB('DRCINFO ',LUINTD,LUPRI) 2893 READ (LUINTD) LBINTD, LVLDRC, NCMOT_OLD 2894 IF (NCMOT_OLD .NE. NCMOT) THEN 2895 WRITE (LUPRI,*) 2896 & 'NXTH2D error, wrong dimension of CMO array' 2897 WRITE (LUPRI,*) 'dim(CMO) on MODRCINT ',NCMOT_OLD 2898 WRITE (LUPRI,*) 'dim(CMO) now ',NCMOT 2899 CALL QUIT('NXTH2D error, old MODRCINT not correct') 2900 END IF 2901 LVLMIN = 0 2902 IF (NEEDTP(1).NE.0 .OR. NEEDTP(2).NE.0) LVLMIN = 1 2903 IF (NEEDTP(4).NE.0 .OR. NEEDTP(5).NE.0 .OR. 2904 & NEEDTP(6).NE.0) LVLMIN = 2 2905 IF (LVLDRC .LT. LVLMIN) THEN 2906 WRITE (LUPRI,*) 2907 & 'NXTH2D error, LVLDRC on LUINTD too small' 2908 WRITE (LUPRI,*) 'LVLDRC =',LVLDRC 2909 WRITE (LUPRI,*) 'need :',LVLMIN 2910 CALL QTRACE(LUPRI) 2911 CALL QUIT('NXTH2D error, LVLDRC too small') 2912 END IF 2913 CALL MEMGET2('REAL','BUF' ,KBUF ,LBINTD,WRK,KFREE,LFREE) 2914 CALL MEMGET2('INTE','IBUF',KIBUF,LBINTD,WRK,KFREE,LFREE) 2915 KNEXT = KFREE 2916 ELSE 2917C ... check that BUF,IBUF has not been destroyed by calling 2918C routine. 2919 IF (KNEXT.EQ. -1 ) THEN 2920 WRITE (LUPRI,*) 2921 & 'NXTH2D error, IDIST must be zero in first call' 2922 WRITE (LUPRI,*) 'IDIST =',IDIST 2923 CALL QTRACE(LUPRI) 2924 CALL QUIT('NXTH2D error, IDIST must be zero in first call') 2925 END IF 2926 IF (KFREE.LT.KNEXT) THEN 2927 WRITE (LUPRI,*) 2928 & 'NXTH2D error, KFREE lower than buffer allocation' 2929 WRITE (LUPRI,*) 'KFREE ',KFREE 2930 WRITE (LUPRI,*) 'KBUF ',KBUF 2931 WRITE (LUPRI,*) 'KIBUF ',KIBUF 2932 WRITE (LUPRI,*) 'KNEXT ',KNEXT, 2933 & ' ( next avail. address after buffers)' 2934 END IF 2935 CALL MEMCHK('IDIST .ne. 0 MEMCHK in NXTH2D',WRK,KBUF) 2936 IF (KFREE.LT.KNEXT) THEN 2937 CALL QTRACE(LUPRI) 2938 CALL QUIT('NXTH2D error:KFREE lower than buffer allocation') 2939 END IF 2940 END IF 2941 CALL NX2H2D(LUINTD,IC,ID,H2CD,NEEDTP,WRK(KBUF),WRK(KIBUF),IDIST) 2942C 2943C If no more buffers (IDIST .lt. 0) then release buffer space 2944C 2945 IF (IDIST .LT. 0) THEN 2946 CALL MEMREL('Releasing buffer space in NXTH2D',WRK,KBUF, 2947 & KBUF,KFREE,LFREE) 2948 KNEXT = -1 2949 END IF 2950 CALL QEXIT('NXTH2D') 2951 RETURN 2952 END 2953C /* Deck nx2h2d */ 2954 SUBROUTINE NX2H2D(LUINTD,IC,ID,H2CD,NEEDTP,BUF,IBUF,IDIST) 2955C 2956C Written by Hans Joergen Aa. Jensen October 1989 2957C 2958C Purpose: 2959C 2960C Read next Dirac two-electron integral distribution <**|cd> 2961C where (cd) distribution is needed according to NEEDTP(ITYPCD) 2962C 2963C Input: 2964C NEEDTP(i); non-zero for needed (cd) distribution types 2965C if .gt. 0 all distributions of that type needed 2966C IDIST; .eq. 0 first read 2967C .gt. 1 intermediate read 2968C .lt. 0 end-of-file has been reached previously 2969C Output: 2970C H2CD(NORBT,NORBT); H2CD(a,b) = <ab|cd> 2971C IC,ID; value of c and d 2972C IDIST; .gt. 0 when next distribution IC,ID available in H2CD 2973C = -1 when no more distributions 2974C Scratch: 2975C BUF, IBUF 2976C 2977C **************************************************************** 2978C 2979#include "implicit.h" 2980#include "priunit.h" 2981 REAL*8 H2CD(NORBT,NORBT),BUF(LBINTD) 2982 INTEGER*2 IBUF(2,LBINTD) 2983 INTEGER NEEDTP(-4:6) 2984C 2985C 2986C Used from common blocks: 2987C INFORB : NORBT 2988C INFIND : JTACT,JTSEC,ISW(*),... 2989C INFTAP : LUINTD,LBINTD 2990C 2991#include "maxash.h" 2992#include "maxorb.h" 2993#include "inforb.h" 2994#include "infind.h" 2995#include "inftap.h" 2996C 2997#include "orbtypdef.h" 2998C 2999 INTEGER*4 INDCDN, INDCD, LENGTH 3000 INTEGER*2 INDCD2(2) 3001 EQUIVALENCE (INDCD, INDCD2) 3002 SAVE INDCDN, LENGTH 3003C 3004#include "ibtdef.h" 3005C 3006C **************************************************************** 3007C 3008C If first read (IDIST .EQ. 0) 3009C then setup for reading MO integrals ... 3010C 3011 IF (IDIST .EQ. 0) THEN 3012 REWIND LUINTD 3013 CALL MOLLAB('DRCTWOEL',LUINTD,LUPRI) 3014 100 READ (LUINTD) BUF,IBUF,LENGTH,INDCDN 3015 IF (INDCDN.LT.0) THEN 3016 WRITE (LUPRI,'(/A)') 3017 & ' **** NXTH2D-ERROR *** no MO integrals on LUINTD' 3018 CALL QTRACE(LUPRI) 3019 CALL QUIT('*** ERROR ***(NXTH2D) no mo integrals on LUINTD') 3020 END IF 3021 IF (LENGTH.EQ.0) GOTO 100 3022 END IF 3023C 3024C *** Initialize H2CD 3025C 3026 CALL DZERO(H2CD,N2ORBX) 3027C 3028C *** Read next distribution which is needed according to NEEDTP(*) 3029C into H2CD 3030C 3031C ITYPCD values: 1=i*i : 2=t*i : 3=t*t : 4=a*i : 5=a*t : 6=a*a 3032C 0 for not wanted type. 3033C 3034C The CD distributions are stored by the present transformation 3035C program with IC.ge.ID 3036C 3037 200 CONTINUE 3038 IF (INDCDN .EQ. -1) THEN 3039C Last distribution has been read 3040 IDIST = -1 3041 GO TO 9999 3042 END IF 3043 INDCD = INDCDN 3044 IDIST = IDIST + 1 3045 IC = INDCD2(1) 3046 ID = INDCD2(2) 3047 ITYPC = IOBTYP(IC) 3048 ITYPD = IOBTYP(ID) 3049 ITYPCD = IDBTYP(ITYPC,ITYPD) 3050 IF ( NEEDTP(ITYPCD) .EQ. 0 ) THEN 3051 ITYPCD = 0 3052 ELSE 3053 GO TO 400 3054C ... first buffer is already in BUF,IBUF 3055 END IF 3056C 3057 300 READ (LUINTD) BUF,IBUF,LENGTH,INDCDN 3058 IF (INDCDN.NE.INDCD) THEN 3059C ... finished reading the INDCD distribution 3060C if not needed (ITYPCD .eq. 0) go find type of 3061C new INDCDN distribution 3062 IF (ITYPCD.EQ.0) THEN 3063 GO TO 200 3064 ELSE 3065 GO TO 9999 3066 END IF 3067 END IF 3068 IF (LENGTH.EQ.0) GO TO 300 3069 IF (ITYPCD.EQ.0) GO TO 300 3070C 3071 400 DO 450 I = 1,LENGTH 3072 IA = IBUF(1,I) 3073 IB = IBUF(2,I) 3074 H2CD(IA,IB) = BUF(I) 3075 450 CONTINUE 3076 GO TO 300 3077C 3078C******************************************************************* 3079C 3080C End of subroutine NX2H2D 3081C 3082 9999 CONTINUE 3083 RETURN 3084 END 3085! -- end of sirtra.F -- 3086