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 /* deck ccder1 */ 20*=====================================================================* 21 SUBROUTINE CCDER1(IATOM,LABELOP,LDERINT,WORK,LWORK) 22*---------------------------------------------------------------------* 23* 24* Purpose: get first derivative integrals from ABACUS 25* and sort them into the ordering used in CC 26* 27* Written by Christof Haettig, 06-May-1998 28* 29*=====================================================================* 30#if defined (IMPLICIT_NONE) 31 IMPLICIT NONE 32#else 33# include "implicit.h" 34#endif 35#include "priunit.h" 36#include "mxcent.h" 37#include "nuclei.h" 38#include "ccorb.h" 39 40* local parameters: 41 CHARACTER*(19) MSGDBG 42 PARAMETER (MSGDBG = '[debug] CCDER1> ') 43 LOGICAL LOCDBG 44 PARAMETER (LOCDBG = .false.) 45 46 LOGICAL LDERINT(8,3) 47 CHARACTER*8 LABELOP 48 INTEGER LWORK 49 INTEGER IATOM, NDMAT, MAXDIF, MXCOMP 50 INTEGER KDMAT, KEND, KFMAT, KISYMDM, KIFCTYP, LEND 51 52#if defined (SYS_CRAY) 53 REAL WORK(LWORK) 54 REAL TIM0 55#else 56 DOUBLE PRECISION WORK(LWORK) 57#endif 58 59 60 CALL QENTER('CCDER1') 61*---------------------------------------------------------------------* 62* check, if not the same integrals have been calculated the last time: 63*---------------------------------------------------------------------* 64C IF (LABELOP .EQ. LAST_LABELOP) THEN 65C CALL QEXIT('CCDER1') 66C RETURN 67C END IF 68 69*---------------------------------------------------------------------* 70* some intializations 71*---------------------------------------------------------------------* 72 IF (LOCDBG) WRITE (LUPRI,*) MSGDBG, 'entered CCDER1.' 73 74 KEND = 1 ! work space 75 76 NDMAT = 0 ! # densityies for fock matrices 77 78 MAXDIF = 1 ! order for derivatives --> first derivatives 79 80 MXCOMP = 3 ! max. comp. per atom 81 82 CALL RHSINI ! initialize some ABACUS common blocks 83*---------------------------------------------------------------------* 84* begin: 85*---------------------------------------------------------------------* 86 KDMAT = KEND 87 KFMAT = KDMAT + NDMAT*N2BASX 88 KISYMDM = KFMAT + NDMAT*MXCOMP*NUCDEG(IATOM)*N2BASX 89 KIFCTYP = KISYMDM + NDMAT*MXCOMP*NUCDEG(IATOM) 90 KEND = KIFCTYP + NDMAT*MXCOMP*NUCDEG(IATOM) 91 LEND = LWORK - KEND + 1 92 93 CALL CCGETH2D(IATOM,MAXDIF,LDERINT,LABELOP, 94 & WORK(KDMAT),WORK(KFMAT), NDMAT, 95 & WORK(KISYMDM),WORK(KIFCTYP), MXCOMP, 96 & WORK(KEND),LEND) 97 98C LAST_LABELOP = LABELOP 99 100*---------------------------------------------------------------------* 101* return: 102*---------------------------------------------------------------------* 103 CALL FLSHFO(LUPRI) 104 105 CALL QEXIT('CCDER1') 106 RETURN 107 END 108 109*=====================================================================* 110* END OF SUBROUTINE CCDER1 * 111*=====================================================================* 112*======================================================================* 113C /* Deck ccsortderao */ 114 SUBROUTINE CCSORTDERAO(WORK,LWORK,MXCOMP,LDERINT,ITYPE,IPRINT) 115*----------------------------------------------------------------------* 116C 117C Purpose: sort derivative 2-el. AO integrals. 118C 119C MXCOMP: maximum number of 'components' 120C --> three (x,y,z) for first derivative integrals 121C 122C ITYPE : 1 - geometric first derivatives 123C 5 - magnetic first derivatives 124C 125C Written by Christof Haettig 06-May-1998 126C based on Henrik Koch's CCSD_SORTAO routine 127C magnetic derivatives added Sep-1999 128C 129*----------------------------------------------------------------------* 130 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 131#include "priunit.h" 132#include "iratdef.h" 133#include "inftap.h" 134C 135 LOGICAL LOCDBG ! local debug flag 136 PARAMETER (LOCDBG = .false.) 137C 138 DIMENSION WORK(LWORK) 139 LOGICAL LDERINT(8,MXCOMP) 140C 141 INTEGER LENGTH(8), LUAODR(8) 142C 143 CHARACTER*8 NAME(8) 144C 145 DATA NAME /'CCAODER1','CCAODER2','CCAODER3','CCAODER4', 146 * 'CCAODER5','CCAODER6','CCAODER7','CCAODER8'/ 147#include "ccorb.h" 148#include "ccsdsym.h" 149#include "eribuf.h" 150 151*----------------------------------------------------------------------* 152* open files for sorted integrals: 153*----------------------------------------------------------------------* 154 DO ISYM = 1,NSYM 155 LUAODR(ISYM) = -1 156 CALL WOPEN2(LUAODR(ISYM),NAME(ISYM),64,0) 157 END DO 158 159*----------------------------------------------------------------------* 160* set buffer information. 161*----------------------------------------------------------------------* 162 CALL ERIBUF_INI ! set NIBUF, NBITS, IBIT1, IBIT2 163 LBUF = 600 164 165*----------------------------------------------------------------------* 166* Buffer allocation. 167*----------------------------------------------------------------------* 168 KRBUF = 1 169 KIBUF = KRBUF + LBUF 170 KAOAB = KIBUF + (NIBUF*LBUF + 1)/2 + 1 ! IBUF always integer*4 171 KAOG = KAOAB + (N2BASX + 1)/IRAT + 1 172 KEND1 = KAOG + (NBAST*NSYM*NSYM*MXCOMP + 1)/IRAT + 1 173 LWRK1 = LWORK - KEND1 174C 175 IF (LWRK1 .LT. 0) 176 & CALL QUIT('Insufficient work space in CCSORTDERAO.') 177*----------------------------------------------------------------------* 178* Calculate in the index arrays needed in the sort. 179*----------------------------------------------------------------------* 180 181 CALL CCSD_INIT2B(WORK(KAOAB),WORK(KAOG),WORK(KAOG),ITYPE, 182 & .FALSE.,MXCOMP,LDERINT) 183 184*----------------------------------------------------------------------* 185* precalculate length of integral (* *|* del) distributions: 186*----------------------------------------------------------------------* 187 DO ISYMD = 1, NSYM 188 LENGTH(ISYMD) = 0 189 DO ICOOR = 1, MXCOMP 190 DO ICORSY = 1, NSYM 191 IF ( LDERINT(ICORSY,ICOOR) ) THEN 192 ISYDIS = MULD2H(ISYMD,ICORSY) 193 194 IF (ITYPE.EQ.0 .OR. ITYPE.EQ.1) THEN 195 LENGTH(ISYMD) = LENGTH(ISYMD) + NDISAO(ISYDIS) 196 ELSE IF (ITYPE.EQ.5) THEN 197 LENGTH(ISYMD) = LENGTH(ISYMD) + NDISAOSQ(ISYDIS) 198 ELSE 199 CALL QUIT('Illegal ITYPE in CCSORTDERAO.') 200 END IF 201 202 END IF 203 END DO 204 END DO 205 206 END DO 207 208*----------------------------------------------------------------------* 209* Loop over batches of integrals. 210*----------------------------------------------------------------------* 211C 212 DO 100 ISYMD = 1,NSYM 213C 214 IOFF2 = 1 215C 216 NTOTD = NBAS(ISYMD) 217 IF (NTOTD .EQ. 0) GOTO 100 218 219 NUMBAT = MIN(NTOTD,LWRK1/LENGTH(ISYMD)) 220C 221 IF (NUMBAT .EQ. 0) THEN 222 WRITE(LUPRI,*) 'In CCSD_SORTAO NUMBAT is zero' 223 CALL QUIT('Insufiicient work space in CCRDAO') 224 ENDIF 225C 226 ITOTBA = (NTOTD-1)/NUMBAT + 1 227C 228 ID1 = IBAS(ISYMD) + 1 229 ID2 = IBAS(ISYMD) 230 IOFF1 = IBAS(ISYMD) 231C 232 DO 200 I = 1,ITOTBA 233C 234 INUMBA = NUMBAT 235 IF (NUMBAT*I .GT. NTOTD) THEN 236 INUMBA = NTOTD - NUMBAT*(I-1) 237 ENDIF 238C 239 ID2 = ID2 + INUMBA 240C 241 CALL DZERO(WORK(KEND1),LENGTH(ISYMD)*INUMBA) 242C 243 CALL CCSORTDER1(WORK(KEND1),WORK(KIBUF),WORK(KRBUF), 244 * WORK(KAOAB),WORK(KAOG),ISYMD,LENGTH(ISYMD), 245 * IOFF1,ID1,ID2,NIBUF,LBUF,NBITS,IBIT1,ITYPE) 246C 247c CALL AROUND('Norm of integral distributions:') 248c IPRC = KEND1 249c DO IPRD = ID1,ID2 250c WRITE (LUPRI,*) 'D distribution',IPRD 251c DO ICOOR = 1, MXCOMP 252c DO ICORSY = 1, NSYM 253c IF ( LDERINT(ICORSY,ICOOR) ) THEN 254c WRITE (LUPRI,*) 'coordinate/symmetry',ICOOR,ICORSY 255c ISYDIS = MULD2H(ISYMD,ICORSY) 256c XNORM = 0.0d0 257c DO IPSYMG = 1,NSYM 258c ISYMAB = MULD2H(IPSYMG,ISYDIS) 259c IF (ITYPE.EQ.0 .OR. ITYPE.EQ.1) THEN 260c LEN = NNBST(ISYMAB) * NBAS(IPSYMG) 261c XNORM=XNORM+DDOT(LEN,WORK(IPRC),1,WORK(IPRC),1) 262c IPRC = IPRC + NNBST(ISYMAB)*NBAS(IPSYMG) 263c ELSE IF (ITYPE.EQ.5) THEN 264c LEN = N2BST(ISYMAB) * NBAS(IPSYMG) 265c XNORM=XNORM+DDOT(LEN,WORK(IPRC),1,WORK(IPRC),1) 266c DO G = 1, NBAS(IPSYMG) 267c IGAM = G + IBAS(IPSYMG) 268c KOFFX = IPRC + N2BST(ISYMAB)*(G-1) 269c YNORM = DDOT(N2BST(ISYMAB),WORK(KOFFX),1, 270c * WORK(KOFFX),1) 271c WRITE (LUPRI,'(a,5i5,f10.5)'),'CCSORTDERAO> ',IPRD 272c * IGAM,ICOOR,ICORSY,IOFF2+KOFFX-KEND1,YNORM 273c END DO 274c IPRC = IPRC + N2BST(ISYMAB)*NBAS(IPSYMG) 275c ELSE 276c CALL QUIT('Unknown ITYPE in CCSORTDERAO.') 277c END IF 278c END DO 279c WRITE (LUPRI,*) 'Norm of distribution:',XNORM 280c 281c END IF 282c END DO 283c END DO 284c END DO 285C 286 CALL CCSORTDER2(WORK(KEND1),IOFF2,INUMBA, 287 * LENGTH(ISYMD),ISYMD,NAME,LUAODR) 288C 289 IF ( (IPRINT.GT.50) .OR. LOCDBG ) THEN 290 CALL AROUND('Integral distribution') 291 IPRC = KEND1 292 DO IPRD = ID1,ID2 293 WRITE (LUPRI,*) 'D distribution',IPRD 294 DO ICOOR = 1, MXCOMP 295 DO ICORSY = 1, NSYM 296 IF ( LDERINT(ICORSY,ICOOR) ) THEN 297 WRITE (LUPRI,*) 'coordinate/symmetry', 298 & ICOOR,ICORSY 299 ISYDIS = MULD2H(ISYMD,ICORSY) 300 DO IPSYMG = 1,NSYM 301 WRITE(LUPRI,*) 'Gamma symmetry',IPSYMG 302 ISYMAB = MULD2H(IPSYMG,ISYDIS) 303 304 IF (LOCDBG) THEN 305 306 DO G = 1, NBAS(IPSYMG) 307 IG = IBAS(IPSYMG) + G 308 DO ISYMB = 1, NSYM 309 ISYMA = MULD2H(ISYMB,ISYMAB) 310 WRITE (LUPRI,*) 'symmet.:', 311 * ISYMA,ISYMB,IPSYMG 312 IF (((ITYPE.EQ.0 .OR. ITYPE.EQ.1) 313 * .AND. (ISYMB .GT. ISYMA) ) 314 * .OR. (ITYPE.EQ.5) )THEN 315 DO B = 1, NBAS(ISYMB) 316 IB = IBAS(ISYMB) + B 317 DO A = 1, NBAS(ISYMA) 318 IA = IBAS(ISYMA) + A 319 WRITE (LUPRI,'(4I5,G20.10,I10)') 320 * IA,IB,IG,IPRD,WORK(IPRC), 321 * IPRC-KEND1+1 322 IPRC = IPRC + 1 323 END DO 324 END DO 325 ELSE IF ((ITYPE.EQ.0 .OR. ITYPE.EQ.1) 326 * .AND. (ISYMB.EQ.ISYMA)) THEN 327 DO B = 1, NBAS(ISYMB) 328 IB = IBAS(ISYMB) + B 329 DO A = 1, B 330 IA = IBAS(ISYMA) + A 331 WRITE (LUPRI,'(4I5,G20.10,I10)') 332 * IA,IB,IG,IPRD,WORK(IPRC), 333 * IPRC-KEND1+1 334 IPRC = IPRC + 1 335 END DO 336 END DO 337 ELSE IF ((ITYPE.EQ.0 .OR. ITYPE.EQ.1) 338 * .AND. (ISYMB.LT.ISYMA)) THEN 339 CONTINUE 340 ELSE 341 CALL QUIT( 342 * 'Unknown ITYPE in CCSORTDERAO.') 343 END IF 344 END DO 345 END DO 346 347 ELSE 348 IF (ITYPE.EQ.0 .OR. ITYPE.EQ.1) THEN 349 CALL OUTPUT(WORK(IPRC),1,NNBST(ISYMAB),1, 350 * NBAS(IPSYMG),NNBST(ISYMAB), 351 * NBAS(IPSYMG),1,LUPRI) 352 IPRC = IPRC + NNBST(ISYMAB)*NBAS(IPSYMG) 353 ELSE IF (ITYPE.EQ.5) THEN 354 CALL OUTPUT(WORK(IPRC),1,N2BST(ISYMAB),1, 355 * NBAS(IPSYMG),N2BST(ISYMAB), 356 * NBAS(IPSYMG),1,LUPRI) 357 IPRC = IPRC + N2BST(ISYMAB)*NBAS(IPSYMG) 358 ELSE 359 CALL QUIT( 360 * 'Unknown ITYPE in CCSORTDERAO.') 361 END IF 362 END IF 363 END DO 364 END IF 365 END DO 366 END DO 367 END DO 368 END IF 369C 370 ID1 = ID1 + INUMBA 371 IOFF1 = IOFF1 + INUMBA 372C 373 200 CONTINUE 374C 375 100 CONTINUE 376C 377 DO ISYM = 1,NSYM 378 CALL WCLOSE2(LUAODR(ISYM),NAME(ISYM),'KEEP') 379 END DO 380C 381 CALL GPCLOSE(LU2DER,'DELETE') 382C 383 RETURN 384 END 385*======================================================================* 386C /* Deck ccsortder1 */ 387 SUBROUTINE CCSORTDER1(XINT,IBUF4,RBUF,KAOAB,KAOGDER,ISYMD,LENGTH, 388 * IOFF,ID1,ID2,NIBUF,LBUF,NBITS,IBIT1,ITYPE) 389*----------------------------------------------------------------------* 390C 391C Written by Christof Haettig 06-May-1998 392C based on Henrik Kochs CCSD_SORT1 routine 393C generalizations for London integrals Sep-1999 394C 395*----------------------------------------------------------------------* 396#include "implicit.h" 397#include "priunit.h" 398#include "dummy.h" 399#include "ibtpar.h" 400#include "ccorb.h" 401 DIMENSION XINT(LENGTH,*),RBUF(LBUF) 402 INTEGER*4 IBUF4(LBUF*NIBUF), LENGTH4 403 INTEGER INDX4(4,LBUF) 404 DIMENSION KAOAB(NBAST,NBAST),KAOGDER(NBAST,NSYM,NSYM,*) 405 INTEGER ITYPE 406C 407#include "inftap.h" 408 409C 410 CALL REWSPL(LU2DER) 411C 412 IF (ITYPE.EQ.5) THEN 413 CALL MOLLAB('AO2MGINT',LU2DER,LUPRI) 414 END IF 415C 416 IF (ITYPE.EQ.0 .OR. ITYPE.EQ.1) THEN 417 SIGN = +1.0D0 418 ELSE IF (ITYPE.EQ.5) THEN 419 SIGN = -1.0D0 420 ELSE 421 CALL QUIT('Unknown ITYPE in CCSORTDER1.') 422 END IF 423C 424C 425 10 READ(LU2DER,ERR=2000) RBUF,IBUF4,LENGTH4 426 LENGTH = LENGTH4 427C 428 IF (LENGTH .GT. 0) THEN 429 CALL AOLAB4_cc(IBUF4,NIBUF,NBITS,INDX4,LENGTH) 430 431 DO I = 1, LENGTH 432 433 IP = INDX4(4,I) 434 435 IF (IP .EQ. 0) THEN 436 437* ... FLAG : ICOOR,ICORSY... 438 ICOOR = INDX4(3,I) 439 ICORSY = INDX4(2,I) + 1 440 441 ELSE 442 443* ... INTEGRAL : (IP IQ | IR IS) ... 444 IQ = INDX4(3,I) 445 IR = INDX4(2,I) 446 IS = INDX4(1,I) 447C 448 IF ((IS .GE. ID1) .AND. (IS .LE. ID2)) THEN 449 ! store as (PQ|RS) 450 IADR = KAOGDER(IR,ISYMD,ICORSY,ICOOR) + KAOAB(IP,IQ) 451 XINT(IADR,IS-IOFF) = SIGN * RBUF(I) 452 !IF(ITYPE.EQ.5)WRITE (LUPRI,*)IP,IQ,IR,IS,XINT(IADR,IS-IOFF) 453 ENDIF 454 IF ((IR .GE. ID1) .AND. (IR .LE. ID2)) THEN 455 ! store as (QP|SR) 456 IADR = KAOGDER(IS,ISYMD,ICORSY,ICOOR) + KAOAB(IQ,IP) 457 XINT(IADR,IR-IOFF) = RBUF(I) 458 !IF(ITYPE.EQ.5)WRITE (LUPRI,*)IQ,IP,IS,IR,XINT(IADR,IR-IOFF) 459 ENDIF 460 IF ((IP .GE. ID1) .AND. (IP .LE. ID2)) THEN 461 ! store as (SR|QP) 462 IADR = KAOGDER(IQ,ISYMD,ICORSY,ICOOR) + KAOAB(IS,IR) 463 XINT(IADR,IP-IOFF) = RBUF(I) 464 !IF(ITYPE.EQ.5)WRITE (LUPRI,*)IS,IR,IQ,IP,XINT(IADR,IP-IOFF) 465 ENDIF 466 IF ((IQ .GE. ID1) .AND. (IQ .LE. ID2)) THEN 467 ! store as (RS|PQ) 468 IADR = KAOGDER(IP,ISYMD,ICORSY,ICOOR) + KAOAB(IR,IS) 469 XINT(IADR,IQ-IOFF) = SIGN * RBUF(I) 470 !IF(ITYPE.EQ.5)WRITE (LUPRI,*)IR,IS,IP,IQ,XINT(IADR,IQ-IOFF) 471 ENDIF 472 473 END IF 474C 475 END DO 476 477 ELSE IF (LENGTH .LT. 0) THEN 478 GOTO 100 479 END IF 480C 481 GOTO 10 482C 483 100 CONTINUE 484C 485 RETURN 486 487 2000 CALL QUIT('Error reading derivative integrals in CCSORTDER1') 488 489 END 490*======================================================================* 491C /* Deck ccsortder2 */ 492 SUBROUTINE CCSORTDER2(XINT,IOFF,INUMBA,LENGTH,ISYM,NAME,LUAODR) 493*----------------------------------------------------------------------* 494C 495C Written by Christof Haettig 06-May-1998 496C based on Henrik Kochs CCSD_SORT2 routine 497C 498*----------------------------------------------------------------------* 499#include "implicit.h" 500 DIMENSION XINT(*) 501C 502 INTEGER LUAODR(8) 503 CHARACTER*8 NAME(8) 504C 505 CALL PUTWA2(LUAODR(ISYM),NAME(ISYM),XINT,IOFF,LENGTH*INUMBA) 506 IOFF = IOFF + LENGTH*INUMBA 507C 508 RETURN 509 END 510*======================================================================* 511C /* Deck ccsd_init2b */ 512 SUBROUTINE CCSD_INIT2B(KAOAB,KAOG,KAOGDER,ITYPE,LDISTRIB, 513 & MXCOMP,LDERINT) 514*----------------------------------------------------------------------* 515C 516C Henrik Koch and Alfredo Sanchez. 29-Jun-1994 517C derivative option by Christof Haettig 06-May-1998 518C 519C Set up indexing arrays for integral sort 520C 521C ITYPE = 0 : set up KAOG for undiff. 2-el integrals 522C (KAOGDER ignored) 523C 524C ITYPE = 1 : set up KAOGDER for geom. 1. derivative 2-el integrals 525C (KAOG ignored) 526C 527C ITYPE = 5 : set up KAOGDER for magn. 1. derivative 2-el integrals 528C (KAOG ignored) 529C 530C LDISTRIB : if true only KAOAB is set up 531C (used to sort distributions) 532C 533C LDERINT : array of logicals that tell if there are any 534C derivative integrals for a given component 535C and given symmetry 536C 537*----------------------------------------------------------------------* 538#include "implicit.h" 539#include "ccorb.h" 540 DIMENSION KAOAB(NBAST,NBAST) 541 DIMENSION KAOG(NBAST,NSYM) 542 DIMENSION KAOGDER(NBAST,NSYM,NSYM,MXCOMP) 543 LOGICAL LDISTRIB, LDERINT(8,MXCOMP) 544 DIMENSION NDIMAB(8) 545#include "ccsdsym.h" 546 547*----------------------------------------------------------------------* 548* set up KAOAB index array for the leading indeces alpha,beta : 549*----------------------------------------------------------------------* 550 551 IF ( ITYPE.EQ.0 .OR. ITYPE.EQ.1 ) THEN 552C -------------------------------------------------------- 553C for undiff. Hamiltonian or geometric derivatives pack 554C fast running indeces alpha,beta in triangular form: 555C -------------------------------------------------------- 556 DO ISYMAB = 1,NSYM 557 NDIMAB(ISYMAB) = NNBST(ISYMAB) 558 ICOUNT = 0 559 DO ISYMB = 1,NSYM 560 ISYMA = MULD2H(ISYMB,ISYMAB) 561 IF (ISYMB .GT. ISYMA) THEN 562 DO B = 1,NBAS(ISYMB) 563 IB = IBAS(ISYMB) + B 564 DO A = 1,NBAS(ISYMA) 565 IA = IBAS(ISYMA) + A 566 ICOUNT = ICOUNT + 1 567 KAOAB(IA,IB) = ICOUNT 568 KAOAB(IB,IA) = ICOUNT 569 END DO 570 END DO 571 ELSE IF (ISYMA .EQ. ISYMB) THEN 572 DO B = 1,NBAS(ISYMB) 573 IB = IBAS(ISYMB) + B 574 DO A = 1,B 575 IA = IBAS(ISYMA) + A 576 ICOUNT = ICOUNT + 1 577 KAOAB(IA,IB) = ICOUNT 578 KAOAB(IB,IA) = ICOUNT 579 END DO 580 END DO 581 END IF 582 END DO 583 END DO 584 ELSE 585C -------------------------------------------------------- 586C if the Hamiltonian is not real (London orbitials) pack 587C fast running indeces alpha,beta in squared form: 588C -------------------------------------------------------- 589 DO ISYMAB = 1,NSYM 590 NDIMAB(ISYMAB) = N2BST(ISYMAB) 591 DO ISYMB = 1,NSYM 592 ISYMA = MULD2H(ISYMB,ISYMAB) 593 ICOUNT = IAODIS(ISYMA,ISYMB) 594 DO B = 1,NBAS(ISYMB) 595 IB = IBAS(ISYMB) + B 596 DO A = 1,NBAS(ISYMA) 597 IA = IBAS(ISYMA) + A 598 ICOUNT = ICOUNT + 1 599 KAOAB(IA,IB) = ICOUNT 600 END DO 601 END DO 602 END DO 603 END DO 604 END IF 605 606 607*----------------------------------------------------------------------* 608* set up KAOG/KAOGDER index array for the third index gamma: 609*----------------------------------------------------------------------* 610 IF (.NOT.LDISTRIB) THEN 611 612 IF (ITYPE.EQ.0) THEN 613C -------------------------------------------------------- 614C for undifferentiated integrals set up KAOG array: 615C -------------------------------------------------------- 616 DO ISYMD = 1,NSYM 617 ISYDIS = MULD2H(ISYMD,ISYMOP) 618 ICOUNT = 0 619 DO ISYMG = 1,NSYM 620 ISYMAB = MULD2H(ISYMG,ISYDIS) 621 DO G = 1,NBAS(ISYMG) 622 IG = IBAS(ISYMG) + G 623 KAOG(IG,ISYMD) = ICOUNT 624 ICOUNT = ICOUNT + NDIMAB(ISYMAB) 625 END DO 626 END DO 627 END DO 628C 629 ELSE IF (ITYPE.EQ.1 .OR. ITYPE.EQ.5) THEN 630C -------------------------------------------------------- 631C for differentiated integrals set up KAOGDER array: 632C -------------------------------------------------------- 633 DO ISYMD = 1,NSYM 634 ICOUNT = 0 635 DO ICOOR = 1, MXCOMP 636 DO ICORSY = 1, NSYM 637 IF ( LDERINT(ICORSY,ICOOR) ) THEN 638 ISYDIS = MULD2H(ISYMD,ICORSY) 639 DO ISYMG = 1,NSYM 640 ISYMAB = MULD2H(ISYMG,ISYDIS) 641 DO G = 1,NBAS(ISYMG) 642 IG = IBAS(ISYMG) + G 643 KAOGDER(IG,ISYMD,ICORSY,ICOOR) = ICOUNT 644 ICOUNT = ICOUNT + NDIMAB(ISYMAB) 645 END DO 646 END DO 647 END IF 648 END DO 649 END DO 650 END DO 651C 652 ELSE 653 CALL QUIT('Unknown ITYPE in CCSD_INIT2B.') 654 END IF 655C 656 END IF 657 658 RETURN 659 END 660*======================================================================* 661c /* deck ccgeth2d */ 662*=====================================================================* 663 SUBROUTINE CCGETH2D(IATOM, MAXDIF, LDERINT, LABELOP, 664 & DMAT, FMAT, NDMAT, ISYMDM, IFCTYP, MXCOMP, 665 & WORK, LWORK) 666*---------------------------------------------------------------------* 667* 668* Purpose: call ABACUS to calculate derivatives of 2-el integrals 669* and sort them to CC storage scheme 670* 671* MAXDIF : max. derivative order 672* LABELOP : operator label 673* 674* Christof Haettig, May 1998 675*=====================================================================* 676#if defined (IMPLICIT_NONE) 677 IMPLICIT NONE 678#else 679# include "implicit.h" 680#endif 681#include "priunit.h" 682#include "ccsdinp.h" 683#include "iratdef.h" 684#include "aovec.h" 685#include "mxcent.h" 686#include "maxaqn.h" 687#include "maxorb.h" 688#include "nuclei.h" 689#include "nodint.h" 690#include "ccorb.h" 691#include "cch2d.h" 692#include "inftap.h" 693#include "dummy.h" 694#include "second.h" 695 696* local parameters: 697 CHARACTER*(19) MSGDBG 698 PARAMETER (MSGDBG = '[debug] CCGETH2D> ') 699 LOGICAL LOCDBG 700 PARAMETER (LOCDBG = .false.) 701 INTEGER NBUFS 702 PARAMETER ( NBUFS = 600 ) 703 704 LOGICAL LDERINT(8,3) 705 CHARACTER*8 LABELOP 706 INTEGER IATOM, MAXDIF, NUMDIS, NDMAT 707 INTEGER LWORK 708 INTEGER ISYMDM(*), IFCTYP(*) 709 710 711 LOGICAL RELCAL 712 INTEGER MAXDIS 713 INTEGER IBUF(NBUFS) 714 INTEGER KJSTRS, KNPRIM, KNCONT, KIORBS, KJORBS, KKORBS, KLAST 715 INTEGER NFMAT, LFMAT, K2INT, L2INT, ITYPE, I2TYP 716 INTEGER ICOOR, ICORSY, JR, JS, JP, JQ, JRS, JPQ, I 717 INTEGER MXCOMP, ISYM, IPRHER 718 719#if defined (SYS_CRAY) 720 REAL WORK(LWORK) 721 REAL DMAT(*), FMAT(*) 722 REAL TIM0, TIM1, TIM2, TIM3 723#else 724 DOUBLE PRECISION WORK(LWORK) 725 DOUBLE PRECISION DMAT(*), FMAT(*) 726 DOUBLE PRECISION TIM0, TIM1, TIM2, TIM3 727#endif 728 729* external functions 730#if defined (SYS_CRAY) 731 REAL BUF(NBUFS) 732#else 733 DOUBLE PRECISION BUF(NBUFS) 734#endif 735 736 CALL QENTER('CCGETH2D') 737*---------------------------------------------------------------------* 738* begin 739*---------------------------------------------------------------------* 740 IF (DEBUG.OR.LOCDBG) WRITE (LUPRI,*) MSGDBG, 'entered CCGETH2D.' 741 742* initialize timing: 743 TIM0 = SECOND() 744 745* initialize address for next free element on work: 746 KLAST = 1 747 748* set print level for integral part: 749C IPRHER = IPRINT / 5 750 IPRHER = 0 751 752*---------------------------------------------------------------------* 753* allocate & initialize some vectors for TWOINT 754*---------------------------------------------------------------------* 755 KJSTRS = KLAST 756 KNPRIM = KJSTRS + (MXSHEL*MXAOVC*2 + 1)/IRAT 757 KNCONT = KNPRIM + (MXSHEL*MXAOVC*2 + 1)/IRAT 758 KIORBS = KNCONT + (MXSHEL*MXAOVC*2 + 1)/IRAT 759 KJORBS = KIORBS + (MXSHEL*MXAOVC + 1)/IRAT 760 KKORBS = KJORBS + (MXSHEL*MXAOVC + 1)/IRAT 761 KLAST = KKORBS + (MXSHEL*MXAOVC + 1)/IRAT 762 IF (KLAST .GT. LWORK) THEN 763 CALL QUIT('Insufficient work space in CCGETH2D.') 764 END IF 765 766 CALL PAOVEC(WORK(KJSTRS),WORK(KNPRIM),WORK(KNCONT), 767 & WORK(KIORBS),WORK(KJORBS),WORK(KKORBS), 768 & IATOM,.FALSE.,IPRHER) 769 770 KLAST = KJORBS 771 772 773*---------------------------------------------------------------------* 774* initialize Fock matrices and related integer arrays ISYMDM,IFCTYP: 775*---------------------------------------------------------------------* 776 NFMAT = MXCOMP*NDMAT*NUCDEG(IATOM) 777 LFMAT = NFMAT * N2BASX 778 779 IF (NFMAT.GT.0) THEN 780 IF (LFMAT.GT.0) CALL DZERO(FMAT,LFMAT) 781 DO I = 1, NFMAT 782 ISYMDM(I) = 0 783 IFCTYP(I) = 13 784 END DO 785 END IF 786 787*---------------------------------------------------------------------* 788* calculate the derivatives of the twoelectron integrals and write 789* them to file: 790*---------------------------------------------------------------------* 791 K2INT = KLAST 792 L2INT = LWORK - K2INT + 1 793 794 IF ( LABELOP(1:5).EQ.'1DHAM' ) THEN 795 ITYPE = 1 796 ! open file for derivative integrals here... 797 LU2DER = -1 798 CALL GPOPEN(LU2DER,'AO2DRINT','UNKNOWN','SEQUENTIAL', 799 & 'UNFORMATTED',IDUMMY,.FALSE.) 800 REWIND (LU2DER) 801 ELSE IF ( LABELOP(1:5).EQ.'dh/dB' ) THEN 802 ITYPE = 5 803 ! file for derivative integrals will be opened in DR2WRT 804C LU2DER = -1 805C CALL GPOPEN(LU2DER,'AO2MGINT','UNKNOWN','SEQUENTIAL', 806C & 'UNFORMATTED',IDUMMY,.FALSE.) 807C REWIND (LU2DER) 808 ELSE 809 WRITE (LUPRI,*) 'Unknown operator label in CCGETH2D.' 810 CALL QUIT('Unknown operator label in CCGETH2D.') 811 END IF 812 813 IF (LOCDBG) THEN 814 WRITE (LUPRI,*) 'CCGETH2D> LABELOP,ITYPE,IATOM:', 815 & LABELOP(1:5),ITYPE,IATOM 816 END IF 817 818 MAXDIS = 1 819 I2TYP = 0 820 RELCAL = .FALSE. 821 TKTIME = .FALSE. 822 823 IF (LOCDBG) WRITE(LUPRI,*) 'CALL NOW TWOINT...' 824 825 CALL TWOINT(WORK(K2INT),L2INT,DUMMY, 826 & FMAT, DMAT, NDMAT, ISYMDM, IFCTYP, 827 & DUMMY,IDUMMY,NUMDIS,MAXDIS, 828 & ITYPE,MAXDIF,IATOM,NODV,NOPV,NOCONT, 829 & TKTIME,IPRHER,IPRNTA,IPRNTB,IPRNTC,IPRNTD, 830 & RETUR,IDUMMY,I2TYP,WORK(KJSTRS), 831 & WORK(KNPRIM),WORK(KNCONT),WORK(KIORBS), 832 & IDUMMY,IDUMMY,DUMMY,DUMMY,DUMMY, 833 & DUMMY,RELCAL,.false.) 834 835 TIM1 = SECOND() 836 837 IF (LOCDBG) WRITE(LUPRI,*) 'RETURNED FROM TWOINT...' 838 839 840*---------------------------------------------------------------------* 841* sort the derivative integrals: 842*---------------------------------------------------------------------* 843 TIM2 = SECOND() 844 845 IF (MAXDIF .EQ. 1) THEN 846 847 IF (LOCDBG) WRITE (LUPRI,*) 'Setting up LDERINT array:' 848 DO ICOOR = 1, 3 849 DO ICORSY = 1, NSYM 850 IF ( NDSINT(ICOOR,ICORSY-1) .GT. 0 ) THEN 851 LDERINT(ICORSY,ICOOR) = .TRUE. 852 ELSE 853 LDERINT(ICORSY,ICOOR) = .FALSE. 854 END IF 855 IF (LOCDBG) WRITE (LUPRI,'(2X,3I5,L5)') ICOOR,ICORSY, 856 & NDSINT(ICOOR,ICORSY-1),LDERINT(ICORSY,ICOOR) 857 END DO 858 END DO 859 860 ELSE 861 CALL QUIT('MAXDIF <> 1 not implemented in CCGETH2D.') 862 END IF 863 864 IF (LOCDBG) CALL FLSHFO(LUPRI) 865 866 CALL CCSORTDERAO(WORK,LWORK,3,LDERINT,ITYPE,IPRHER) 867 868 TIM3 = SECOND() 869 870*---------------------------------------------------------------------* 871* print timing & return: 872*---------------------------------------------------------------------* 873 IF (IPRINT.GT.1 .OR. LOCDBG) THEN 874 WRITE (LUPRI,'(/A,A,F12.2," seconds.")') 875 & ' Time used in TWO', 876 & 'INT: ', TIM1 - TIM0 877 WRITE (LUPRI,'(A,A,F12.2," seconds.")') 878 & ' Time used in CCS', 879 & 'SORTDERAO: ', TIM3 - TIM2 880 WRITE (LUPRI,'(/A,A,F12.2," seconds.")') 881 & ' Total time used ', 882 & 'in CCGETH2D :', SECOND() - TIM0 883 END IF 884 885 CALL FLSHFO(LUPRI) 886 887 CALL QEXIT('CCGETH2D') 888 RETURN 889 END 890*=====================================================================* 891* END OF SUBROUTINE CCGETH2D * 892*=====================================================================* 893