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 19C FILE: rspoit.F 20C 21C "oit" : One-Index Transformation 22C Everything is in MO basis in the routines in this file. 23C 24C /* Deck rspoit */ 25 SUBROUTINE RSPOIT(IDAXO,IDAXN,JTRLVL,LSYOIT,OITMAT, 26 * WRK,KFREE,LFREE,IPROIT) 27C 28C Copyright 25-Sep-1990 Hans Joergen Aa. Jensen 29C 30C Input: 31C IDAXO : file index for integrals to be one-index transformed 32C (IDAXO = 0 means untransformed integrals) 33C JTRLVL : range 0-4, number of general indices in transformed ints. 34C LSYOIT : Symmetry of one-index transformation 35C OITMAT : The transformation matrix 36C 37C Output: 38C IDAXN : file index for the transformed integrals 39C 40#include "implicit.h" 41 DIMENSION OITMAT(NORBT,NORBT), WRK(*) 42C 43C Used from common blocks: 44C INFORB : NORBT,N2ORBX 45C 46#include "priunit.h" 47#include "inforb.h" 48#include "infpri.h" 49C 50 EXTERNAL OITBD 51C 52 CALL QENTER('RSPOIT') 53C 54 IF (IPROIT .GE. 10) THEN 55 WRITE (LUPRI,'(//A/A)') 56 * ' Output from RSPOIT (one-index transformation)', 57 * ' ---------------------------------------------' 58 WRITE (LUPRI,'(/A,I5)') ' Print level :',IPROIT 59 CALL GETTIM(TSTRT,WSTRT) 60 END IF 61C 62C Check input 63C 64 NERR = 0 65 IF (IDAXO .LT. 0) THEN 66 NERR = NERR + 1 67 WRITE (LUPRI,*) 'RSPOIT: Illegal IDAXO =',IDAXO 68 END IF 69 IF (JTRLVL .LT. 0 .OR. JTRLVL .GT. 4) THEN 70 NERR = NERR + 1 71 WRITE (LUPRI,*) 'RSPOIT: Illegal JTRLVL =',JTRLVL 72 END IF 73 IF (LSYOIT .LT. 1 .OR. LSYOIT .GT. NSYM) THEN 74 NERR = NERR + 1 75 WRITE (LUPRI,*) 'RSPOIT: Illegal LSYOIT =',LSYOIT 76 END IF 77 IF (NERR .GT. 0) THEN 78 CALL QTRACE(LUERR) 79 CALL QUIT('Input errors in RSPOIT') 80 END IF 81C 82C Open integral files and assign IDAXN 83C 84 LRDAX = N2ORBT 85 IDAXN = 999 86 CALL OITOPN(IDAXO,IDAXN,LUDAXO,LUDAXN,JTRLVO,JTRLVL, 87 * LSYMXO,LSYOIT,LRDAX) 88C 89 IF (IPROIT .GE. 2) CALL FLSHFO(LUPRI) 90C 91C Do the one-index transformation 92C 93 IF (IDAXO .EQ. 0) THEN 94 CALL OIT2M(LUDAXN,JTRLVL,LSYOIT,OITMAT,WRK,KFREE,LFREE,IPROIT) 95 ELSE 96 CALL OIT2X(IDAXO,LUDAXN,JTRLVL,LSYMXO,LSYOIT,OITMAT, 97 * WRK,KFREE,LFREE,IPROIT) 98 END IF 99C 100 IF (IPROIT .GE. 10) THEN 101 CALL GETTIM(TEND,WEND) 102 WRITE (LUPRI,'(//A,F12.2,A)') 103 * ' CPU time used in RSPOIT :',TEND-TSTRT 104 WRITE (LUPRI,'(A,F12.2,A)') 105 * ' Elapsed time in RSPOIT :',WEND-WSTRT 106 END IF 107C 108 CALL QEXIT('RSPOIT') 109 RETURN 110 END 111C /* Deck oitbd */ 112 BLOCK DATA OITBD 113#include "cbdax.h" 114 DATA LUDAX /51,52,53,54,55/ 115 DATA IOPEN, IUSED, LSYDAX, LVLDAX 116 * /MXDAX*0, MXDAX*0, MXDAX*0, MXDAX*0/ 117 END 118C /* Deck oitopn */ 119 SUBROUTINE OITOPN(IDAXO,IDAXN,LUDAXO,LUDAXN,JTRLVO,JTRLVN, 120 * LSYMXO,LSYOIT,LRDAX) 121C 122C 11-Nov-1990 Hans Joergen Aa. Jensen 123C 124C If IDAXO .gt. 0 then open old unit. 125C IDAXO .eq. 0 refers to untransformed integrals. 126C If IDAXN .gt. 0 then open new unit. 127C 128#include "implicit.h" 129#include "iratdef.h" 130C 131C Used from common blocks: 132C CBDAX : * 133C INFORB : MULD2H 134C INFTAP : LUINTM 135C 136#include "cbdax.h" 137#include "priunit.h" 138#include "inforb.h" 139#include "inftap.h" 140C 141C 142 IDAXN1 = IDAXN 143 GO TO 5 144 ENTRY OITOPO(IDAXO,LUDAXO,JTRLVO,LSYMXO,LRDAX) 145 IDAXN1 = -1 146 5 CONTINUE 147 IF (IDAXO .GT. MXDAX) THEN 148 WRITE (LUPRI,*) 'OITOPN ERROR: IDAXO =',IDAXO 149 CALL QTRACE(LUPRI) 150 CALL QUIT('OITOPN ERROR: ILLEGAL IDAXO') 151 END IF 152 IF (IDAXN1 .LE. 0) GO TO 11 153 DO 10 I = 1,MXDAX 154 IF (IUSED(I) .EQ. 0) THEN 155 IDAXN = I 156 IUSED(IDAXN) = 1 157 GO TO 11 158 END IF 159 10 CONTINUE 160 WRITE (LUPRI,*) 'OITOPN ERROR: No more available units' 161 CALL QTRACE(LUPRI) 162 CALL QUIT('OITOPN : no more available units') 163 11 CONTINUE 164C 165 IF (IDAXO .GT. 0) THEN 166 LUDAXO = LUDAX(IDAXO) 167 LSYMXO = LSYDAX(IDAXO) 168 JTRLVO = LVLDAX(IDAXO) 169 IF (IUSED(IDAXO) .EQ. 0) THEN 170 WRITE (LUPRI,*) 'OITOPN ERROR: File with IDAXO =',IDAXO, 171 * ' has not been made.' 172 CALL QTRACE(LUPRI) 173 CALL QUIT('OITOPN : Old unit does not exist') 174 END IF 175 IF (IOPEN(IDAXO) .EQ. 0) THEN 176 CALL GPOPEN(LUDAX(IDAXO),' ','OLD','DIRECT','UNFORMATTED', 177 & IRAT*LRDAX,OLDDX) 178 IOPEN(IDAXO) = 1 179 END IF 180 ELSE IF (IDAXO .EQ. 0) THEN 181 LUDAXO = LUINTM 182 LSYMXO = 1 183 JTRLVO = 4 184 END IF 185C 186 IF (IDAXN1 .GT. 0) THEN 187 LUDAXN = -9999 188 CALL GPOPEN(LUDAXN,' ','NEW','DIRECT','UNFORMATTED',IRAT*LRDAX, 189 & OLDDX) 190 LUDAX(IDAXN) = LUDAXN 191 IOPEN(IDAXN) = 1 192 LSYDAX(IDAXN) = MULD2H(LSYMXO,LSYOIT) 193 LVLDAX(IDAXN) = JTRLVN 194 IF (IDAXO .GE. 0) THEN 195 JTRCHK = MIN(4,JTRLVN+1) 196 IF (JTRLVO .LT. JTRCHK) THEN 197 WRITE (LUPRI,*) 'OITOPN ERROR: File with IDAXO =', 198 * IDAXO,' has too low level.' 199 WRITE (LUPRI,*) ' Old level :',JTRLVO 200 WRITE (LUPRI,*) ' New level requested :',JTRLVN 201 CALL QTRACE(LUPRI) 202 CALL QUIT('OITOPN error,too low level on old LUDAX file') 203 END IF 204 END IF 205 END IF 206 RETURN 207 END 208C /* Deck oitclo */ 209 SUBROUTINE OITCLO(IDAXN,DISPOS) 210C 211C 16-Nov-1990 Hans Joergen Aa. Jensen 212C 213C Close unit with index IDAXN, if DISPOS(1:3) = 'DEL' then 214C delete the file. 215C 216#include "implicit.h" 217 CHARACTER*(*) DISPOS 218C 219C Used from common blocks: 220C CBDAX : LUDAX(*), IOPEN(*), ISUED(*) 221C 222#include "cbdax.h" 223#include "priunit.h" 224C 225 IF (IDAXN .LT. 1 .OR. IDAXN .GT. MXDAX) THEN 226 WRITE (LUPRI,*) 'OITCLO ERROR: IDAXN =',IDAXN 227 CALL QTRACE(LUPRI) 228 CALL QUIT('OITCLO ERROR: ILLEGAL IDAXN') 229 END IF 230 IF (IUSED(IDAXN) .EQ. 0) GO TO 9999 231 LUDAXN = LUDAX(IDAXN) 232 IF (DISPOS(1:3) .EQ. 'DEL') THEN 233 IF (IOPEN(IDAXN) .EQ. 0) THEN 234 CALL OITOPO(IDAXN,LUDAXN,JTRLVN,LSYMXN,LRDAX) 235C CALL OITOPO(IDAXO,LUDAXO,JTRLVO,LSYMXO,LRDAX) 236 END IF 237 CALL GPCLOSE(LUDAXN,'DELETE') 238 IUSED(IDAXN) = 0 239 ELSE IF (IOPEN(IDAXN) .NE. 0) THEN 240 CALL GPCLOSE(LUDAXN,'KEEP') 241 END IF 242 9999 IOPEN(IDAXN) = 0 243 CONTINUE 244 RETURN 245 END 246C /* Deck oith1 */ 247 SUBROUTINE OITH1(LSYOIT,OITMAT,H1,H1X,IH1SYM) 248C 249C Copyright 16-Nov-1990 Hans Joergen Aa. Jensen 250C 251C One-index transform H1(norbt,norbt) of symmetry IH1SYM 252C to H1X(norbt,norbt) using OITMAT(norbt,norbt) of symmetry 253C LSYOIT. 254C 255C Input : H1(a,b) of symmetry IH1SYM 256C Output: H1X(a,b) = H1X(a,b) 257C + H1(a,b) one-index transformed with OITMAT 258C = H1X(a,b) 259C + sum(c) [ OITMAT(a,c) H1(c,b) 260C - H1(a,c) OITMAT(c,b) ] 261C 262#include "implicit.h" 263 DIMENSION OITMAT(NORBT,NORBT), H1(NORBT,NORBT), H1X(NORBT,NORBT) 264C 265C Used from common blocks: 266C INFORB : NSYM, MULD2H, NORBT, NORB(*), IORB(*) 267C 268#include "inforb.h" 269C 270C 271C One-index transform both indices and add to previous content. 272C 273C (a~ b~) = SUM(c) [ OITMAT(a,c)*(c b) 274C - (a c) OITMAT(c,b) ] 275C 276 DO 100 ICSYM = 1,NSYM 277 IBSYM = MULD2H(ICSYM,IH1SYM) 278 IASYM = MULD2H(ICSYM,LSYOIT) 279 NORBC = NORB(ICSYM) 280 NORBB = NORB(IBSYM) 281 NORBA = NORB(IASYM) 282 IF ((NORBA.NE.0) .AND. (NORBC.NE.0) .AND. (NORBB.NE.0)) THEN 283 ICST = IORB(ICSYM) + 1 284 IBST = IORB(IBSYM) + 1 285 IAST = IORB(IASYM) + 1 286 CALL DGEMM('N','N',NORBA,NORBB,NORBC,1.D0, 287 & OITMAT(IAST,ICST),NORBT, 288 & H1(ICST,IBST),NORBT,1.D0, 289 & H1X(IAST,IBST),NORBT) 290 CALL DGEMM('N','N',NORBB,NORBA,NORBC,-1.D0, 291 & H1(IBST,ICST),NORBT, 292 & OITMAT(ICST,IAST),NORBT,1.D0, 293 & H1X(IBST,IAST),NORBT) 294 END IF 295 100 CONTINUE 296C 297C End of OITH1 298C 299 RETURN 300 END 301C /* Deck oitd1 */ 302 SUBROUTINE OITD1(LSYOIT,OITMAT,D1,D1X,ID1SYM) 303C 304C Copyright 8-Oct-2013 Hans Joergen Aa. Jensen 305C Purpose: One-index backtransformation of the one-electron density 306C matrix D1 using OITMAT. 307C Based on OITH1 308C Note that the transposed OITMAT is used here compared to in OITH1, 309C this is the only internal difference between the two subroutines. 310C 311C One-index transform D1(norbt,norbt) of symmetry ID1SYM 312C to D1X(norbt,norbt) using OITMAT(norbt,norbt) of symmetry LSYOIT. 313C 314C Input : D1(a,b) of symmetry ID1SYM 315C Output: D1X(a,b) = D1X(a,b) 316C + D1(a,b) one-index transformed with OITMAT 317C = D1X(a,b) 318C + sum(c) [ OITMAT(c,a) D1(c,b) 319C - D1(a,c) OITMAT(b,c) ] 320C 321C 322#include "implicit.h" 323 DIMENSION OITMAT(NORBT,NORBT), D1(NORBT,NORBT), D1X(NORBT,NORBT) 324C 325C Used from common blocks: 326C INFORB : NSYM, MULD2H, NORBT, NORB(*), IORB(*) 327C 328#include "inforb.h" 329C 330C 331C One-index transform both indices and add to previous content. 332C 333C (a~ b~) = SUM(c) [ OITMAT(c,a)*(c b) 334C - (a c) OITMAT(b,c) ] 335C 336 DO 100 ICSYM = 1,NSYM 337 IBSYM = MULD2H(ICSYM,ID1SYM) 338 IASYM = MULD2H(ICSYM,LSYOIT) 339 NORBC = NORB(ICSYM) 340 NORBB = NORB(IBSYM) 341 NORBA = NORB(IASYM) 342 IF ((NORBA.NE.0) .AND. (NORBC.NE.0) .AND. (NORBB.NE.0)) THEN 343 ICST = IORB(ICSYM) + 1 344 IBST = IORB(IBSYM) + 1 345 IAST = IORB(IASYM) + 1 346 CALL DGEMM('T','N',NORBA,NORBB,NORBC,1.D0, 347 & OITMAT(IAST,ICST),NORBT, 348 & D1(ICST,IBST),NORBT,1.D0, 349 & D1X(IAST,IBST),NORBT) 350 CALL DGEMM('N','T',NORBB,NORBA,NORBC,-1.D0, 351 & D1(IBST,ICST),NORBT, 352 & OITMAT(ICST,IAST),NORBT,1.D0, 353 & D1X(IBST,IAST),NORBT) 354 END IF 355 100 CONTINUE 356C 357C End of OITD1 358C 359 RETURN 360 END 361C /* Deck oit2m */ 362 SUBROUTINE OIT2M (LUDAXN,JTRLVL,LSYOIT,OITMAT, 363 * WRK,KFRSAV,LFRSAV,IPROIT) 364C 365C Copyright 15. Nov 1990 by Hans Jorgen Aa. Jensen. 366C 367C Purpose: 368C Construct one-index transformed 2-electron integrals. 369C 370C Input: 371C LUDAXN: unit number for output direct access file. 372C JTRLVL: Transformation level of transformed integrals. 373C LSYOIT: Symmetry of OITMAT 374C OITMAT: One-index transformation matrix 375C 376C Scratch: 377C WRK(KFRSAV:KFRSAV-1+LFRSAV) 378C 379C ********************************************************************* 380C 381#include "implicit.h" 382 DIMENSION OITMAT(N2ORBX), WRK(LFRSAV) 383C 384C Used from common blocks 385C INFORB : NSYM,N2ORBX,N2ORBT,NNORBX,NORBT,... 386C INFDIM : MWORK,NORBMA,... 387C 388#include "maxorb.h" 389#include "inforb.h" 390#include "infdim.h" 391#include "infpri.h" 392C 393C Local variables 394C 395 DIMENSION NEEDMU(-4:6), N2DIS(8) 396 397 CALL QENTER('OIT2M') 398 KFREE = KFRSAV 399 LFREE = LFRSAV 400C 401C Set NEEDMU array 402C ... occupied-occupied distributions always needed 403C ... unocc-occ distrib. needed for level 1 and higher 404C ... all distributions are needed for level 2 and higher 405C 406 NEEDMU(-4:6) = 0 407 IF (JTRLVL .GE. 2) THEN 408 NEEDMU(1:6) = 1 409 ELSE IF (JTRLVL .GE. 1) THEN 410 NEEDMU(1:5) = 1 411 ELSE 412 NEEDMU(1:3) = 1 413 END IF 414C 415C Determine if in core or out of core: 416C 417 MWORK1 = MIN(MWORK,LFREE-10*N2ORBX) 418C ... 10 is an arbitrary number, which is hoped to be sufficient 419 NH2XCD = MWORK1 / N2ORBT 420 IF (JTRLVL .EQ. 0 .OR. JTRLVL .EQ. 1) THEN 421 IF (NH2XCD .GE. N2OCCX) THEN 422 ICTOIT = 2 423 NH2XCD = N2OCCX 424 ELSE 425 ICTOIT = 3 426 IF (JTRLVL .EQ. 0) THEN 427 NH2XCD = NSYM*((NISHMA+NASHMA+1)*(NISHMA+NASHMA))/2 428 NH2XCD = MIN(NH2XCD,NNOCCX) 429C MAERKE need NNOCCT but that is not defined yet 430 ELSE 431 NH2XCD = NSYM*(NISHMA + NASHMA)*NORBMA 432 NH2XCD = MIN(NH2XCD,NNOCCX + NOCCT*NSSHT) 433C MAERKE need max of all symmetries but that is not defined 434 END IF 435 END IF 436 ELSE IF (JTRLVL .EQ. 2) THEN 437 IF (NH2XCD .GE. N2OCCX + 2*NOCCT*(NORBT-NOCCT) ) THEN 438 ICTOIT = 2 439 NH2XCD = N2OCCX + 2*NOCCT*(NORBT-NOCCT) 440 ELSE 441 ICTOIT = 3 442 NH2XCD = NNORBT 443 END IF 444 ELSE 445 IF (NH2XCD .GE. NNORBX) THEN 446 ICTOIT = 1 447 NH2XCD = NNORBX 448 ELSE 449 ICTOIT = 3 450 NH2XCD = NNORBT 451 END IF 452 END IF 453C 454C Allocate work memory 455C 456 IF (ICTOIT .EQ. 3) THEN 457 LH2X = MIN(N2ORBT,MWORK1/NH2XCD) 458 LH2X = MAX(NORBMA,LH2X) 459 ELSE 460 LH2X = N2ORBT 461 END IF 462 LH2XT = LH2X*NH2XCD 463 CALL MEMGET('REAL',KH2CD ,N2ORBX,WRK,KFREE,LFREE) 464 CALL MEMGET('REAL',KH2XCD,N2ORBX,WRK,KFREE,LFREE) 465 CALL MEMGET('INTE',KINDAB,N2ORBX,WRK,KFREE,LFREE) 466 CALL MEMGET('INTE',KINDCD,N2ORBX,WRK,KFREE,LFREE) 467 CALL MEMGET('INTE',KIN2CD,N2ORBX,WRK,KFREE,LFREE) 468 CALL MEMGET('INTE',KICDTR,N2ORBX,WRK,KFREE,LFREE) 469 CALL MEMGET('REAL',KH2X ,LH2XT ,WRK,KFREE,LFREE) 470C 471C Open auxiliary DA file 472C 473 IF (ICTOIT .EQ. 3) THEN 474 IDAXO = -1 475 IDAXT = 999 476 JTRLVT= 4 477 LSYMXO= 1 478 LRDAX = N2ORBT 479 CALL OITOPN(IDAXO,IDAXT,LUDAXO,LUDAXT,JTRLVO,JTRLVT, 480 * LSYMXO,LSYOIT,LRDAX) 481C CALL OITOPN(IDAXO,IDAXN,LUDAXO,LUDAXN,JTRLVO,JTRLVN, 482C * LSYMXO,LSYOIT,LRDAX) 483 ELSE 484 LUDAXT = -999 485 END IF 486C 487C Set INDAB and ICDTRA arrays 488C 489 CALL OITIND(WRK(KINDAB),N2DIS) 490 CALL OITICD(JTRLVL,WRK(KICDTR),WRK(KINDCD),IPROIT) 491C CALL OITICD(JTRLVL,ICDTRA,ITRTYP,IPROIT) 492C 493 CALL OIT2M2(LUDAXN,JTRLVL,ICTOIT,LUDAXT,LSYOIT, 494 * OITMAT,WRK(KH2CD),WRK(KH2XCD),WRK(KICDTR), 495 * NEEDMU,N2DIS,WRK(KINDAB),WRK(KINDCD),WRK(KIN2CD), 496 * WRK(KH2X),LH2X,NH2XCD, WRK,KFREE,LFREE,IPROIT) 497C 498 IF (ICTOIT .EQ. 3) CALL OITCLO(IDAXT,'DELETE') 499C 500C *** end of subroutine OIT2M 501C 502 CALL MEMREL('OIT2M',WRK,KFRSAV,KFRSAV,KFREE,LFREE) 503 CALL QEXIT('OIT2M') 504 RETURN 505 END 506C /* Deck oit2m2 */ 507 SUBROUTINE OIT2M2(LUDAXN,JTRLVL,ICTOIT,LUDAXT,LSYOIT, 508 * OITMAT,H2CD,H2XCD, 509 * ICDTRA,NEEDMU,N2DIS,INDAB,INDCD,IN2CD, 510 * H2X,LH2X,NH2XCD,WRK,KFRSAV,LFRSAV,IPROIT) 511C 512C Copyright 15. Nov 1990 by Hans Jorgen Aa. Jensen. 513C 514C Purpose: 515C Construct one-index transformed 2-electron integrals. 516C 517C Input: 518C LUDAXN: unit number for output direct access file. 519C JTRLVL: Transformation level of transformed integrals. 520C ICTOIT: control parameter for transformation 521C LUDAXT: temporary file for out-of-core transformation, if ICTOIT=3 522C LSYOIT: Symmetry of OITMAT 523C OITMAT: One-index transformation matrix 524C N2DIS : number of distributions of each symmetry 525C INDAB : index for (AB) symmetry packed integrals. 526C NH2XCD: Number of allocated buffers in H2X 527C 528C Scratch: 529C The rest, including 530C WRK(KFRSAV:KFRSAV-1+LFRSAV) 531C 532C ********************************************************************* 533C 534#include "implicit.h" 535 DIMENSION OITMAT(N2ORBX), H2CD(N2ORBX), H2XCD(N2ORBX) 536 DIMENSION H2X(LH2X,NH2XCD), WRK(*), N2DIS(8) 537 DIMENSION ICDTRA(NORBT,NORBT), NEEDMU(-4:6), INDAB(NORBT,NORBT) 538 DIMENSION INDCD(NORBT,NORBT), IN2CD(NORBT,NORBT) 539 PARAMETER ( D1 = 1.0D0 ) 540C 541C Used from common blocks 542C INFORB : NSYM,N2ORBX,N2ORBT,NNORBX,NORBT,... 543C INFIND : ISMO(*) 544C 545#include "maxash.h" 546#include "maxorb.h" 547#include "priunit.h" 548#include "inforb.h" 549#include "infind.h" 550#include "infpri.h" 551C 552 LOGICAL CDEQDC 553 554 CALL QENTER('OIT2M2') 555 KFREE = KFRSAV 556 LFREE = LFRSAV 557C 558C Allocate work memory 559C 560C 561 IF (ICTOIT .LT. 1 .OR. ICTOIT .GT. 3) IPROIT = 5 562 IF (IPROIT .GE. 5) THEN 563 WRITE (LUPRI,'(/A)') 564 & ' Test output from OIT2M2.' 565 WRITE (LUPRI,*) 'LUDAXN :',LUDAXN 566 WRITE (LUPRI,*) 'JTRLVL :',JTRLVL 567 WRITE (LUPRI,*) 'ICTOIT :',ICTOIT 568 WRITE (LUPRI,*) 'LSYOIT :',LSYOIT 569 WRITE (LUPRI,*) 'NH2XCD :',NH2XCD 570 WRITE (LUPRI,*) 'LH2X :',LH2X 571 WRITE (LUPRI,*) 'N2ORBT :',N2ORBT 572 WRITE (LUPRI,*) 'NNORBX :',NNORBX 573 WRITE (LUPRI,*) 'N2DIS :',(N2DIS(I),I=1,NSYM) 574 CALL FLSHFO(LUPRI) 575 END IF 576 IF (ICTOIT .LT. 1 .OR. ICTOIT .GT. 3) THEN 577 WRITE (LUPRI,*) ' *** ERROR: Illegal ICTOIT' 578 CALL QTRACE(LUPRI) 579 CALL QUIT('*** ERROR: Illegal ICTOIT in OIT2M2.') 580 END IF 581C 582C **************************************************************** 583C Loop over Mulliken distributions allowed in NEEDMU(*) 584C 585 CALL IZERO(INDCD,N2ORBX) 586 IF (ICTOIT .EQ. 2) THEN 587 CALL DZERO(H2X,NH2XCD*N2ORBT) 588 END IF 589 JDIST = 0 590 IDIST = 0 591 100 CALL NXTH2M(IC,ID,H2CD,NEEDMU,WRK,KFREE,LFREE,IDIST) 592 IF (IDIST .LT. 0) GO TO 800 593C ... if idist .lt. 0 then no more distributions 594 JDIST = JDIST + 1 595 IF (ICTOIT .EQ. 1 .AND. JDIST .GT. NH2XCD) THEN 596 WRITE (LUPRI,*) 597 * 'OIT2M2 error, insufficient allocation for H2X' 598 WRITE (LUPRI,*) ' -- Allocated :',NH2XCD 599 CALL QTRACE(LUPRI) 600 CALL QUIT('OIT2M2 error, insufficient allocation for H2X') 601 END IF 602 INDCD(IC,ID) = JDIST 603 INDCD(ID,IC) = JDIST 604 IF (IC .LT. ID) THEN 605 IA = IC 606 IC = ID 607 ID = IA 608 END IF 609C 610 ICSYM = ISMO(IC) 611 IDSYM = ISMO(ID) 612 ICDSYM = MULD2H(ICSYM,IDSYM) 613 IABSYM = ICDSYM 614C 615 IF (IPROIT .GE. 40) THEN 616 WRITE (LUPRI,'(/A,I5,A,2I5)') 617 & ' Mulliken distribution no.',JDIST,', IC and ID:',IC,ID 618 WRITE (LUPRI,*) ' IDIST =',IDIST 619 WRITE (LUPRI,*) 'ICDSYM =',ICDSYM 620 CALL OUTPUT(H2CD,1,NORBT,1,NORBT,NORBT,NORBT,-1,LUPRI) 621 END IF 622C 623 CALL DZERO(H2XCD,N2ORBX) 624 CALL OITH1(LSYOIT,OITMAT,H2CD,H2XCD,IABSYM) 625C CALL OITH1(LSYOIT,OITMAT,H1,H1X,IH1SYM) 626C 627 IF (IPROIT .GE. 50) THEN 628 WRITE (LUPRI,'(/A,I5,A,2I5)') 629 & ' OIT2M2 distribution no.',JDIST, 630 & ', IC and ID:',IC,ID 631 WRITE (LUPRI,*) 'One-index transformed in IA and IB:' 632 WRITE (LUPRI,*) 'IABSYM =',IABSYM 633 WRITE (LUPRI,*) 'LSYOIT =',LSYOIT 634 CALL OUTPUT(H2XCD,1,NORBT,1,NORBT,NORBT,NORBT,-1,LUPRI) 635 END IF 636C 637C We now have the one-index transformed integrals in H2XCD(*,*): 638C (X Y / C D) = (X B / C D) + (A Y / C D) 639C where X = A~1 and Y = B~2 640C 641 IXYSYM = MULD2H(IABSYM,LSYOIT) 642 N2DXY = N2DIS(IXYSYM) 643 IF (ICTOIT .EQ. 1) THEN 644C 645C Save half-transformed integrals in H2X 646C 647 CALL OITPAK(H2XCD,H2X(1,JDIST),IXYSYM,1) 648C CALL OITPAK(H2XCD,H2X,IXYSYM,IWAY) 649 ELSE IF (ICTOIT .EQ. 2) THEN 650C 651C Distribute half-transformed integrals in H2X 652C using H2CD as temporary storage for packed integrals 653C 654 CALL OITPAK(H2XCD,H2CD,IXYSYM,1) 655 IRECCD = ICDTRA(IC,ID) 656 IRECDC = ICDTRA(ID,IC) 657 ICDX = INDAB (IC,ID) 658 IDCX = INDAB (ID,IC) 659 IF (IRECCD .GT. 0) THEN 660 CALL DAXPY(N2DXY,D1,H2CD,1,H2X(1,IRECCD),1) 661 END IF 662 IF (IRECDC .GT. 0 .AND. IRECDC .NE. IRECCD) THEN 663 CALL DAXPY(N2DXY,D1,H2CD,1,H2X(1,IRECDC),1) 664 END IF 665 DO 1200 IY = 1,NORBT 666 ISYMY = ISMO(IY) 667 ISYMX = MULD2H(IXYSYM,ISYMY) 668 IXST = IORB(ISYMX) + 1 669 IXEND = IORB(ISYMX) + NORB(ISYMX) 670 DO 1100 IX = IXST,IXEND 671 IRECXY = ICDTRA(IX,IY) 672 IF (IRECXY .GT. 0) THEN 673 IXYX = INDAB(IX,IY) 674 H2X(ICDX,IRECXY) = H2X(ICDX,IRECXY) + H2CD(IXYX) 675 IF (ICDX .NE. IDCX) H2X(IDCX,IRECXY) 676 * = H2X(IDCX,IRECXY) + H2CD(IXYX) 677 END IF 678 1100 CONTINUE 679 1200 CONTINUE 680 ELSE 681C 682C ICTOIT = 3: out of core 683C Save half-transformed integrals on disk 684C 685 CALL OITPAK(H2XCD,H2CD,IXYSYM,1) 686 WRITE (LUDAXT,REC=JDIST) (H2CD(I),I=1,N2DXY) 687 END IF 688C 689C Go to 100 to get next needed Mulliken distribution 690C 691 GO TO 100 692C 693C arrive at 800 when finished with all needed Mulliken distributions 694C 695 800 CONTINUE 696 IF (IPROIT .GE. 5) THEN 697 WRITE (LUPRI,'(//A/,3(/A,I5))') 698 & ' End of test output of Mulliken distributions from OIT2M2.', 699 & ' Total number of distributions treated :',JDIST, 700 & ' Total number of distributions (NNORBX) :',NNORBX, 701 & ' Total number allocated in H2X :',NH2XCD 702 CALL FLSHFO(LUPRI) 703 END IF 704 IF (KFREE .NE. KFRSAV) 705 & CALL MEMREL('OIT2M2-1',WRK,KFRSAV,KFRSAV,KFREE,LFREE) 706C 707C 708C **************************************************************** 709C 710C Now H2X(ij,kl) = (it jt / k l), 711C finish H2X by symmetrizing 712C (remember (it jt / kt lt) = (it jt / k l) + (kt lt / i j) ). 713C 714C If requested, print H2X 715C 716C 717 IF (ICTOIT .EQ. 1 .OR. ICTOIT .EQ. 3) THEN 718 CALL ICOPY(N2ORBX,INDCD,1,IN2CD,1) 719 END IF 720 ISYH2X = LSYOIT 721 CDEQDC = .TRUE. 722 DO 2800 ISYMCD = 1,NSYM 723 ISYMAB = MULD2H(ISYH2X,ISYMCD) 724 N2DAB = N2DIS(ISYMAB) 725 IF (N2DAB .EQ. 0) GO TO 2800 726 IDEND = 0 727 2000 CONTINUE 728 IDST = IDEND + 1 729 IF (ICTOIT .EQ. 3) THEN 730 CALL OITCOR(ISYMAB,ISYMCD,IDST,IDEND,ICDXOF,CDEQDC, 731 * INDCD,IN2CD,INDAB,LUDAXT,H2XCD,H2X,LH2X,NH2XCD) 732 ELSE 733 ICDXOF = 0 734 IDEND = NORBT 735 END IF 736 DO 2400 ID = IDST,IDEND 737 ISYMD = ISMO(ID) 738 ISYMC = MULD2H(ISYMD,ISYMCD) 739 ICST = IORB(ISYMC) + 1 740 ICEND = IORB(ISYMC) + NORB(ISYMC) 741 DO 2300 IC = ICST,ICEND 742 IRECCD = ICDTRA(IC,ID) 743 IF (IRECCD .EQ. 0) GO TO 2300 744 IF (ICTOIT .EQ. 1 .OR. ICTOIT .EQ. 3) THEN 745 ICDX = INDAB(IC,ID) 746 ICDDIS = INDCD(IC,ID) 747 IF (ICDDIS .EQ. 0) THEN 748 CALL QUIT('ERROR OIT2M2,INDCD .eq. 0 when ICDTRA .ne. 0') 749 END IF 750 IF (ICTOIT .EQ. 1) THEN 751 CALL DCOPY(N2DAB,H2X(1,ICDDIS),1,H2XCD,1) 752 ELSE 753 IF (ICDX .LE. ICDXOF) THEN 754 CALL QUIT('ERROR OIT2M2, ICDX .le. '// 755 & 'ICDXOF for ICTOIT = 3') 756 END IF 757 READ (LUDAXT,REC=ICDDIS) (H2XCD(I),I=1,N2DAB) 758 END IF 759 DO 2200 IB = 1,NORBT 760 ISYMB = ISMO(IB) 761 ISYMA = MULD2H(ISYMB,ISYMAB) 762 IAST = IORB(ISYMA) + 1 763 IAEND = IORB(ISYMA) + NORB(ISYMA) 764 DO 2100 IA = IAST,IAEND 765 IABDIS = IN2CD(IA,IB) 766 IF (IABDIS .EQ. 0) GO TO 2100 767 IABX = INDAB(IA,IB) 768 IF (IABDIS .GT. 0) THEN 769 H2XCD(IABX) = H2XCD(IABX) + H2X(ICDX-ICDXOF,IABDIS) 770 ELSE 771 READ(LUDAXT,REC=-IABDIS) (H2CD(I),I=1,ICDX) 772 H2XCD(IABX) = H2XCD(IABX) + H2CD(ICDX) 773 END IF 774 2100 CONTINUE 775 2200 CONTINUE 776 END IF 777 IF (IPROIT .GE. 40) THEN 778 WRITE (LUPRI,'(/A,I8,A,2I5)') 779 & ' OIT2M2 record no.',IRECCD,', IC and ID:',IC,ID 780 WRITE (LUPRI,*) 'ISYMCD =',ISYMCD 781 WRITE (LUPRI,*) 'ISYMAB =',ISYMAB 782 WRITE (LUPRI,*) 'LSYOIT =',LSYOIT 783 IF (ICTOIT .EQ. 1 .OR. ICTOIT .EQ. 3) THEN 784 CALL OITPAK(H2CD,H2XCD,ISYMAB,-1) 785 ELSE 786 CALL OITPAK(H2CD,H2X(1,IRECCD),ISYMAB,-1) 787 END IF 788 CALL OUTPUT(H2CD,1,NORBT,1,NORBT,NORBT,NORBT,-1,LUPRI) 789 END IF 790 IF (ICTOIT .EQ. 1 .OR. ICTOIT .EQ. 3) THEN 791 WRITE (LUDAXN,REC=IRECCD) (H2XCD(I),I=1,N2DAB) 792 ELSE 793 WRITE (LUDAXN,REC=IRECCD) (H2X(I,IRECCD),I=1,N2DAB) 794 END IF 795 2300 CONTINUE 796 2400 CONTINUE 797 IF (IDEND .LT. NORBT) GO TO 2000 798 2800 CONTINUE 799C 800 IF (KFREE .NE. KFRSAV) 801 & CALL MEMREL('OIT2M2-2',WRK,KFRSAV,KFRSAV,KFREE,LFREE) 802C 803C *** end of subroutine OIT2M2 804C 805 9999 CALL QEXIT('OIT2M2') 806 RETURN 807 END 808C /* Deck oit2x */ 809 SUBROUTINE OIT2X (IDAXO,LUDAXN,JTRLVL,LSYMXO,LSYOIT,OITMAT, 810 * WRK,KFRSAV,LFRSAV,IPROIT) 811C 812C Copyright 21. Nov 1990 by Hans Jorgen Aa. Jensen. 813C 814C Purpose: 815C Construct one-index transformed 2-electron integrals. 816C 817C Input: 818C IDAXO : Identification of input direct access file 819C LUDAXN: unit number for output direct access file. 820C JTRLVL: Transformation level of transformed integrals. 821C LSYMXO: Symmetry of "old" integrals 822C LSYOIT: Symmetry of OITMAT 823C OITMAT: One-index transformation matrix 824C 825C Scratch: 826C WRK(KFRSAV:KFRSAV-1+LFRSAV) 827C 828C ********************************************************************* 829C 830#include "implicit.h" 831 DIMENSION OITMAT(N2ORBX), WRK(LFRSAV) 832C 833C Used from common blocks 834C INFORB : NSYM,N2ORBX,N2ORBT,NNORBX,NORBT,... 835C INFDIM : MWORK, NORBMA, ... 836C 837#include "maxorb.h" 838#include "inforb.h" 839#include "infdim.h" 840#include "infpri.h" 841C 842C Local variables 843C 844 DIMENSION NEEDMU(-4:6), N2DIS(8) 845 846 CALL QENTER('OIT2X') 847 KFREE = KFRSAV 848 LFREE = LFRSAV 849C 850C Set NEEDMU array 851C ... occupied-occupied distributions always needed 852C ... unocc-occ distrib. needed for level 1 and higher 853C ... all distributions are needed for level 2 and higher 854C 855 NEEDMU(-4:6) = 0 856 IF (JTRLVL .GE. 2) THEN 857 NEEDMU(1:6) = 1 858 ELSE IF (JTRLVL .GE. 1) THEN 859 NEEDMU(1:5) = 1 860 ELSE 861 NEEDMU(1:3) = 1 862 END IF 863C 864C Determine if in core or out of core 865C 866 MWORK1 = MIN(MWORK,LFREE-10*N2ORBX) 867C ... 10 is an arbitrary number, which is hoped to be sufficient 868 NH2XCD = MWORK1 / N2ORBT 869 IF (JTRLVL .EQ. 0 .OR. JTRLVL .EQ. 1) THEN 870 IF (NH2XCD .GE. N2OCCX) THEN 871 ICTOIT = 2 872 NH2XCD = N2OCCX 873 ELSE 874 ICTOIT = 3 875 IF (JTRLVL .EQ. 0) THEN 876 NH2XCD = NSYM*(NISHMA+NASHMA)*(NISHMA+NASHMA) 877 NH2XCD = MIN(NH2XCD,N2OCCX) 878C MAERKE need N2OCCT but that is not defined 879 ELSE 880 NH2XCD = NSYM*(NISHMA + NASHMA)*(NSSHMA + NORBMA) 881 NH2XCD = MIN(NH2XCD,N2OCCX + 2*NOCCT*NSSHT) 882C MAERKE need max of all symmetries but that is not defined 883 END IF 884 END IF 885 ELSE IF (JTRLVL .EQ. 2) THEN 886 IF (NH2XCD .GE. N2OCCX + 2*NOCCT*(NORBT-NOCCT)) THEN 887 ICTOIT = 2 888 NH2XCD = N2OCCX + 2*NOCCT*(NORBT-NOCCT) 889 ELSE 890 ICTOIT = 3 891 NH2XCD = N2ORBT 892 END IF 893 ELSE 894 IF (NH2XCD .GE. N2ORBX) THEN 895 ICTOIT = 1 896 NH2XCD = N2ORBX 897 ELSE 898 ICTOIT = 3 899 NH2XCD = N2ORBT 900 END IF 901 END IF 902C 903C Allocate work memory 904C 905 IF (ICTOIT .EQ. 3) THEN 906 LH2X = MIN(N2ORBT,MWORK1/NH2XCD) 907 LH2X = MAX(NORBMA,LH2X) 908 ELSE 909 LH2X = N2ORBT 910 END IF 911 LH2XT = LH2X*NH2XCD 912 CALL MEMGET('REAL',KH2CD ,N2ORBX,WRK,KFREE,LFREE) 913 CALL MEMGET('REAL',KH2XCD,N2ORBX,WRK,KFREE,LFREE) 914 CALL MEMGET('INTE',KINDAB,N2ORBX,WRK,KFREE,LFREE) 915 CALL MEMGET('INTE',KINDCD,N2ORBX,WRK,KFREE,LFREE) 916 CALL MEMGET('INTE',KIN2CD,N2ORBX,WRK,KFREE,LFREE) 917 CALL MEMGET('INTE',KICDTR,N2ORBX,WRK,KFREE,LFREE) 918 CALL MEMGET('REAL',KH2X ,LH2XT ,WRK,KFREE,LFREE) 919C 920C Open auxiliary DA file 921C 922 IF (ICTOIT .EQ. 3) THEN 923 IDAXX = -1 924 IDAXT = 999 925 JTRLVT= 4 926 LRDAX = N2ORBT 927 CALL OITOPN(IDAXX,IDAXT,LUDAXX,LUDAXT,JTRLVX,JTRLVT, 928 * LSYMXO,LSYOIT,LRDAX) 929C CALL OITOPN(IDAXO,IDAXN,LUDAXO,LUDAXN,JTRLVO,JTRLVN, 930C * LSYMXO,LSYOIT,LRDAX) 931 ELSE 932 LUDAXT = -999 933 END IF 934C 935C Set INDAB and ICDTRA arrays 936C 937 CALL OITIND(WRK(KINDAB),N2DIS) 938 CALL OITICD(JTRLVL,WRK(KICDTR),WRK(KINDCD),IPROIT) 939C CALL OITICD(JTRLVL,ICDTRA,ITRTYP,IPROIT) 940 CALL OIT2X2(IDAXO,LUDAXN,JTRLVL,ICTOIT,LUDAXT,LSYMXO,LSYOIT, 941 * OITMAT,WRK(KH2CD),WRK(KH2XCD), WRK(KICDTR), 942 * NEEDMU,N2DIS,WRK(KINDAB),WRK(KINDCD),WRK(KIN2CD), 943 * WRK(KH2X),LH2X,NH2XCD, WRK,KFREE,LFREE,IPROIT) 944C 945 IF (ICTOIT .EQ. 3) CALL OITCLO(IDAXT,'DELETE') 946C 947C *** end of subroutine OIT2X 948C 949 CALL MEMREL('OIT2X',WRK,KFRSAV,KFRSAV,KFREE,LFREE) 950 CALL QEXIT('OIT2X') 951 RETURN 952 END 953C /* Deck oit2x2 */ 954 SUBROUTINE OIT2X2(IDAXO,LUDAXN,JTRLVL,ICTOIT,LUDAXT,LSYMXO,LSYOIT, 955 * OITMAT,H2CD,H2XCD, 956 * ICDTRA,NEEDMU,N2DIS,INDAB,INDCD,IN2CD, 957 * H2X,LH2X,NH2XCD,WRK,KFRSAV,LFRSAV,IPROIT) 958C 959C Copyright 15. Nov 1990 by Hans Jorgen Aa. Jensen. 960C 961C Purpose: 962C Construct one-index transformed 2-electron integrals. 963C 964C Input: 965C IDAXO : Identification of "old" integrals 966C LUDAXN: unit number for output direct access file. 967C JTRLVL: Transformation level of transformed integrals. 968C LSYMXO: Symmetry of "old" integrals 969C LSYOIT: Symmetry of OITMAT 970C ICTOIT: Control parameter for transformation 971C LUDAXT: temporary file for out-of-core transformation, if ICTOIT=3 972C OITMAT: One-index transformation matrix 973C N2DIS : number of distributions of each symmetry 974C INDAB : index for (AB) symmetry packed integrals. 975C NH2XCD: Number of allocated buffers in H2X 976C 977C Scratch: 978C The rest, including 979C WRK(KFRSAV:KFRSAV-1+LFRSAV) 980C 981C ********************************************************************* 982C 983#include "implicit.h" 984 DIMENSION OITMAT(N2ORBX), H2CD(N2ORBX), H2XCD(N2ORBX) 985 DIMENSION H2X(LH2X,NH2XCD), WRK(*), N2DIS(8) 986 DIMENSION ICDTRA(NORBT,NORBT), NEEDMU(-4:6), INDAB(NORBT,NORBT) 987 DIMENSION INDCD(NORBT,NORBT), IN2CD(NORBT,NORBT) 988 PARAMETER ( D1 = 1.0D0 ) 989C 990C Used from common blocks 991C INFORB : NSYM,N2ORBX,N2ORBT,NNORBX,NORBT,... 992C INFIND : ISMO(*) 993C 994#include "maxash.h" 995#include "maxorb.h" 996#include "priunit.h" 997#include "inforb.h" 998#include "infind.h" 999#include "infpri.h" 1000 1001 LOGICAL CDEQDC 1002 1003 CALL QENTER('OIT2X2') 1004 KFREE = KFRSAV 1005 LFREE = LFRSAV 1006C 1007C Allocate work memory 1008C 1009C 1010 IF (ICTOIT .LT. 1 .OR. ICTOIT .GT. 3) IPROIT = 5 1011 IF (IPROIT .GE. 5) THEN 1012 WRITE (LUPRI,'(///A)') 1013 & ' Test output from OIT2X2.' 1014 WRITE (LUPRI,*) ' IDAXO :',IDAXO 1015 WRITE (LUPRI,*) 'LUDAXN :',LUDAXN 1016 WRITE (LUPRI,*) 'JTRLVL :',JTRLVL 1017 WRITE (LUPRI,*) 'LSYMXO :',LSYMXO 1018 WRITE (LUPRI,*) 'LSYOIT :',LSYOIT 1019 WRITE (LUPRI,*) 'ICTOIT :',ICTOIT 1020 WRITE (LUPRI,*) 'NH2XCD :',NH2XCD 1021 WRITE (LUPRI,*) 'LH2X :',LH2X 1022 WRITE (LUPRI,*) 'N2ORBT :',N2ORBT 1023 WRITE (LUPRI,*) 'N2OCCX :',N2OCCX 1024 WRITE (LUPRI,*) 'N2ORBX :',N2ORBX 1025 WRITE (LUPRI,*) 'N2DIS :',(N2DIS(I),I=1,NSYM) 1026 CALL FLSHFO(LUPRI) 1027 END IF 1028 IF (ICTOIT .LT. 1 .OR. ICTOIT .GT. 3) THEN 1029 WRITE (LUPRI,*) ' *** ERROR: Illegal ICTOIT' 1030 CALL QTRACE(LUPRI) 1031 CALL QUIT('*** ERROR: Illegal ICTOIT in OIT2X2') 1032 END IF 1033C 1034C **************************************************************** 1035C Loop over Mulliken distributions allowed in NEEDMU(*) 1036C 1037 CALL IZERO(INDCD,N2ORBX) 1038 IF (ICTOIT .EQ. 2) THEN 1039 CALL DZERO(H2X,NH2XCD*N2ORBT) 1040 END IF 1041 JDIST = 0 1042 IDIST = 0 1043 100 CALL NXTH2X(IC,ID,H2CD,IDAXO,NEEDMU,WRK,KFREE,LFREE,IDIST) 1044 IF (IDIST .LT. 0) GO TO 800 1045C ... if idist .lt. 0 then no more distributions 1046 JDIST = JDIST + 1 1047 IF (ICTOIT .EQ. 1 .AND. JDIST .GT. NH2XCD) THEN 1048 WRITE (LUPRI,*) 1049 * 'OIT2X2 error, insufficient allocation for H2X' 1050 WRITE (LUPRI,*) ' -- Allocated :',NH2XCD 1051 CALL QTRACE(LUPRI) 1052 CALL QUIT('OIT2X2 error, insufficient allocation for H2X') 1053 END IF 1054 INDCD(IC,ID) = JDIST 1055C 1056 ICSYM = ISMO(IC) 1057 IDSYM = ISMO(ID) 1058 ICDSYM = MULD2H(ICSYM,IDSYM) 1059 IABSYM = MULD2H(LSYMXO,ICDSYM) 1060C 1061 IF (IPROIT .GE. 40) THEN 1062 WRITE (LUPRI,'(/A,I5,A,2I5)') 1063 & ' OIT2X2 Distribution no.',JDIST,', IC and ID:',IC,ID 1064 WRITE (LUPRI,*) ' IDIST =',IDIST 1065 WRITE (LUPRI,*) 'ICDSYM =',ICDSYM 1066 WRITE (LUPRI,*) 'IABSYM =',IABSYM 1067 CALL OUTPUT(H2CD,1,NORBT,1,NORBT,NORBT,NORBT,-1,LUPRI) 1068 END IF 1069C 1070 CALL DZERO(H2XCD,N2ORBX) 1071 CALL OITH1(LSYOIT,OITMAT,H2CD,H2XCD,IABSYM) 1072C CALL OITH1(LSYOIT,OITMAT,H1,H1X,IH1SYM) 1073C 1074 IF (IPROIT .GE. 50) THEN 1075 WRITE (LUPRI,'(/A,I5,A,2I5)') 1076 & ' OIT2X2 distribution no.',JDIST,', IC and ID:',IC,ID 1077 WRITE (LUPRI,*) 'One-index transformed in IA and IB:' 1078 WRITE (LUPRI,*) 'IABSYM =',IABSYM 1079 WRITE (LUPRI,*) 'LSYOIT =',LSYOIT 1080 CALL OUTPUT(H2XCD,1,NORBT,1,NORBT,NORBT,NORBT,-1,LUPRI) 1081 END IF 1082C 1083C We now have the one-index transformed integrals in H2XCD(*,*): 1084C (X Y / C D) = (X B / C D) + (A Y / C D) 1085C where X = A~1 and Y = B~2 1086C 1087 IXYSYM = MULD2H(IABSYM,LSYOIT) 1088 N2DXY = N2DIS(IXYSYM) 1089 IF (ICTOIT .EQ. 1) THEN 1090C 1091C Save half-transformed integrals in H2X 1092C 1093 CALL OITPAK(H2XCD,H2X(1,JDIST),IXYSYM,1) 1094C CALL OITPAK(H2XCD,H2X,IXYSYM,IWAY) 1095 ELSE IF (ICTOIT .EQ. 2) THEN 1096C 1097C Distribute half-transformed integrals in H2X 1098C using H2CD as temporary storage for packed integrals 1099C 1100 CALL OITPAK(H2XCD,H2CD,IXYSYM,1) 1101 IRECCD = ICDTRA(IC,ID) 1102 ICDX = INDAB (IC,ID) 1103 IF (IRECCD .GT. 0) THEN 1104 CALL DAXPY(N2DXY,D1,H2CD,1,H2X(1,IRECCD),1) 1105 END IF 1106 DO 1200 IY = 1,NORBT 1107 ISYMY = ISMO(IY) 1108 ISYMX = MULD2H(IXYSYM,ISYMY) 1109 IXST = IORB(ISYMX) + 1 1110 IXEND = IORB(ISYMX) + NORB(ISYMX) 1111 DO 1100 IX = IXST,IXEND 1112 IRECXY = ICDTRA(IX,IY) 1113 IF (IRECXY .GT. 0) THEN 1114 IXYX = INDAB(IX,IY) 1115 H2X(ICDX,IRECXY) = H2X(ICDX,IRECXY) + H2CD(IXYX) 1116 END IF 1117 1100 CONTINUE 1118 1200 CONTINUE 1119 ELSE 1120C 1121C ICTOIT = 3: out of core 1122C Save half-transformed integrals on disk 1123C 1124 CALL OITPAK(H2XCD,H2CD,IXYSYM,1) 1125 WRITE (LUDAXT,REC=JDIST) (H2CD(I),I=1,N2DXY) 1126 END IF 1127C 1128C Go to 100 to get next needed Mulliken distribution 1129C 1130 GO TO 100 1131C 1132C arrive at 800 when finished with all needed Mulliken distributions 1133C 1134 800 CONTINUE 1135 IF (IPROIT .GE. 5) THEN 1136 WRITE (LUPRI,'(//A/,3(/A,I5))') 1137 & ' End of test output of input distributions from OIT2X2.', 1138 & ' Total number of distributions treated :',JDIST, 1139 & ' Total number of distributions (N2ORBX) :',N2ORBX, 1140 & ' Total number allocated in H2X :',NH2XCD 1141 CALL FLSHFO(LUPRI) 1142 END IF 1143 IF (KFREE .NE. KFRSAV) 1144 & CALL MEMREL('OIT2X2-1',WRK,KFRSAV,KFRSAV,KFREE,LFREE) 1145C 1146C 1147C **************************************************************** 1148C 1149C Now H2X(ij,kl) = (it jt / k l), 1150C finish H2X by symmetrizing 1151C (remember (it jt / kt lt) = (it jt / k l) + (kt lt / i j) ). 1152C 1153C If requested, print H2X 1154C 1155C 1156 IF (ICTOIT .EQ. 1 .OR. ICTOIT .EQ. 3) THEN 1157 CALL ICOPY(N2ORBX,INDCD,1,IN2CD,1) 1158 END IF 1159 ISYH2X = MULD2H(LSYOIT,LSYMXO) 1160 CDEQDC = .FALSE. 1161 DO 2800 ISYMCD = 1,NSYM 1162 ISYMAB = MULD2H(ISYH2X,ISYMCD) 1163 N2DAB = N2DIS(ISYMAB) 1164 IF (N2DAB .EQ. 0) GO TO 2800 1165 IDEND = 0 1166 2000 CONTINUE 1167 IDST = IDEND + 1 1168 IF (ICTOIT .EQ. 3) THEN 1169 CALL OITCOR(ISYMAB,ISYMCD,IDST,IDEND,ICDXOF,CDEQDC, 1170 * INDCD,IN2CD,INDAB,LUDAXT,H2XCD,H2X,LH2X,NH2XCD) 1171 ELSE 1172 ICDXOF = 0 1173 IDEND = NORBT 1174 END IF 1175 DO 2400 ID = IDST,IDEND 1176 ISYMD = ISMO(ID) 1177 ISYMC = MULD2H(ISYMD,ISYMCD) 1178 ICST = IORB(ISYMC) + 1 1179 ICEND = IORB(ISYMC) + NORB(ISYMC) 1180 DO 2300 IC = ICST,ICEND 1181 IRECCD = ICDTRA(IC,ID) 1182 IF (IRECCD .EQ. 0) GO TO 2300 1183 IF (ICTOIT .EQ. 1 .OR. ICTOIT .EQ. 3) THEN 1184 ICDX = INDAB(IC,ID) 1185 ICDDIS = INDCD(IC,ID) 1186 IF (ICDDIS .EQ. 0) THEN 1187 CALL QUIT('ERROR OIT2X2, INDCD .eq. 0 when ICDTRA .ne.0') 1188 END IF 1189 IF (ICTOIT .EQ. 1) THEN 1190 CALL DCOPY(N2DAB,H2X(1,ICDDIS),1,H2XCD,1) 1191 ELSE 1192 IF (ICDX .LE. ICDXOF) THEN 1193 CALL QUIT('ERROR OIT2X2, ICDX .le. '// 1194 & 'ICDXOF for ICTOIT = 3') 1195 END IF 1196 READ (LUDAXT,REC=ICDDIS) (H2XCD(I),I=1,N2DAB) 1197 END IF 1198 DO 2200 IB = 1,NORBT 1199 ISYMB = ISMO(IB) 1200 ISYMA = MULD2H(ISYMB,ISYMAB) 1201 IAST = IORB(ISYMA) + 1 1202 IAEND = IORB(ISYMA) + NORB(ISYMA) 1203 DO 2100 IA = IAST,IAEND 1204 IABDIS = IN2CD(IA,IB) 1205 IF (IABDIS .EQ. 0) GO TO 2100 1206 IABX = INDAB(IA,IB) 1207 IF (IABDIS .GT. 0) THEN 1208 H2XCD(IABX) = H2XCD(IABX) + H2X(ICDX-ICDXOF,IABDIS) 1209 ELSE 1210 READ(LUDAXT,REC=-IABDIS) (H2CD(I),I=1,ICDX) 1211 H2XCD(IABX) = H2XCD(IABX) + H2CD(ICDX) 1212 END IF 1213 2100 CONTINUE 1214 2200 CONTINUE 1215 END IF 1216 IF (IPROIT .GE. 40) THEN 1217 WRITE (LUPRI,'(/A,I8,A,2I5)') 1218 & ' OIT2X2 record no.',IRECCD,', IC and ID:',IC,ID 1219 WRITE (LUPRI,*) 'ISYMCD =',ISYMCD 1220 WRITE (LUPRI,*) 'ISYMAB =',ISYMAB 1221 WRITE (LUPRI,*) 'LSYOIT =',LSYOIT 1222 IF (ICTOIT .EQ. 1 .OR. ICTOIT .EQ. 3) THEN 1223 CALL OITPAK(H2CD,H2XCD,ISYMAB,-1) 1224 ELSE 1225 CALL OITPAK(H2CD,H2X(1,IRECCD),ISYMAB,-1) 1226 END IF 1227 CALL OUTPUT(H2CD,1,NORBT,1,NORBT,NORBT,NORBT,-1,LUPRI) 1228 END IF 1229 IF (ICTOIT .EQ. 1 .OR. ICTOIT .EQ. 3) THEN 1230 WRITE (LUDAXN,REC=IRECCD) (H2XCD(I),I=1,N2DAB) 1231 ELSE 1232 WRITE (LUDAXN,REC=IRECCD) (H2X(I,IRECCD),I=1,N2DAB) 1233 END IF 1234 2300 CONTINUE 1235 2400 CONTINUE 1236 IF (IDEND .LT. NORBT) GO TO 2000 1237 2800 CONTINUE 1238C 1239 IF (KFREE .NE. KFRSAV) 1240 & CALL MEMREL('OIT2X2-2',WRK,KFRSAV,KFRSAV,KFREE,LFREE) 1241C 1242C *** end of subroutine OIT2X2 1243C 1244 9999 CALL QEXIT('OIT2X2') 1245 RETURN 1246 END 1247C /* Deck oitpak */ 1248 SUBROUTINE OITPAK(H2XCD,H2X,IXYSYM,IWAY) 1249C 1250C Copyright 16-Nov-1990 Hans Joergen Aa. Jensen 1251C 1252C IWAY .ge. 0 : Pack H2XCD in H2X 1253C else : Unpack H2X in H2XCD 1254C 1255#include "implicit.h" 1256 DIMENSION H2XCD(NORBT,NORBT), H2X(N2ORBT) 1257C 1258C Used from common blocks: 1259C INFORB : NORBT,N2ORBT,N2ORBX,NSYM,MULD2H 1260C 1261#include "inforb.h" 1262C 1263 IF (IWAY .LT. 0) CALL DZERO(H2XCD,N2ORBX) 1264 ISTH2X = 1 1265 DO 300 IBSYM = 1,NSYM 1266 IASYM = MULD2H(IBSYM,IXYSYM) 1267 NORBA = NORB(IASYM) 1268 NORBB = NORB(IBSYM) 1269 IAST = IORB(IASYM) + 1 1270 IBST = IORB(IBSYM) + 1 1271 IF (IWAY .GE. 0) THEN 1272 CALL MCOPY(NORBA,NORBB,H2XCD(IAST,IBST),NORBT, 1273 * H2X(ISTH2X),NORBA) 1274 ELSE 1275 CALL MCOPY(NORBA,NORBB,H2X(ISTH2X),NORBA, 1276 * H2XCD(IAST,IBST),NORBT) 1277 END IF 1278 ISTH2X = ISTH2X + NORBA*NORBB 1279 300 CONTINUE 1280C 1281C MCOPY(nrowa,ncola,A,nrdima,B,nrdimb) 1282C 1283 RETURN 1284 END 1285C /* Deck oitind */ 1286 SUBROUTINE OITIND(INDAB,N2DIS) 1287C 1288C Copyright 16-Nov-1990 Hans Joergen Aa. Jensen 1289C 1290C Find index for integral (AB) in symmetry packed integral list 1291C 1292#include "implicit.h" 1293 DIMENSION INDAB(NORBT,NORBT), N2DIS(8) 1294C 1295C Used from common blocks: 1296C INFORB : NORBT,N2ORBT,N2ORBX,NSYM,MULD2H 1297C 1298#include "inforb.h" 1299C 1300 LOGICAL NOIND 1301C 1302 NOIND = .FALSE. 1303 GO TO 10 1304 ENTRY OITDIS(N2DIS) 1305 NOIND = .TRUE. 1306 10 CONTINUE 1307C 1308 DO 400 IABSYM = 1,NSYM 1309 INDXAB = 0 1310 DO 300 IBSYM = 1,NSYM 1311 IASYM = MULD2H(IBSYM,IABSYM) 1312 IF (NOIND) THEN 1313 INDXAB = INDXAB + NORB(IASYM)*NORB(IBSYM) 1314 ELSE 1315 IAST = IORB(IASYM) + 1 1316 IAEND = IORB(IASYM) + NORB(IASYM) 1317 IBST = IORB(IBSYM) + 1 1318 IBEND = IORB(IBSYM) + NORB(IBSYM) 1319 DO 200 IB = IBST,IBEND 1320 DO 100 IA = IAST,IAEND 1321 INDXAB = INDXAB + 1 1322 INDAB(IA,IB) = INDXAB 1323 100 CONTINUE 1324 200 CONTINUE 1325 END IF 1326 300 CONTINUE 1327 N2DIS(IABSYM) = INDXAB 1328 400 CONTINUE 1329C 1330 RETURN 1331 END 1332C /* Deck oiticd */ 1333 SUBROUTINE OITICD(JTRLVL,ICDTRA,ITRTYP,IPROIT) 1334C 1335C Copyright 19-Nov-1990 Hans Joergen Aa. Jensen 1336C 1337C Find file index for (C D) distribution of final integrals. 1338C 1339C JTRLVL = general transformation level 1340C KTRLVL = 2 in this version 1341C 1342C KTRLVL = 0 : a designates inactive+active orbitals 1343C g designates general non-frozen orbitals 1344C KTRLVL = 1 : a designates active orbitals 1345C g designates general non-frozen orbitals 1346C KTRLVL = 2 : a designates frozen+inactive+active orbitals 1347C g designates all orbitals 1348C 1349C JTRLVL INDEX 1 INDEX 2 INDEX 3 INDEX 4 comment 1350C 0 a a a a CI calc. (0. ord.) 1351C 1 a a a g first-order 1352C 2 a a g g second-order 1353C a g a g 1354C 3 a g g g third-order 1355C g g a g 1356C 4 g g g g full transf. 1357C and permutations between indices 1 and 2 and between 3 and 4. 1358C 1359C Indices 1 and 2 define the distributions needed (stored in 1360C ICDTRA(1:NORBT,1:NORBT)). 1361C The (g g) distributions are needed at level 3 in order to perform 1362C a one-index transformation leading to (aa/gg) integrals at level 2 1363C 1364C ITRTYP(1:NORBT) = number of integral indices in which this orbital 1365C enters (i.e. 0,1,2,3, or 4) 1366C 1367#include "implicit.h" 1368 DIMENSION ICDTRA(NORBT,NORBT), ITRTYP(NORBT) 1369 PARAMETER (KTRLVL = 2) 1370C 1371C Used from common blocks: 1372C INFORB : NSYM,NORBT,... 1373C INFIND : ISMO(*) 1374C 1375#include "maxash.h" 1376#include "maxorb.h" 1377#include "priunit.h" 1378#include "inforb.h" 1379#include "infind.h" 1380#include "infpri.h" 1381C 1382 CHARACTER*35 EXTPRM(3) 1383 DATA EXTPRM/'inactive + active orbitals ', 1384 & 'only active orbitals ', 1385 & 'frozen + inactive + active orbitals'/ 1386C 1387 CALL QENTER('OITICD') 1388C 1389 IF (KTRLVL .LT. 0 .OR. KTRLVL .GT. 2) THEN 1390 WRITE (LUPRI,*) 'Illegal KTRLVL in OITICD, KTRLVL =',KTRLVL 1391 CALL QTRACE(LUPRI) 1392 CALL QUIT('FATAL ERROR OITICD: Illegal KTRLVL') 1393 END IF 1394 IF (IPROIT .GE. 10) THEN 1395 WRITE (LUPRI,'(/1X,A,I5/1X,2A/)') 1396 & 'Integral transformation order : ',JTRLVL, 1397 & 'Extent of primary space : ',EXTPRM(KTRLVL+1) 1398 END IF 1399C 1400C Calculate ITRTYP. 1401C 1402 CALL IZERO(ITRTYP,NORBT) 1403 DO 300 ISYM = 1, NSYM 1404 IF (KTRLVL .EQ. 0) THEN 1405 IMOG = IORB(ISYM) + NFRO(ISYM) + 1 1406 IMOA = IMOG 1407 NMOA = NISH(ISYM) + NASH(ISYM) 1408 ELSE IF (KTRLVL .EQ. 1) THEN 1409 IMOG = IORB(ISYM) + NFRO(ISYM) + 1 1410 IMOA = IMOG + NISH(ISYM) 1411 NMOA = NASH(ISYM) 1412 ELSE IF (KTRLVL .EQ. 2) THEN 1413 IMOG = IORB(ISYM) + 1 1414 IMOA = IMOG 1415 NMOA = NFRO(ISYM) + NISH(ISYM) + NASH(ISYM) 1416 END IF 1417 IMOL = IORB(ISYM) + NORB(ISYM) 1418C 1419 NG = MIN(4,JTRLVL) 1420 DO 240 I = IMOG ,IMOL 1421 ITRTYP(I) = NG 1422 240 CONTINUE 1423 DO 250 I = IMOA ,IMOA - 1 + NMOA 1424 ITRTYP(I) = 4 1425 250 CONTINUE 1426 300 CONTINUE 1427C 1428C Determine ICDTRA matrix 1429C 1430 CALL IZERO(ICDTRA,N2ORBX) 1431 IDIST = 0 1432 DO 790 IJSYM = 1,NSYM 1433 DO 780 I = 1,NORBT 1434 IF (ITRTYP(I) .GE. 2) THEN 1435 ISYM = ISMO(I) 1436 JSYM = MULD2H(ISYM,IJSYM) 1437 JST = IORB(JSYM) + 1 1438 JEND = IORB(JSYM) + NORB(JSYM) 1439 DO 770 J = JST,JEND 1440 IF ( (ITRTYP(I)+ITRTYP(J)) .GE. 6 ) THEN 1441 IDIST = IDIST + 1 1442 ICDTRA(J,I) = IDIST 1443 END IF 1444 770 CONTINUE 1445 END IF 1446 780 CONTINUE 1447 790 CONTINUE 1448 IF (IPROIT .GE. 25) THEN 1449 WRITE (LUPRI,'(/A)') 1450 & ' ICDTRA: list of distributions included :' 1451 DO 810 I = 1,NORBT 1452 WRITE(LUPRI,'(I5,A,(T8,10I7))') 1453 * I,' :',(ICDTRA(I,J),J=1,NORBT) 1454 810 CONTINUE 1455 END IF 1456 CALL QEXIT ('OITICD') 1457 RETURN 1458 END 1459C /* Deck oitcor */ 1460 SUBROUTINE OITCOR(ISYMAB,ISYMCD,IDST,IDEND,ICDXOF,CDEQDC, 1461 * INDCD,IN2CD,INDAB,LUDAXT,H2XCD,H2X,LH2X,NH2XCD) 1462C 1463C Copyright Hans Joergen Aa. Jensen December 1990 1464C 1465C Purpose: read as many half-transformed integrals of symmetry ISYMAB 1466C into core as possible. 1467C 1468C Input: 1469C CDEQDC: if true then INDCD(i,j) = INDCD(j,i) 1470C IDST : Starting ID 1471C 1472#include "implicit.h" 1473#include "priunit.h" 1474 DIMENSION INDCD(NORBT,NORBT),IN2CD(NORBT,NORBT),INDAB(NORBT,NORBT) 1475 DIMENSION H2XCD(N2ORBX),H2X(LH2X,NH2XCD) 1476 LOGICAL CDEQDC 1477C 1478C Used from common blocks: 1479C INFORB: NORBT, N2ORBT,... 1480C INFIND: ISMO(*) 1481C 1482#include "maxash.h" 1483#include "maxorb.h" 1484#include "inforb.h" 1485#include "infind.h" 1486C 1487 CALL QENTER('OITCOR') 1488C 1489C Copy and change sign, negative sign will signify that 1490C the distribution is only on disk 1491C 1492 DO 2 J = 1,NORBT 1493 DO 2 I = 1,NORBT 1494 IN2CD(I,J) = -INDCD(I,J) 1495 2 CONTINUE 1496C 1497C Find ICDXOF and IDEND 1498C 1499 ISYMD = ISMO(IDST) 1500 100 ISYMC = MULD2H(ISYMD,ISYMCD) 1501 IF (NORB(ISYMC) .EQ. 0 .OR. NORB(ISYMD) .EQ. 0) THEN 1502 ISYMD = ISYMD + 1 1503 IF (ISYMD .GT. NSYM) THEN 1504 ICDXOF = 0 1505 IDST = NORBT + 1 1506 IDEND = NORBT 1507 GO TO 9999 1508 END IF 1509 IDST = IORB(ISYMD) + 1 1510 GO TO 100 1511 END IF 1512 ICST = IORB(ISYMC) + 1 1513 ICDST = INDAB(ICST,IDST) 1514 ICDXOF = ICDST - 1 1515C 1516 IDEND = 0 1517 LNEED = 0 1518 DO 120 ID = IDST, NORBT 1519 ISYMD = ISMO(ID) 1520 ISYMC = MULD2H(ISYMD,ISYMCD) 1521 IF (NORB(ISYMC) .EQ. 0) GO TO 120 1522 LNEED = LNEED + NORB(ISYMC) 1523 IF (LNEED .GT. LH2X) GO TO 130 1524 IDEND = ID 1525 ICEND = IORB(ISYMC) + NORB(ISYMC) 1526 120 CONTINUE 1527 130 CONTINUE 1528 IF (IDEND .EQ. 0) THEN 1529 WRITE(LUERR,*) 'OITCOR ERROR: IDEND .eq. 0' 1530 WRITE(LUERR,*) ' IDST, IDEND :',IDST,IDEND 1531 WRITE(LUERR,*) ' LH2X :',LH2X 1532 CALL QTRACE(LUERR) 1533 CALL QUIT('OITCOR ERROR: IDEND .eq. 0') 1534 END IF 1535 ICDEND = INDAB(ICEND,IDEND) 1536 NCD = ICDEND - ICDST + 1 1537 IF (NCD .GT. LH2X) THEN 1538 WRITE(LUERR,*) 'OITCOR ERROR: NCD .gt. LH2X' 1539 WRITE(LUERR,*) ' ICST, ICEND :',ICST,ICEND 1540 WRITE(LUERR,*) ' IDST, IDEND :',IDST,IDEND 1541 WRITE(LUERR,*) ' ICDST, ICDEND :',ICDST,ICDEND 1542 WRITE(LUERR,*) ' NCD :',NCD 1543 WRITE(LUERR,*) ' LH2X :',LH2X 1544 CALL QTRACE(LUERR) 1545 CALL QUIT('OITCOR ERROR: NCD .gt. LH2X') 1546 END IF 1547C 1548 IDISXY = 0 1549 ISYMXY = ISYMAB 1550C 1551 DO 220 IY = 1,NORBT 1552 ISYMY = ISMO(IY) 1553 ISYMX = MULD2H(ISYMY,ISYMXY) 1554 IXST = IORB(ISYMX) + 1 1555 IXEND = IORB(ISYMX) + NORB(ISYMX) 1556 DO 210 IX = IXST,IXEND 1557 IRECXY = -IN2CD(IX,IY) 1558 IF (IRECXY .LE. 0) GO TO 210 1559 IDISXY = IDISXY + 1 1560 IF (ICDST .EQ. 1) THEN 1561 READ(LUDAXT,REC=IRECXY) (H2X(I,IDISXY),I=1,ICDEND) 1562 ELSE 1563 READ(LUDAXT,REC=IRECXY) (H2XCD(I),I=1,ICDEND) 1564 CALL DCOPY(NCD,H2XCD(ICDST),1,H2X(1,IDISXY),1) 1565 END IF 1566 IN2CD(IX,IY) = IDISXY 1567 IF (CDEQDC) IN2CD(IY,IX) = IDISXY 1568 IF (IDISXY .EQ. NH2XCD) GO TO 9999 1569C .............. EXIT FROM LOOP 1570 210 CONTINUE 1571 220 CONTINUE 1572C 1573 9999 CALL QEXIT('OITCOR') 1574 RETURN 1575 END 1576C /* Deck nxth2x */ 1577 SUBROUTINE NXTH2X(IC,ID,H2CD,IDAX,NEEDTP,WRK,KFREE,LFREE,IDIST) 1578C 1579C Copyright Hans Joergen Aa. Jensen November 1990 1580C 1581C NOTE: The space allocated in WRK must not be touched outside 1582C until all desired distributions have been read. 1583C 1584C Purpose: 1585C Read next Mulliken two-electron integral distribution (**|cd) 1586C where (cd) distribution is needed according to NEEDTP(ITYPCD) 1587C (if needtp(itypcd) .gt. 0 all distributions of that type needed; 1588C if needtp(itypcd) .lt. 0 at least one distribution needed; 1589C if needtp(itypcd) .eq. 0 no distributions of that type needed). 1590C 1591C Usage: 1592C Set IDIST = 0 before first call of NXTH2X. 1593C DO NOT CHANGE IDIST or WRK(KFREE1:KFREE2-1) in calling routine 1594C until last distribution has been read (signalled by IDIST .eq. -1) 1595C Prototype code: 1596C IDIST = 0 1597C define NEEDTP(-4:6) 1598C 100 CALL NXTH2X(IC,ID,H2CD,NEEDTP,WRK,KFREE,LFREE,IDIST) 1599C IF (IDIST .GT. 0) THEN 1600C KW1 = KFREE 1601C LW1 = LFREE 1602C use (**|cd) distribution in H2CD as desired 1603C WRK(KW1:KW1-1+LW1) may be used 1604C GO TO 100 1605C END IF 1606C 1607C 1608#include "implicit.h" 1609#include "priunit.h" 1610 DIMENSION H2CD(*),NEEDTP(-4:6),WRK(LFREE) 1611C 1612C Used from common blocks: 1613C INFORB : N2ORBX 1614C INFPRI : LUERR 1615C 1616#include "inforb.h" 1617#include "infpri.h" 1618C 1619 DIMENSION N2DIS(8) 1620 LOGICAL CDSWAP 1621 SAVE N2DIS, CDSWAP, IC1, ID1, LUDAX, LSYMX 1622 SAVE KICDTR, KH2CDP, KNEXT 1623 DATA KNEXT /-1/ 1624C 1625 CALL QENTER('NXTH2X') 1626C 1627C MAERKE : print level: 1628C 1629 IPRINT = 0 1630C 1631C Read untransformed integrals if IDAX = 0 1632C 1633 IF (IDAX .EQ. 0) THEN 1634 IF (IDIST .EQ. 0) CDSWAP = .FALSE. 1635 IF (CDSWAP) THEN 1636 IC = ID1 1637 ID = IC1 1638 CDSWAP = .FALSE. 1639 ELSE 1640 CALL NXTH2M(IC1,ID1,H2CD,NEEDTP,WRK,KFREE,LFREE,IDIST) 1641 IC = IC1 1642 ID = ID1 1643 CDSWAP = (IC1 .NE. ID1) 1644 END IF 1645 GO TO 9999 1646 END IF 1647C 1648C If first read (IDIST .eq. 0) then 1649C Open IDAX file and get information about unit and level 1650C allocate space for ICDTRA. 1651C recalculate ICDTRA 1652C 1653 IF (IDIST .EQ. 0) THEN 1654 CALL OITOPO(IDAX,LUDAX,JTRLVO,LSYMX,LRDAX) 1655C CALL OITOPO(IDAXO,LUDAXO,JTRLVO,LSYMXO,LRDAX) 1656C 1657 CALL MEMGET('INTE',KICDTR,N2ORBX,WRK,KFREE,LFREE) 1658 CALL MEMGET('REAL',KH2CDP,N2ORBT,WRK,KFREE,LFREE) 1659 KNEXT = KFREE 1660C Recalculate index vector for distributions 1661 CALL OITICD(JTRLVO,WRK(KICDTR),WRK(KFREE),IPRINT) 1662C CALL OITICD(JTRLVL,ICDTRA,ITRTYP,IPROIT) 1663 CALL OITDIS(N2DIS) 1664 ELSE 1665C ... check that work allocation has not been destroyed by 1666C calling routine. 1667 IF (KNEXT.EQ. -1 ) THEN 1668 WRITE (LUERR,*) 1669 & 'NXTH2X error, IDIST must be zero in first call' 1670 WRITE (LUERR,*) 'IDIST =',IDIST 1671 CALL QTRACE(LUERR) 1672 CALL QUIT('NXTH2X error, IDIST must be zero in first call') 1673 END IF 1674 IF (KFREE.LT.KNEXT) THEN 1675 WRITE (LUERR,*) 1676 & 'NXTH2X error, KFREE lower than internal allocation' 1677 WRITE (LUERR,*) 'KFREE ',KFREE 1678 WRITE (LUERR,*) 'KICDTR',KICDTR 1679 WRITE (LUERR,*) 'KNEXT ',KNEXT, 1680 & ' ( next avail. address after internal allocation)' 1681 END IF 1682 CALL MEMCHK('IDIST .ne. 0 MEMCHK in NXTH2X',WRK,KICDTR) 1683 IF (KFREE.LT.KNEXT) THEN 1684 CALL QTRACE(LUERR) 1685 CALL QUIT('NXTH2X error: KFREE '// 1686 & 'lower than internal allocation') 1687 END IF 1688 END IF 1689C 1690 CALL NX2H2X(IC,ID,H2CD,NEEDTP,IDIST, LUDAX,LSYMX,N2DIS, 1691 * WRK(KICDTR),WRK(KH2CDP)) 1692C CALL NX2H2X(IC,ID,H2CD,NEEDTP,IDIST, LUDAX,LSYMX,N2DIS, 1693C * ICDTRA,H2CDPK) 1694C 1695C If no more integrals (IDIST .lt. 0) then release internal space 1696C 1697 IF (IDIST .LT. 0) THEN 1698 CALL MEMREL('Releasing internal space in NXTH2X',WRK,KICDTR, 1699 & KICDTR,KFREE,LFREE) 1700 KNEXT = -1 1701 END IF 1702 9999 CALL QEXIT('NXTH2X') 1703 RETURN 1704 END 1705C /* Deck nx2h2x */ 1706 SUBROUTINE NX2H2X(IC,ID,H2CD,NEEDTP,IDIST, 1707 * LUDAX,LSYMX,N2DIS,ICDTRA,H2CDPK) 1708C 1709C Copyright November 1990 by Hans Joergen Aa. Jensen 1710C 1711C Purpose: 1712C 1713C Read next Mulliken two-electron integral distribution (**|cd) 1714C where (cd) distribution is needed according to NEEDTP(ITYPCD) 1715C 1716C Input: 1717C NEEDTP(i); positive for needed (cd) distribution types 1718C negative if not all distributions needed 1719C zero if no distributions needed for this type 1720C IDIST; .eq. 0 first read 1721C .gt. 1 intermediate read 1722C .lt. 0 end-of-file has been reached previously 1723C Output: 1724C H2CD(NORBT,NORBT); H2CD(a,b) = (ab|cd) 1725C IC,ID; value of c and d 1726C IDIST; .gt. 0 when next distribution IC,ID available in H2CD 1727C = -1 when no more distributions 1728C Internal: 1729C ICDTRA(ICD) .ne. 0 if (**|cd) distribution has been transformed. 1730C H2CDPK(*) is used for packed integrals on disk. 1731C 1732C **************************************************************** 1733C 1734#include "implicit.h" 1735 DIMENSION H2CD(NORBT,NORBT), H2CDPK(N2ORBT) 1736 INTEGER NEEDTP(-4:6), ICDTRA(NORBT,NORBT), N2DIS(8) 1737C 1738C 1739C Used from common blocks: 1740C INFORB : N2ORBX,NORBT,NSYM,... 1741C INFIND : IOBTYP(*),? 1742C INFPRI : IPRSIR 1743C 1744#include "maxash.h" 1745#include "maxorb.h" 1746#include "priunit.h" 1747#include "inforb.h" 1748#include "infind.h" 1749#include "infpri.h" 1750C 1751#include "orbtypdef.h" 1752C 1753 SAVE ICOLD,IDOLD 1754C 1755C 1756C **************************************************************** 1757C 1758C If first read (IDIST .EQ. 0) 1759C then initialize ... 1760C 1761 IF (IDIST .EQ. 0) THEN 1762 ICOLD = 0 1763 IDOLD = 1 1764 IF (IPRSIR .GT. 20) THEN 1765 WRITE (LUPRI,'(/A//A,8I8)') 1766 * ' Test output from NX2H2X for IDIST = 0', 1767 * ' N2DIS array :',(N2DIS(I),I=1,NSYM) 1768 END IF 1769 IF (IPRSIR .GT. 50) THEN 1770 WRITE (LUPRI,'(/A)') ' ICDTRA matrix:' 1771 DO 10 I = 1,NORBT 1772 WRITE (LUPRI,'(I5,A,(T8,10I7))') I,' :', 1773 * (ICDTRA(I,J),J=1,NORBT) 1774 10 CONTINUE 1775 END IF 1776 END IF 1777C 1778C *** Read next distribution which is needed according to NEEDTP(*) 1779C into H2CD 1780C 1781C ITYPCD values: 1=i*i : 2=t*i : 3=t*t : 4=a*i : 5=a*t : 6=a*a 1782C 0 for not wanted type. 1783C 1784C 1785 ICNEW = ICOLD 1786 IDNEW = IDOLD 1787 200 CONTINUE 1788 ICNEW = ICNEW + 1 1789 IF (ICNEW .GT. NORBT) THEN 1790 IDNEW = IDNEW + 1 1791 ICNEW = 1 1792 END IF 1793 IF (IDNEW .GT. NORBT) THEN 1794C Last distribution has been read 1795 IDIST = -1 1796 GO TO 9999 1797 END IF 1798C 1799 ITYPC = IOBTYP(ICNEW) 1800 ITYPD = IOBTYP(IDNEW) 1801 ITYPCD = IDBTYP(ITYPC,ITYPD) 1802 IF ( NEEDTP(ITYPCD) .EQ. 0 ) THEN 1803 ITYPCD = 0 1804 ELSE IF (NEEDTP(ITYPCD) .LT. 0) THEN 1805 ITYPCD = -ITYPCD 1806 END IF 1807 IDIST = ICDTRA(ICNEW,IDNEW) 1808C 1809 IF (IPRSIR .GT. 50) THEN 1810 WRITE(LUPRI,*) 'From NX2H2X:' 1811 WRITE(LUPRI,*) 'ICNEW ,IDNEW :',ICNEW,IDNEW 1812 WRITE(LUPRI,*) 'ICDTRA(ICNEW,IDNEW) :',ICDTRA(ICNEW,IDNEW) 1813 WRITE(LUPRI,*) 'ITYPCD :',ITYPCD 1814 END IF 1815C 1816 IF (ITYPCD .GT. 0 .AND. IDIST .EQ. 0) THEN 1817 WRITE (LUERR,*) ' NX2H2X ERROR: needed integral distribution' 1818 WRITE (LUERR,*) ' has not been calculated' 1819 WRITE (LUERR,*) 'IC ,ID :',ICNEW,IDNEW 1820 WRITE (LUERR,*) 'IYPC ,ITYPD :',COBTYP(ITYPC),COBTYP(ITYPD) 1821 CALL QTRACE(LUERR) 1822 CALL QUIT('NX2H2X error: needed integrals not calculated') 1823 END IF 1824C 1825 IF (ITYPCD .EQ. 0) GO TO 200 1826C 1827 ISYMC = ISMO(ICNEW) 1828 ISYMD = ISMO(IDNEW) 1829 ISYMCD = MULD2H(ISYMC,ISYMD) 1830 ISYMAB = MULD2H(ISYMCD,LSYMX) 1831 N2DAB = N2DIS(ISYMAB) 1832C 1833 READ(LUDAX,REC=IDIST) (H2CDPK(I),I=1,N2DAB) 1834C 1835C Unpack integrals 1836C 1837 CALL OITPAK(H2CD,H2CDPK,ISYMAB,-1) 1838C 1839C******************************************************************* 1840C 1841C End of subroutine NX2H2X 1842C 1843 9999 CONTINUE 1844 ICOLD = ICNEW 1845 IDOLD = IDNEW 1846 IC = ICNEW 1847 ID = IDNEW 1848 RETURN 1849 END 1850