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*---------------------------------------------------------------------* 20c/* Deck CC_GMATRIX */ 21*=====================================================================* 22 SUBROUTINE CC_GMATRIX( LISTA, ! inp: A Zeta vector list 23 & LISTB, ! inp: B amplitude list 24 & LISTC, ! inp: C amplitude list 25 & LISTD, ! inp: D amplitude list 26 & NGTRAN, ! inp: nb. of G matrix transf. 27 & MXVEC, ! inp: max nb. of dot products 28 & IGTRAN, ! inp: indeces of G mat. transf. 29 & IGDOTS, ! inp: indeces of dot products 30 & GCON, ! inp: array for dot products 31 & WORK, ! scr: work space 32 & LWORK, ! inp: length of work space 33 & IOPTRES)! inp: output option 34*---------------------------------------------------------------------* 35* 36* Purpose: batched loop over G matrix transformations 37* 38* for more detailed documentation see: CC_GMAT 39* 40* Written by Christof Haettig, April 1999, based on CC_FMATRIX. 41* 42*=====================================================================* 43#if defined (IMPLICIT_NONE) 44 IMPLICIT NONE 45#else 46# include "implicit.h" 47#endif 48#include "priunit.h" 49#include "maxorb.h" 50#include "ccslvinf.h" 51#include "qm3.h" 52 53 LOGICAL LOCDBG 54 PARAMETER (LOCDBG = .FALSE.) 55 56 INTEGER MAXGTRAN 57 PARAMETER (MAXGTRAN = 100) 58 59 CHARACTER*(*) LISTA, LISTB, LISTC, LISTD 60 INTEGER IOPTRES, NGTRAN, MXVEC, LWORK 61 INTEGER IGTRAN(4,NGTRAN) 62 INTEGER IGDOTS(MXVEC,NGTRAN) 63 64#if defined (SYS_CRAY) 65 REAL WORK(LWORK) 66 REAL GCON(MXVEC,NGTRAN) 67#else 68 DOUBLE PRECISION WORK(LWORK) 69 DOUBLE PRECISION GCON(MXVEC,NGTRAN) 70#endif 71 72 INTEGER NTRAN, ISTART, IBATCH, NBATCH 73 INTEGER KTIMTRN, KDELTA1, KT1AMPB, KLAMDHB, KLAMDPB, KT1AMPC 74 INTEGER KLAMDHC, KLAMDPC, KCTR1, KCTR2, KFOCKB, KFOCKC 75 INTEGER KFBCOO, KFBCVV 76 INTEGER KEND, LEND 77 78 CALL QENTER('CC_GMATRIX') 79 80 NBATCH = (NGTRAN+MAXGTRAN-1)/MAXGTRAN 81 82 IF (LOCDBG) THEN 83 WRITE (LUPRI,*) 'Batching over G matrix transformations:' 84 WRITE (LUPRI,*) 'nb. of batches needed:', NBATCH 85 END IF 86 87 DO IBATCH = 1, NBATCH 88 ISTART = (IBATCH-1) * MAXGTRAN + 1 89 NTRAN = MIN(NGTRAN-(ISTART-1),MAXGTRAN) 90 91 KTIMTRN = 1 92 KDELTA1 = KTIMTRN + NTRAN 93 KT1AMPB = KDELTA1 + NTRAN 94 KLAMDHB = KT1AMPB + NTRAN 95 KLAMDPB = KLAMDHB + NTRAN 96 KT1AMPC = KLAMDPB + NTRAN 97 KLAMDHC = KT1AMPC + NTRAN 98 KLAMDPC = KLAMDHC + NTRAN 99 KCTR1 = KLAMDPC + NTRAN 100 KCTR2 = KCTR1 + NTRAN 101 KFOCKB = KCTR2 + NTRAN 102 KFOCKC = KFOCKB + NTRAN 103 KFBCOO = KFOCKC + NTRAN 104 KFBCVV = KFBCOO + NTRAN 105 KEND = KFBCVV + NTRAN 106 LEND = LWORK - KEND 107 108 IF (LEND .LT. 0) THEN 109 WRITE (LUPRI,*) 'Insufficient work space in CC_GMATRIX.' 110 WRITE (LUPRI,*) 'Available :',LWORK,' words,' 111 WRITE (LUPRI,*) 'Need at least:',KEND, ' words.' 112 CALL QUIT('Insufficient work space in CC_GMATRIX.') 113 END IF 114 115 IF (LOCDBG) THEN 116 WRITE (LUPRI,*) 'Batch No.:',IBATCH 117 WRITE (LUPRI,*) 'start at :',ISTART 118 WRITE (LUPRI,*) '# transf.:',NTRAN 119 END IF 120 121 CALL CC_GMAT( LISTA, LISTB, LISTC, LISTD, NTRAN, MXVEC, 122 & IGTRAN(1,ISTART),IGDOTS(1,ISTART),GCON(1,ISTART), 123 & WORK(KTIMTRN), WORK(KDELTA1), WORK(KT1AMPB), 124 & WORK(KLAMDHB), WORK(KLAMDPB), WORK(KT1AMPC), 125 & WORK(KLAMDHC), WORK(KLAMDPC), WORK(KCTR1), 126 & WORK(KCTR2), WORK(KFOCKB), WORK(KFOCKC), 127 & WORK(KFBCOO), WORK(KFBCVV), 128 & WORK(KEND), LEND, IOPTRES) 129 130 END DO 131 132 CALL QEXIT('CC_GMATRIX') 133 134 RETURN 135 END 136 137*---------------------------------------------------------------------* 138* END OF SUBROUTINE CC_GMATRIX * 139*---------------------------------------------------------------------* 140*---------------------------------------------------------------------* 141c/* Deck CC_GMAT */ 142*=====================================================================* 143 SUBROUTINE CC_GMAT( LISTA, ! inp: A Zeta vector list 144 & LISTB, ! inp: B amplitude list 145 & LISTC, ! inp: C amplitude list 146 & LISTD, ! inp: D amplitude list 147 & NGTRAN, ! inp: nb. of G matrix transf. 148 & MXVEC, ! inp: max nb. of dot products 149 & IGTRAN, ! inp: index array of G mat. transf. 150 & IGDOTS, ! inp: index array of dot products 151 & GCON, ! inp: array for dot products 152 & TIMTRN, ! scr: 153 & KDELTA1,! scr: 154 & KT1AMPB,! scr: 155 & KLAMDHB,! scr: 156 & KLAMDPB,! scr: 157 & KT1AMPC,! scr: 158 & KLAMDHC,! scr: 159 & KLAMDPC,! scr: 160 & KCTR1, ! scr: 161 & KCTR2, ! scr: 162 & KFOCKB, ! scr: 163 & KFOCKC, ! scr: 164 & KFBCOO, ! scr: 165 & KFBCVV, ! scr: 166 & WORK, ! scr: work space 167 & LWORK, ! inp: length of work space 168 & IOPTRES)! inp: output option 169*---------------------------------------------------------------------* 170* 171* Purpose: AO-direct calculation of a linear transformation of two 172* CC amplitude vectors, T^B and T^C, with the CC G matrix 173* (second partial derivative of the CC lagrangian) 174* 175* index combinations for L^A, T^B, T^C passed in IGTRAN, 176* indeces for T^D in IGDOTS 177* 178* return of the result vectors: 179* 180* IOPTRES = 0 : all result vectors are written to a direct 181* access file, LISTD is used as file name, 182* the start addresses of the vectors are 183* returned in IGTRAN(4,*) 184* 185* IOPTRES = 1 : the vectors are kept and returned in WORK 186* if possible, start addresses returned in 187* IGTRAN(4,*). N.B.: if WORK is not large 188* enough, IOPTRES is automatically reset to 0 189* 190* IOPTRES = 3 : each result vector is written to its own 191* file by a call to CC_WRRSP, LISTD is used 192* as list type and IGTRAN(4,*) as list index 193* NOTE that IGTRAN(4,*) is in this case input! 194* 195* IOPTRES = 4 : each result vector is added to a vector on 196* file by a call to CC_WARSP, LISTD is used 197* as list type and IGTRAN(4,*) as list index 198* NOTE that IGTRAN(4,*) is in this case input! 199* 200* IOPTRES = 5 : the result vectors are dotted on a array 201* of vectors, the type of the arrays given 202* by LISTD and the indeces from IGDOTS 203* the result of the dot products is returned 204* in the GCON array 205* 206* 207* symmetries/variables: 208* 209* ISYRES : result vector 210* ISYCTR : lagrangian multipliers 211* ISYMTB : B response amplitudes 212* ISYMTC : C response amplitudes 213* 214* DELTA1 : single excitation part (total sum) 215* 216* Note: it is assumed, that the integrals (ia|jb) are on disk 217* 218* used vectors of the size (1/2) V^2 O^2 : 219* CCS: (ia|jb) 220* CC2: (ia|jb), CTR2 221* CCSD: (ia|jb), CTR2, TAmpB2, TAmpC2, DELTA2 222* 223* Written by Christof Haettig, Aug. 1996; restructed, Sept. 1998. 224* 225* included CC-R12: Christian Neiss, november 2005 226* 227*=====================================================================* 228 USE PELIB_INTERFACE, ONLY: USE_PELIB, PELIB_IFC_QRTRANSFORMER 229#if defined (IMPLICIT_NONE) 230 IMPLICIT NONE 231#else 232# include "implicit.h" 233#endif 234#include "priunit.h" 235#include "iratdef.h" 236#include "cbieri.h" 237#include "mxcent.h" 238#include "eritap.h" 239#include "maxorb.h" 240#include "distcl.h" 241#include "ccorb.h" 242#include "ccisao.h" 243#include "ccsdsym.h" 244#include "ccsdinp.h" 245#include "aovec.h" 246#include "blocks.h" 247#include "second.h" 248#include "ccnoddy.h" 249#include "r12int.h" 250#include "ccr12int.h" 251#include "ccsections.h" 252#include "ccslvinf.h" 253#include "qm3.h" 254 255* local parameters: 256 CHARACTER MSGDBG*(16) 257 PARAMETER (MSGDBG='[debug] CC_GMAT>') 258#if defined (SYS_CRAY) 259 REAL ZERO, ONE, TWO, THREE 260#else 261 DOUBLE PRECISION ZERO, ONE, TWO, THREE 262#endif 263 PARAMETER (ZERO = 0.0d0, ONE = 1.0d0, TWO = 2.0d0, THREE = 3.0d0) 264 LOGICAL LOCDBG 265 PARAMETER (LOCDBG = .FALSE.) 266 LOGICAL LDOUBL 267 PARAMETER (LDOUBL = .TRUE.) 268 INTEGER ISYM0 269 PARAMETER (ISYM0 = 1) 270 INTEGER LUGMAT 271 272 CHARACTER*(*) LISTA, LISTB, LISTC, LISTD 273 CHARACTER*6 FILGMA 274 INTEGER MXVEC, NGTRAN, IOPTRES, LWORK 275 INTEGER IGTRAN(4,NGTRAN) 276 INTEGER IGDOTS(MXVEC,NGTRAN) 277 278 CHARACTER*(1) RSPTYP 279 CHARACTER*(8) FNTOC 280 CHARACTER*(7) FN3FOP2X 281 CHARACTER*(6) FNCKJD, FN3VI, FN3FOPX 282 CHARACTER*(5) FNDKBC3 283 INTEGER LUTOC,LU3FOP2X,LUCKJD, LU3VI, LU3FOPX,LUDKBC3 284 PARAMETER (FNTOC ='CCSDT_OC',FN3FOP2X='PTFOP2X', 285 * FNCKJD='CKJDEL',FN3VI ='CC3_VI' , 286 * FN3FOPX='PTFOPX' ,FNDKBC3='DKBC3') 287 288 289#if defined (SYS_CRAY) 290 REAL GCON(MXVEC,NGTRAN) 291 REAL WORK(LWORK) 292 REAL DDOT, FACT, XNORM, TIMTRN(NGTRAN) 293 REAL DTIME, TIMALL, TIMFCK, TIMIO, TIMCCS, TIMRDAO 294 REAL TIM21CD, TIM21SD, TIMDBL, TIMINT, TIMTRBT,CONVRT 295 REAL TIMTRIP 296#else 297 DOUBLE PRECISION GCON(MXVEC,NGTRAN) 298 DOUBLE PRECISION WORK(LWORK) 299 DOUBLE PRECISION DDOT, FACT, XNORM, TIMTRN(NGTRAN) 300 DOUBLE PRECISION DTIME, TIMALL, TIMFCK, TIMIO, TIMCCS, TIMRDAO 301 DOUBLE PRECISION TIM21CD, TIM21SD, TIMDBL, TIMINT, TIMTRBT,CONVRT 302 DOUBLE PRECISION TIMTRIP 303#endif 304 305 INTEGER KT1AMP0, KT1AMPB(NGTRAN), KT1AMPC(NGTRAN), KCTR1(NGTRAN) 306 INTEGER KLAMDP0, KLAMDPB(NGTRAN), KLAMDPC(NGTRAN) 307 INTEGER KLAMDH0, KLAMDHB(NGTRAN), KLAMDHC(NGTRAN) 308 INTEGER KFOCKB(NGTRAN), KFOCKC(NGTRAN) 309 INTEGER KFBCOO(NGTRAN), KFBCVV(NGTRAN) 310 INTEGER KDELTA1(NGTRAN), KCTR2(NGTRAN) 311 312 CHARACTER MODEL*(10), MODELW*(10), APROXR12*(3) 313 LOGICAL LSAME, LLSAME, CCSDR12 314 INTEGER INDEXA(MXCORB_CC) 315 INTEGER IDLSTA, IDLSTB, IDLSTC, IDLSTD 316 INTEGER ISYRES, ISYCTR, ISYMTB, ISYMTC, ISYMTD 317 INTEGER KEND0, KEND1, KEND1A, KEND2, KFREE, KENDSV, KENDE3 318 INTEGER LEND0, LEND1, LEND1A, LEND2, LFREE, LWRKSV, LENDE3 319 INTEGER KENDB2, KENDB3, KENDA2, KENDA3, KEND3, KXIAJB,KC4O,KX4O 320 INTEGER LENDB2, LENDB3, LENDA2, LENDA3, LEND3, KDELTA2, KCDTERM 321 INTEGER ISYTCTB, ISYOVOV, ISYMFB, ISYMFC, ISYFBC, ISYDEL 322 INTEGER KLIAJB, KSCR, LSCR, LWRK0, ITRAN, ISTRT, IEND, IFILE 323 INTEGER KTTZET, KTZBOO, KTZBVV, KXINT, KDUM, KZETA2, IOPTW 324 INTEGER KT2AMPB, KT2AMPC, KEMAT1, KEMAT2 325 INTEGER KDELTBC, KDELTCB, KINDXB, KCCFB1 326 INTEGER KODCL1, KODCL2, KODBC1, KODBC2, KRDBC1, KRDBC2 327 INTEGER KODPP1, KODPP2, KRDPP1, KRDPP2, KRECNR 328 INTEGER NTOSYM,ISYMD1,NTOT,ILLL,NUMDIS,IDEL,IDEL2, LENALL,IVEC 329 INTEGER ISYDIS,ISYC4O,ISYX4O,IOPT,ADR,IERR,IADRDEL,LEN1,LEN2 330 INTEGER ISYMI, ISYMJ, ISYMA, ISYMB, NAB, NBA, NIJ, NJI 331 INTEGER KDSRHF0,KDSRHFBC,ISYMXB,ISYMXC,ISYCTB,ISYCTC,IOPTWE 332 INTEGER ITAIKJ(8,8), NTAIKJ(8), ICOUNT, ISYM, ISYMAIK 333 INTEGER KXBIAJK, KEBIAJK, KZBIAJK, KCBIAJK, KXLMDHB, KXLMDPB 334 INTEGER KXCIAJK, KECIAJK, KZCIAJK, KCCIAJK, KXLMDHC, KXLMDPC 335 INTEGER MDISAO,MDSRHF,M2BST,MT2BGD,MT2SQ,MT2AM,MSCRATCH,MFREE 336 337 INTEGER ILSTSYM, IOPTTCME, IOPTWR12, LENMOD 338 INTEGER IDUM 339 INTEGER KDELTAR12, KR12SCR, KVABKL, LUNIT, LENR12 340 341 CALL QENTER('CC_GMAT') 342 343*---------------------------------------------------------------------* 344* begin: 345*---------------------------------------------------------------------* 346 IF (LOCDBG) THEN 347 Call AROUND('ENTERED CC_GMAT') 348 IF (DIRECT) WRITE(LUPRI,'(/1X,A)') 'AO direct transformation' 349 CALL FLSHFO(LUPRI) 350 END IF 351 352 IF ( .NOT. (CCS .OR. CC2 .OR. CCSD .OR. CC3) ) THEN 353 WRITE(LUPRI,'(/1x,a)') 'CC_GMAT called for a Coupled Cluster ' 354 & //'method not implemented in CC_GMAT...' 355 CALL QUIT('Unknown CC method in CC_GMAT.') 356 END IF 357 358 IF (LISTB(1:1).NE.'R' .OR. LISTC(1:1).NE.'R') THEN 359 WRITE(LUPRI,'(2A)') 'LISTB and LISTC must refer to', 360 & ' t-amplitude vectors in CC_GMAT.' 361 CALL QUIT('Illegal LISTB or LISTC in CC_GMAT.') 362 END IF 363 364 !set CCSDR12 365 CCSDR12 = CCR12 .AND. .NOT.(CC2.OR.CCS) 366 367 IF (IPRINT.GT.0) THEN 368 369 WRITE (LUPRI,'(//1X,A1,50("="),A1)')'+','+' 370 371 IF (LISTA(1:2).EQ.'L0') THEN 372 WRITE (LUPRI,'(1x,A52)') 373 & '| G MATRIX TRANSFORMATION SECTION |' 374 ELSE 375 WRITE (LUPRI,'(1X,A52)') 376 & '| GENERALIZED G MATRIX TRANSFORMATION SECTION |' 377 END IF 378 379 IF (IOPTRES.EQ.3) THEN 380 WRITE (LUPRI,'(1X,A52)') 381 & '| (result is written to file) |' 382 ELSE IF (IOPTRES.EQ.4) THEN 383 WRITE (LUPRI,'(1X,A52)') 384 & '| (result is added to a vector on file) |' 385 ELSE IF (IOPTRES.EQ.5) THEN 386 WRITE (LUPRI,'(1X,A52)') 387 & '| (result used to calculate dot products) |' 388 END IF 389 390 IF (IOPTRES.EQ.5) THEN 391 WRITE (LUPRI,'(1X,A1,50("-"),A1)') '+','+' 392 WRITE (LUPRI,'(1X,A52)') 393 & '| L vec. | R vec. | R vec. | # prod. | |' 394 WRITE (LUPRI,'(1X,4(A,A3),A)') '| ',LISTA(1:3), ' No.| ', 395 & LISTB(1:3),' No.| ',LISTC(1:3),' No.| with ', 396 & LISTD(1:3),' | time/secs |' 397 ELSE 398 WRITE (LUPRI,'(1X,A1,50("-"),A1)') '+','+' 399 WRITE (LUPRI,'(1X,A52)') 400 & '| L vec. | R vec. | R vec. | result | |' 401 WRITE (LUPRI,'(1X,A2,A3,3(A,A3),A)') '| ',LISTA(1:3), 402 & ' No.| ',LISTB(1:3),' No.| ',LISTC(1:3),' No.| ', 403 & LISTD(1:3),' No. | time/secs |' 404 END IF 405 406 WRITE (LUPRI,'(1X,A1,50("-"),A1)') '+','+' 407 408 END IF 409 410*---------------------------------------------------------------------* 411* initialize timings: 412*---------------------------------------------------------------------* 413 TIMALL = SECOND() 414 TIMIO = ZERO 415 TIMFCK = ZERO 416 TIMCCS = ZERO 417 TIMINT = ZERO 418 TIMRDAO = ZERO 419 TIM21CD = ZERO 420 TIM21SD = ZERO 421 TIMDBL = ZERO 422 TIMTRBT = ZERO 423 TIMTRIP = ZERO 424 425*---------------------------------------------------------------------* 426* flush print unit, set dummy work space address to minus a large numb. 427*---------------------------------------------------------------------* 428 Call FLSHFO(LUPRI) 429 430 KDUM = -99 999 999 ! dummy address 431 432 IF (LOCDBG) THEN 433 WRITE(LUPRI,'(/1x,a,i15)') 'work space in CC_GMAT:',LWORK 434 Call FLSHFO(LUPRI) 435 END IF 436 437*---------------------------------------------------------------------* 438* if IOPTRES=0,1 open direct access file for result vectors: 439* if IOPTRES=5 initialize GCON array 440*---------------------------------------------------------------------* 441 IF ( IOPTRES.EQ.0 .OR. IOPTRES.EQ.1 ) THEN 442 443 FILGMA = LISTD(1:6) 444 445 LUGMAT = -1 446 CALL WOPEN2(LUGMAT,FILGMA,64,0) 447 IADRDEL = 1 448 449 ELSE IF ( IOPTRES.EQ.3 .OR. IOPTRES.EQ.4 ) THEN 450 CONTINUE 451 ELSE IF ( IOPTRES.EQ.5 ) THEN 452 IF (MXVEC*NGTRAN.NE.0) CALL DZERO(GCON,MXVEC*NGTRAN) 453 ELSE 454 CALL QUIT('Illegal value of IOPTRES in CC_GMAT.') 455 END IF 456 457*---------------------------------------------------------------------* 458* set special symmetry array for one-index transformed integrals 459* and several maximum dimensions: 460*---------------------------------------------------------------------* 461 MDISAO = 0 462 MDSRHF = 0 463 M2BST = 0 464 MT2BGD = 0 465 MT2SQ = 0 466 MT2AM = 0 467 DO ISYM = 1, NSYM 468 ICOUNT = 0 469 DO ISYMAIK = 1, NSYM 470 ISYMJ = MULD2H(ISYMAIK,ISYM) 471 ITAIKJ(ISYMAIK,ISYMJ) = ICOUNT 472 ICOUNT = ICOUNT + NT2BCD(ISYMAIK)*NRHF(ISYMJ) 473 END DO 474 NTAIKJ(ISYM) = ICOUNT 475 MDISAO = MAX(MDISAO,NDISAO(ISYM)) 476 MDSRHF = MAX(MDSRHF,NDSRHF(ISYM)) 477 M2BST = MAX(M2BST ,N2BST(ISYM) ) 478 MT2BGD = MAX(MT2BGD,NT2BGD(ISYM)) 479 MT2SQ = MAX(MT2SQ ,NT2SQ(ISYM) ) 480 MT2AM = MAX(MT2AM ,NT2AM(ISYM) ) 481 END DO 482 483*---------------------------------------------------------------------* 484* allocate work space for (ia|jb) at the very end of the work space 485* and read them in and calculate L(ia|jb) in place: 486*---------------------------------------------------------------------* 487 ISYOVOV = ISYM0 488 KLIAJB = LWORK - NT2AM(ISYM0) 489 LWRK0 = KLIAJB - 1 490 491 IF (LWRK0 .LT. 0) THEN 492 CALL QUIT('Insufficient work space in CC_GMAT. (0)') 493 END IF 494 495 Call CCG_RDIAJB(WORK(KLIAJB),NT2AM(ISYOVOV)) 496 497 IOPTTCME = 1 498 Call CCSD_TCMEPK(WORK(KLIAJB),ONE,ISYOVOV,IOPTTCME) 499 500*---------------------------------------------------------------------* 501* allocate work space for zeroth-order amplitudes and Lambda matrices, 502* read amplitudes and calculate Lambda matrices: 503*---------------------------------------------------------------------* 504 KT1AMP0 = 1 505 KLAMDP0 = KT1AMP0 + NT1AM(ISYM0) 506 KLAMDH0 = KLAMDP0 + NGLMDT(ISYM0) 507 KEND1 = KLAMDH0 + NGLMDT(ISYM0) 508 LEND1 = LWRK0 - KEND1 509 510 IF (LEND1 .LT. 0) THEN 511 CALL QUIT('Insufficient work space in CC_GMAT. (1a)') 512 END IF 513 514* read zeroth-order amplitudes: 515 IOPT = 1 516 Call CC_RDRSP('R0',0,ISYM0,IOPT,MODEL,WORK(KT1AMP0),WORK(KDUM)) 517 518* get zeroth-order Lambda matrices: 519 Call LAMMAT(WORK(KLAMDP0),WORK(KLAMDH0),WORK(KT1AMP0), 520 & WORK(KEND1),LEND1) 521 522*=====================================================================* 523* global intermediate & CCS section: 524* allocate work space for global two-index intermediates 525* and calculate them in a loop over transformations, 526* then calculate the `CCS' contributions to the result vector 527*=====================================================================* 528 529 DO ITRAN = 1, NGTRAN 530 531 TIMTRN(ITRAN) = - SECOND() 532C 533C ------------------------------------------ 534C set vector indeces and several symmetries: 535C ------------------------------------------ 536C 537 IDLSTA = IGTRAN(1,ITRAN) 538 IDLSTB = IGTRAN(2,ITRAN) 539 IDLSTC = IGTRAN(3,ITRAN) 540 541 ISYCTR = ILSTSYM(LISTA,IDLSTA) 542 ISYMTB = ILSTSYM(LISTB,IDLSTB) 543 ISYMTC = ILSTSYM(LISTC,IDLSTC) 544 545 ISYTCTB = MULD2H(ISYMTB,ISYMTC) ! TAmpB x TAmpC 546 ISYRES = MULD2H(ISYTCTB,ISYCTR) ! result vector 547 548 ISYMFB = MULD2H(ISYMTB,ISYOVOV) 549 ISYMFC = MULD2H(ISYMTC,ISYOVOV) 550 ISYFBC = MULD2H(ISYTCTB,ISYOVOV) 551C 552 IF (LOCDBG) THEN 553 WRITE (LUPRI,*) LISTA,LISTB,LISTC,LISTD 554 WRITE (LUPRI,*) IDLSTA,IDLSTB,IDLSTC,ITRAN 555 WRITE (LUPRI,*) ISYCTR,ISYMTB,ISYMTC,ISYRES,ISYTCTB 556 END IF 557C 558C ------------------------------------------------------------ 559C allocate single parts of the vectors and Fock intermediates: 560C ------------------------------------------------------------ 561C 562 KDELTA1(ITRAN) = KEND1 563 KT1AMPB(ITRAN) = KDELTA1(ITRAN) + NT1AM(ISYRES) 564 KLAMDHB(ITRAN) = KT1AMPB(ITRAN) + NT1AM(ISYMTB) 565 KLAMDPB(ITRAN) = KLAMDHB(ITRAN) + NGLMDT(ISYMTB) 566 KT1AMPC(ITRAN) = KLAMDPB(ITRAN) + NGLMDT(ISYMTB) 567 KLAMDHC(ITRAN) = KT1AMPC(ITRAN) + NT1AM(ISYMTC) 568 KLAMDPC(ITRAN) = KLAMDHC(ITRAN) + NGLMDT(ISYMTC) 569 KCTR1(ITRAN) = KLAMDPC(ITRAN) + NGLMDT(ISYMTC) 570 KFOCKB(ITRAN) = KCTR1(ITRAN) + NT1AM(ISYCTR) 571 KFOCKC(ITRAN) = KFOCKB(ITRAN) + NT1AM(ISYMFB) 572 KFBCOO(ITRAN) = KFOCKC(ITRAN) + NT1AM(ISYMFC) 573 KFBCVV(ITRAN) = KFBCOO(ITRAN) + NMATIJ(ISYFBC) 574 KEND1 = KFBCVV(ITRAN) + NMATAB(ISYFBC) 575 576 LEND1 = LWRK0 - KEND1 577 578* allocate scratch vector needed in Fock matrices calculations 579* and for the first CCS term: 580 KSCR = KEND1 581 LSCR = MAX(NMATIJ(ISYFBC),NMATAB(ISYFBC),NT1AM(ISYRES)) 582 583 IF (LEND1 .LT. LSCR) THEN 584 CALL QUIT('Insufficient work space in CC_GMAT. (1b)') 585 END IF 586C 587C --------------------------------------------- 588C initialize single parts of the result vector: 589C --------------------------------------------- 590C 591* initialize singles excitation result vector: 592 Call DZERO(WORK(KDELTA1(ITRAN)),NT1AM(ISYRES)) 593 594 DTIME = SECOND() 595C 596C -------------------------------------------------------------- 597C for CCS and zeroth order zeta vector all contributions vanish: 598C ---> we skip the rest of the loop 599C -------------------------------------------------------------- 600C 601 IF (CCS .AND. LISTA(1:2).EQ.'L0') GOTO 666 602C 603* read A response zeta vector: 604 IOPT = 1 605 Call CC_RDRSP(LISTA,IDLSTA,ISYCTR,IOPT,MODEL, 606 & WORK(KCTR1(ITRAN)),WORK(KDUM)) 607 608* read B response amplitudes: 609 IOPT = 1 610 Call CC_RDRSP(LISTB,IDLSTB,ISYMTB,IOPT,MODEL, 611 & WORK(KT1AMPB(ITRAN)),WORK(KDUM)) 612 613* read C response amplitudes: 614 IOPT = 1 615 Call CC_RDRSP(LISTC,IDLSTC,ISYMTC,IOPT,MODEL, 616 & WORK(KT1AMPC(ITRAN)),WORK(KDUM)) 617 618 TIMIO = TIMIO + SECOND() - DTIME 619 620 IF (LOCDBG) THEN 621 Call AROUND('zeroth order T1 amplitudes:') 622 Call CC_PRP(WORK(KT1AMP0),WORK(KDUM),ISYM0,1,0) 623 624 Call AROUND('ZETA1 multipliers:') 625 WRITE (LUPRI,*) 626 * 'List, Index, Symmetry:',LISTA, IDLSTA, ISYCTR 627 Call CC_PRP(WORK(KCTR1(ITRAN)),WORK(KDUM),ISYCTR,1,0) 628 629 Call AROUND('response T1 amplitudes B:') 630 WRITE (LUPRI,*) 631 * 'List, Index, Symmetry:',LISTB, IDLSTB, ISYMTB 632 Call CC_PRP(WORK(KT1AMPB(ITRAN)),WORK(KDUM),ISYMTB,1,0) 633 634 Call AROUND('response T1 amplitudes C:') 635 WRITE (LUPRI,*) 636 * 'List, Index, Symmetry:',LISTC, IDLSTC, ISYMTC 637 Call CC_PRP(WORK(KT1AMPC(ITRAN)),WORK(KDUM),ISYMTC,1,0) 638 639 Call FLSHFO(LUPRI) 640 END IF 641C 642C ----------------------------------- 643C calculate response Lambda matrices: 644C ----------------------------------- 645C 646* calculate B response Lambda matrices: 647 Call CCLR_LAMTRA(WORK(KLAMDP0),WORK(KLAMDPB(ITRAN)), 648 & WORK(KLAMDH0),WORK(KLAMDHB(ITRAN)), 649 & WORK(KT1AMPB(ITRAN)),ISYMTB) 650 651* calculate C response Lambda matrices: 652 Call CCLR_LAMTRA(WORK(KLAMDP0),WORK(KLAMDPC(ITRAN)), 653 & WORK(KLAMDH0),WORK(KLAMDHC(ITRAN)), 654 & WORK(KT1AMPC(ITRAN)),ISYMTC) 655C 656C ------------------------------------------------------------- 657C calculate occupied/virtual part of first order Fock matrices: 658C ------------------------------------------------------------- 659C 660 DTIME = SECOND() 661 662* calculate tF^B_ia: 663 Call CCG_LXD(WORK(KFOCKB(ITRAN)), ISYMFB, 664 & WORK(KT1AMPB(ITRAN)),ISYMTB, 665 & WORK(KLIAJB), ISYOVOV,0) 666 667* calculate tF^C_ia: 668 Call CCG_LXD(WORK(KFOCKC(ITRAN)), ISYMFC, 669 & WORK(KT1AMPC(ITRAN)),ISYMTC, 670 & WORK(KLIAJB), ISYOVOV,0) 671 672 IF (LOCDBG) THEN 673 Call AROUND(' one-index B transformed Fock marix ') 674 Call CC_PRP(WORK(KFOCKB(ITRAN)),WORK(KDUM),ISYMFB,1,0) 675 676 Call AROUND(' one-index C transformed Fock marix ') 677 Call CC_PRP(WORK(KFOCKC(ITRAN)),WORK(KDUM),ISYMFC,1,0) 678 END IF 679C 680C ------------------------------------------------ 681C double one-index transformed Fock matrix ttF^BC: 682C ------------------------------------------------ 683C 684* calculate occ/occ block of double one-index 685* transformed Fock matrix: 686 Call CCG_1ITROO(WORK(KFBCOO(ITRAN)), ISYFBC, 687 & WORK(KFOCKB(ITRAN)), ISYMFB, 688 & WORK(KT1AMPC(ITRAN)),ISYMTC) 689 690 Call CCG_1ITROO(WORK(KSCR), ISYFBC, 691 & WORK(KFOCKC(ITRAN)), ISYMFC, 692 & WORK(KT1AMPB(ITRAN)),ISYMTB) 693 694 Call DAXPY(NMATIJ(ISYFBC),ONE,WORK(KSCR),1, 695 & WORK(KFBCOO(ITRAN)),1) 696 697* calculate vir/vir block of double one-index 698* transformed Fock matrix: 699 Call CCG_1ITRVV(WORK(KFBCVV(ITRAN)), ISYFBC, 700 & WORK(KFOCKB(ITRAN)), ISYMFB, 701 & WORK(KT1AMPC(ITRAN)),ISYMTC ) 702 703 Call CCG_1ITRVV(WORK(KSCR), ISYFBC, 704 & WORK(KFOCKC(ITRAN)), ISYMFC, 705 & WORK(KT1AMPB(ITRAN)),ISYMTB ) 706 707 Call DAXPY(NMATAB(ISYFBC),ONE,WORK(KSCR),1, 708 & WORK(KFBCVV(ITRAN)),1) 709 710 TIMFCK = TIMFCK + SECOND() - DTIME 711C 712C ------------------------------------------------- 713C CCS part: < Zeta_1 | [ttH^BC, tau_1] | HF> 714C ------------------------------------------------- 715C 716 DTIME = SECOND() 717 718* do one-index transformation with Zeta vector: 719 IOPT = 2 720 Call CCG_1ITRVO(WORK(KSCR),ISYRES, 721 & WORK(KFBCOO(ITRAN)),WORK(KFBCVV(ITRAN)), 722 & ISYTCTB,WORK(KCTR1(ITRAN)),ISYCTR,IOPT ) 723 724* add contribution to result vector: 725 Call DAXPY(NT1AM(ISYRES),ONE,WORK(KSCR),1, 726 & WORK(KDELTA1(ITRAN)),1) 727 728 IF (LOCDBG) THEN 729 WRITE (LUPRI,*) MSGDBG, 'ISYCTR,-RES,-TCTB:',ISYCTR,ISYRES, 730 & ISYTCTB 731 WRITE (LUPRI,*) MSGDBG, 'CTR1:' 732 Call CC_PRP(WORK(KCTR1(ITRAN)),WORK,ISYCTR,1,0) 733 WRITE (LUPRI,*) MSGDBG, 734 & 'double one-index trans. Fock mat. (o/o):' 735 WRITE (LUPRI,'(5D21.14)') 736 & (WORK(KFBCOO(ITRAN)-1+I),I=1,NMATIJ(ISYFBC)) 737 WRITE (LUPRI,*) MSGDBG, 738 & 'double one-index trans. Fock mat. (v/v):' 739 WRITE (LUPRI,'(5D21.14)') 740 & (WORK(KFBCVV(ITRAN)-1+I),I=1,NMATAB(ISYFBC)) 741 WRITE (LUPRI,*) MSGDBG, 742 & 'contribution of double trans. F to DELTA:' 743 Call CC_PRP(WORK(KSCR),WORK,ISYRES,1,0) 744 END IF 745C 746C ------------------------------------------------------- 747C second CCS contribution (involving L(ia,jb) integrals): 748C ------------------------------------------------------- 749C 750 ISYCTB = MULD2H(ISYCTR,ISYMTB) 751 752* allocate temporay memory for tZeta^B, ttZeta^BC matrices 753* and a scratch vector: 754 KTTZET = KEND1 755 KTZBOO = KTTZET + NT1AM(ISYRES) 756 KTZBVV = KTZBOO + NMATIJ(ISYCTB) 757 KSCR = KTZBVV + NMATAB(ISYCTB) 758 KEND2 = KSCR + NT1AM(ISYRES) 759 LEND2 = LWORK - KEND2 760 761 IF (LEND2 .LT. 0) THEN 762 CALL QUIT('Insufficient work space in CC_GMAT. (2)') 763 END IF 764 765* calculate one-index transformed Zeta matrix: 766 Call CCG_1ITROO(WORK(KTZBOO), ISYCTB, 767 & WORK(KCTR1(ITRAN)), ISYCTR, 768 & WORK(KT1AMPB(ITRAN)), ISYMTB ) 769 770 Call CCG_1ITRVV(WORK(KTZBVV), ISYCTB, 771 & WORK(KCTR1(ITRAN)), ISYCTR, 772 & WORK(KT1AMPB(ITRAN)), ISYMTB ) 773 774* calculate double one-index transformed Zeta matrix: 775 IOPT = 1 776 Call CCG_1ITRVO(WORK(KTTZET),ISYRES, 777 & WORK(KTZBOO),WORK(KTZBVV),ISYCTB, 778 & WORK(KT1AMPC(ITRAN)),ISYMTC,IOPT ) 779 780* contract with L(ia|jb): 781 Call CCG_LXD(WORK(KSCR), ISYRES, WORK(KTTZET), 782 & ISYRES, WORK(KLIAJB), ISYOVOV,0) 783 784 IF (LOCDBG) THEN 785 Call AROUND('double one-index transformed Zeta matrix:') 786 Call CC_PRP(WORK(KTTZET),WORK(KDUM),ISYRES,1,0) 787 Call AROUND('Contribution of L(ia|jb) to CCS part:') 788 Call CC_PRP(WORK(KSCR),WORK(KDUM),ISYRES,1,0) 789 END IF 790 791* add result to delta vector: 792 Call DAXPY(NT1AM(ISYRES),ONE,WORK(KSCR),1, 793 & WORK(KDELTA1(ITRAN)),1) 794 795 IF (LOCDBG) THEN 796 WRITE (LUPRI,*) MSGDBG, 797 & 'Delta vector at the end of the CCS part:' 798 Call CC_PRP(WORK(KDELTA1(ITRAN)),WORK(KDUM),ISYRES,1,0) 799 WRITE (LUPRI,*) MSGDBG, 'Norm of Delta1 after CCS part:', 800 & DSQRT(DDOT(NT1AM(ISYRES),WORK(KDELTA1(ITRAN)),1, 801 & WORK(KDELTA1(ITRAN)),1)) 802 END IF 803C 804 666 CONTINUE 805C 806 TIMCCS = TIMCCS + SECOND() - DTIME 807 808 TIMTRN(ITRAN) = TIMTRN(ITRAN) + SECOND() 809 810 END DO 811*---------------------------------------------------------------------* 812* end of CCS part; regain work space from L(ia|jb) 813* in memory now for each transformation: 814* 815* KT1AMP0: singles cluster amplitudes 816* KLAMDH0: zeroth-order Lambda hole matrix 817* KLAMDP0: zeroth-order Lambda particle matrix 818* KDELTA1: singles result vector 819* KT1AMPB: singles response amplitudes B 820* KLAMDHB: response Lambda hole matrix for B 821* KLAMDPB: response LAmbda particle matrix for B 822* KT1AMPC: singles response amplitudes C 823* KLAMDHC: response Lambda hole matrix for C 824* KLAMDPC: response LAmbda particle matrix for C 825* KCTR1 : singles response multipliers A 826* KFOCKB : Fock matrix tF^B(i a) 827* KFOCKC : Fock matrix tF^C(i a) 828* KFBCOO : Fock matrix ttF^BC(i j) 829* KFBCVV : Fock matrix ttF^BC(a b) 830* 831* KEND1 : free work space 832* 833*---------------------------------------------------------------------* 834 835 LEND1 = LWORK - KEND1 836 837 IF (CCS) GOTO 777 838 839*=====================================================================* 840* 841* CC2 part: C and D terms: <Zeta_2| [ttH^BC, tau_1] |HF> 842* 843* this part needs the AO integrals and the double excitation part 844* of the Lagrangian multipliers (Zeta_2) --> we make a compromise 845* between memory and CPU requirements: 846* 847* -- outermost batched loop over transformations with different Zeta_2 848* read all Zeta_2 vectors for this batch 849* -- loop over integral distributions 850* -- loop over transformations within this batch 851* calculate CCG_21CD contributions 852* -- end loop 853* -- end loop 854* -- end loop 855* 856* >> within the loop over AO integrals no I/O over vectors 857* >> integral recalculation scales only with nb. of Zeta_2 vectors 858* >> only one integral calculation for real G matrix: Zeta = 'L0' 859* 860*=====================================================================* 861 862*---------------------------------------------------------------------* 863* reserve memory for AO integrals and work space ERI & CCG_21CD: 864*---------------------------------------------------------------------* 865 MSCRATCH = MDISAO + 2*MDSRHF + 10*M2BST + 2*MT2BGD 866 & + MXPRIM*MXCONT + (8*MXSHEL*MXCONT + 1)/IRAT 867 MSCRATCH = MAX(MSCRATCH,MT2AM) 868 869 IF ( LEND1 .LT. MSCRATCH+MT2SQ ) THEN 870 CALL QUIT('Insufficient work space in CC_GMAT. (CC2)') 871 END IF 872 873*---------------------------------------------------------------------* 874* start batched loop over transformations with variable batch lengths: 875*---------------------------------------------------------------------* 876 IEND = 0 877 878 DO WHILE ( IEND.LT.NGTRAN ) 879C 880C -------------------------------------------------- 881C determine length of next batch and read the double 882C excitation parts of the Lagrangian multipliers: 883C (read as packed matrix, which is then squared up) 884C -------------------------------------------------- 885C 886 ISTRT = IEND + 1 887 888 IDLSTA = -999 889 KEND2 = KEND1 890 LEND2 = LEND1 891 MFREE = LEND1 - MSCRATCH 892 893 DO WHILE (MFREE.GT.MT2SQ .AND. IEND.LT.NGTRAN) 894 IEND = IEND + 1 895 IF ( IGTRAN(1,IEND) .EQ. IDLSTA ) THEN 896 KCTR2(IEND) = KCTR2(IEND-1) 897 ELSE 898 KCTR2(IEND) = KEND2 899 IDLSTA = IGTRAN(1,IEND) 900 ISYCTR = ILSTSYM(LISTA,IDLSTA) 901 KEND2 = KCTR2(IEND) + NT2SQ(ISYCTR) 902 LEND2 = LWORK - KEND2 903 MFREE = LEND2 - MSCRATCH 904 905 IOPT = 2 906 DTIME = SECOND() 907 Call CC_RDRSP(LISTA,IDLSTA,ISYCTR,IOPT,MODEL, 908 & WORK(KDUM),WORK(KEND2)) 909 TIMIO = TIMIO + SECOND() - DTIME 910 911 Call CC_T2SQ (WORK(KEND2), WORK(KCTR2(IEND)), ISYCTR) 912 913 END IF 914 END DO 915 916*---------------------------------------------------------------------* 917* initialize integral calculation and start loop over calls to ERI: 918*---------------------------------------------------------------------* 919 DTIME = SECOND() 920 921 IF (DIRECT) THEN 922 NTOSYM = 1 923 924 IF (HERDIR) THEN 925 CALL HERDI1(WORK(KEND2),LEND2,IPRERI) 926 ELSE 927 KCCFB1 = KEND2 928 KINDXB = KCCFB1 + MXPRIM*MXCONT 929 KEND2 = KINDXB + (8*MXSHEL*MXCONT + 1)/IRAT 930 LEND2 = LWORK - KEND2 931 CALL ERIDI1(KODCL1,KODCL2,KODBC1,KODBC2,KRDBC1,KRDBC2, 932 * KODPP1,KODPP2,KRDPP1,KRDPP2, 933 * KFREE,LFREE,KEND2,WORK(KCCFB1),WORK(KINDXB), 934 * WORK(KEND2),LEND2,IPRERI) 935 KEND2 = KFREE 936 LEND2 = LFREE 937 END IF 938 ELSE 939 NTOSYM = NSYM 940 END IF 941 942 TIMINT = TIMINT + SECOND() - DTIME 943 944 KENDSV = KEND2 945 LWRKSV = LEND2 946 947 DO ISYMD1 = 1, NTOSYM 948 949 IF (DIRECT) THEN 950 IF (HERDIR) THEN 951 NTOT = MAXSHL 952 ELSE 953 NTOT = MXCALL 954 ENDIF 955 ELSE 956 NTOT = NBAS(ISYMD1) 957 END IF 958 959 DO ILLL = 1, NTOT 960 961 DTIME = SECOND() 962 963 IF (DIRECT) THEN 964 965 KEND2 = KENDSV 966 LEND2 = LWRKSV 967 968 IF (HERDIR) THEN 969 CALL HERDI2(WORK(KEND2),LEND2,INDEXA,ILLL,NUMDIS, 970 & IPRINT) 971 ELSE 972 CALL ERIDI2(ILLL,INDEXA,NUMDIS,0,0D0, 973 * WORK(KODCL1),WORK(KODCL2), 974 * WORK(KODBC1),WORK(KODBC2), 975 * WORK(KRDBC1),WORK(KRDBC2), 976 * WORK(KODPP1),WORK(KODPP2), 977 * WORK(KRDPP1),WORK(KRDPP2), 978 * WORK(KCCFB1),WORK(KINDXB), 979 * WORK(KEND2),LEND2,IPRERI) 980 END IF 981 982 KRECNR = KEND2 983 KEND2 = KRECNR + (NBUFX(0) - 1)/IRAT + 1 984 LEND2 = LWORK - KEND2 985 986 IF (LEND2 .LT. 0) THEN 987 CALL QUIT('Insufficient work space in CC_GMAT.') 988 END IF 989 990 ELSE 991 NUMDIS = 1 992 END IF 993 994 TIMINT = TIMINT + SECOND() - DTIME 995 996*---------------------------------------------------------------------* 997* loop over number of distributions on the disk 998*---------------------------------------------------------------------* 999 DO IDEL2 = 1, NUMDIS 1000 1001 IF(DIRECT) THEN 1002 IDEL = INDEXA(IDEL2) 1003 IF (NOAUXB) THEN 1004 IDUM = 1 1005 CALL IJKAUX(IDEL,IDUM,IDUM,IDUM) 1006 END IF 1007 ISYDEL = ISAO(IDEL) 1008 ELSE 1009 IDEL = IBAS(ISYMD1) + ILLL 1010 ISYDEL = ISYMD1 1011 END IF 1012 1013 ISYDIS = MULD2H(ISYDEL,ISYMOP) 1014C 1015C --------------------------------------------------- 1016C allocate work space for AO integrals, read them in 1017C and transform gamma index with zeroth-order XLAMDH0 1018C --------------------------------------------------- 1019C 1020 KXINT = KEND2 1021 KDSRHF0 = KXINT + NDISAO(ISYDIS) 1022 KDSRHFBC = KDSRHF0 + NDSRHF(ISYDIS) 1023 KEND3 = KDSRHFBC + MDSRHF 1024 LEND3 = LWORK - KEND3 1025 1026 IF (LEND3 .LT. 0) THEN 1027 CALL QUIT('Insufficient work space in CC_GMAT. (3)') 1028 END IF 1029 1030 DTIME = SECOND() 1031 CALL CCRDAO(WORK(KXINT),IDEL,IDEL2,WORK(KEND3),LEND3, 1032 & WORK(KRECNR),DIRECT) 1033 TIMRDAO = TIMRDAO + SECOND() - DTIME 1034 1035 DTIME = SECOND() 1036 CALL CCTRBT(WORK(KXINT),WORK(KDSRHF0),WORK(KLAMDH0),ISYM0, 1037 & WORK(KEND3),LEND3,ISYDIS) 1038 TIMTRBT = TIMTRBT + SECOND() - DTIME 1039 1040C 1041C --------------------------------------- 1042C loop over transformation in this batch: 1043C --------------------------------------- 1044C 1045 DO ITRAN = ISTRT, IEND 1046 1047 TIMTRN(ITRAN) = TIMTRN(ITRAN) - SECOND() 1048 1049 IDLSTA = IGTRAN(1,ITRAN) 1050 IDLSTB = IGTRAN(2,ITRAN) 1051 IDLSTC = IGTRAN(3,ITRAN) 1052 1053 ISYCTR = ILSTSYM(LISTA,IDLSTA) 1054 ISYMTB = ILSTSYM(LISTB,IDLSTB) 1055 ISYMTC = ILSTSYM(LISTC,IDLSTC) 1056 1057 ISYMXB = MULD2H(ISYDIS,ISYMTB) 1058 ISYMXC = MULD2H(ISYDIS,ISYMTC) 1059 1060 ISYTCTB = MULD2H(ISYMTB,ISYMTC) 1061 ISYRES = MULD2H(ISYTCTB,ISYCTR) 1062 1063 KXLMDHB = KLAMDHB(ITRAN) 1064 KXLMDHC = KLAMDHC(ITRAN) 1065 KXLMDPB = KLAMDPB(ITRAN) 1066 KXLMDPC = KLAMDPC(ITRAN) 1067 1068 LSAME = (LISTC.EQ.LISTB .AND. IDLSTC.EQ.IDLSTB) 1069 1070C ----------------------------------------------------- 1071C calculate the first 21CD term (symmetric in B,C): 1072C ----------------------------------------------------- 1073 LLSAME = LSAME 1074 FACT = ONE 1075 DTIME = SECOND() 1076 Call CCG_21CD( WORK(KDELTA1(ITRAN)),ISYRES, 1077 & WORK(KCTR2(ITRAN)), ISYCTR, 1078 & WORK(KLAMDH0), WORK(KLAMDP0), 1079 & WORK(KXLMDHB), WORK(KXLMDPB), ISYMTB, 1080 & WORK(KXLMDHC), WORK(KXLMDPC), ISYMTC, 1081 & WORK(KDSRHF0), ISYDIS, IDEL, ISYDEL, 1082 & LLSAME, FACT, WORK(KEND3), LEND3 ) 1083 TIM21CD = TIM21CD + SECOND() - DTIME 1084 1085 1086C ------------------------------------------------------ 1087C transform gamma to occupied, using the XLAMDHB matrix 1088C and calculate the second 21CD term (symmetric in C,0): 1089C ------------------------------------------------------ 1090 DTIME = SECOND() 1091 CALL CCTRBT(WORK(KXINT),WORK(KDSRHFBC),WORK(KXLMDHB), 1092 & ISYMTB,WORK(KEND3),LEND3,ISYDIS) 1093 TIMTRBT = TIMTRBT + SECOND() - DTIME 1094 1095 FACT = ONE 1096 IF (LSAME) FACT = TWO 1097 LLSAME = .FALSE. 1098 DTIME = SECOND() 1099 Call CCG_21CD( WORK(KDELTA1(ITRAN)),ISYRES, 1100 & WORK(KCTR2(ITRAN)), ISYCTR, 1101 & WORK(KLAMDH0), WORK(KLAMDP0), 1102 & WORK(KXLMDHC), WORK(KXLMDPC), ISYMTC, 1103 & WORK(KLAMDH0), WORK(KLAMDP0), ISYM0, 1104 & WORK(KDSRHFBC),ISYMXB, IDEL, ISYDEL, 1105 & LLSAME, FACT, WORK(KEND3), LEND3 ) 1106 TIM21CD = TIM21CD + SECOND() - DTIME 1107 1108 1109C ----------------------------------------------------- 1110C transform gamma to occupied, using the XLAMDHC matrix 1111C and calculate the third 21CD term (symmetric in B,0): 1112C ----------------------------------------------------- 1113 IF (.NOT. LSAME) THEN 1114 DTIME = SECOND() 1115 CALL CCTRBT(WORK(KXINT),WORK(KDSRHFBC),WORK(KXLMDHC), 1116 & ISYMTC,WORK(KEND3),LEND3,ISYDIS) 1117 TIMTRBT = TIMTRBT + SECOND() - DTIME 1118 1119 FACT = ONE 1120 LLSAME = .FALSE. 1121 DTIME = SECOND() 1122 Call CCG_21CD( WORK(KDELTA1(ITRAN)),ISYRES, 1123 & WORK(KCTR2(ITRAN)), ISYCTR, 1124 & WORK(KLAMDH0), WORK(KLAMDP0), 1125 & WORK(KLAMDH0), WORK(KLAMDP0), ISYM0, 1126 & WORK(KXLMDHB), WORK(KXLMDPB), ISYMTB, 1127 & WORK(KDSRHFBC),ISYMXC, IDEL, ISYDEL, 1128 & LLSAME, FACT, WORK(KEND3), LEND3 ) 1129 TIM21CD = TIM21CD + SECOND() - DTIME 1130 END IF 1131 1132 TIMTRN(ITRAN) = TIMTRN(ITRAN) + SECOND() 1133 1134 END DO ! loop over transformations for this batch 1135 1136 END DO ! IDEL2 = 1, NUMDIS 1137 END DO ! ILLL = 1, NTOT 1138 END DO ! ISYMD1 = 1, NTOSYM 1139 END DO ! loop over batches of transformations 1140 1141 IF (LOCDBG) THEN 1142 WRITE (LUPRI,*) MSGDBG, 1143 & 'Delta vectors at the end of the CC2 part:' 1144 DO ITRAN = 1, NGTRAN 1145 IDLSTA = IGTRAN(1,ITRAN) 1146 IDLSTB = IGTRAN(2,ITRAN) 1147 IDLSTC = IGTRAN(3,ITRAN) 1148 ISYTCTB = MULD2H(ILSTSYM(LISTB,IDLSTB),ILSTSYM(LISTC,IDLSTC)) 1149 ISYRES = MULD2H(ISYTCTB,ILSTSYM(LISTA,IDLSTA)) 1150 WRITE (LUPRI,*) MSGDBG, 'ITRAN = ',ITRAN 1151 WRITE (LUPRI,*) MSGDBG, 'IDLSTA,IDLSTB,IDLSTC=',IDLSTA, 1152 & IDLSTB,IDLSTC 1153 Call CC_PRP(WORK(KDELTA1(ITRAN)),WORK(KDUM),ISYRES,1,0) 1154 WRITE (LUPRI,*) MSGDBG, 'Norm of Delta1 after CC2 part:', 1155 & DSQRT(DDOT(NT1AM(ISYRES),WORK(KDELTA1(ITRAN)),1, 1156 & WORK(KDELTA1(ITRAN)),1)) 1157 END DO 1158 END IF 1159*---------------------------------------------------------------------* 1160* end of CC2 part; recover memory form CTR2 1161* in memory now for each transformation: 1162* 1163* KT1AMP0: singles cluster amplitudes 1164* KLAMDH0: zeroth-order Lambda hole matrix 1165* KLAMDP0: zeroth-order LAmbda particle matrix 1166* KDELTA1: singles excitation part of the result vector 1167* KT1AMPB: single response amplitudes B 1168* KLAMDHB: response Lambda hole matrix for B 1169* KLAMDPB: response LAmbda particle matrix for B 1170* KT1AMPC: single response amplitudes C 1171* KLAMDHC: response Lambda hole matrix for C 1172* KLAMDPC: response LAmbda particle matrix for C 1173* KCTR1 : single response multipliers A 1174* KFOCKB : Fock matrix tF^B(i a) 1175* KFOCKC : Fock matrix tF^C(i a) 1176* KFBCOO : Fock matrix ttF^BC(i j) 1177* KFBCVV : Fock matrix ttF^BC(a b) 1178* 1179* KEND1 : free work space 1180* 1181*---------------------------------------------------------------------* 1182 1183C --------------------------- 1184C here we continue for CCS... 1185C --------------------------- 1186 777 CONTINUE 1187 1188 KEND0 = KEND1 1189 1190*=====================================================================* 1191* the remaining parts are more memory intensive, but to not require 1192* the AO integrals --> we complete the transformations one by one 1193* and save the result one file, or calculate the required dot products 1194*=====================================================================* 1195 1196 DO ITRAN = 1, NGTRAN 1197 1198 TIMTRN(ITRAN) = TIMTRN(ITRAN) - SECOND() 1199 1200 IDLSTA = IGTRAN(1,ITRAN) 1201 IDLSTB = IGTRAN(2,ITRAN) 1202 IDLSTC = IGTRAN(3,ITRAN) 1203 1204 ISYCTR = ILSTSYM(LISTA,IDLSTA) 1205 ISYMTB = ILSTSYM(LISTB,IDLSTB) 1206 ISYMTC = ILSTSYM(LISTC,IDLSTC) 1207 1208 ISYCTB = MULD2H(ISYCTR,ISYMTB) 1209 ISYCTC = MULD2H(ISYCTR,ISYMTC) 1210 1211 ISYTCTB = MULD2H(ISYMTB,ISYMTC) 1212 ISYRES = MULD2H(ISYTCTB,ISYCTR) 1213 1214 ISYMFB = MULD2H(ISYMTB,ISYOVOV) 1215 ISYMFC = MULD2H(ISYMTC,ISYOVOV) 1216 ISYFBC = MULD2H(ISYTCTB,ISYOVOV) 1217C 1218 IF (CCS. OR. CC2) GOTO 888 1219C 1220C ------------------------------------------------------- 1221C allocate work space for one-index transformed integrals 1222C and Lagrangian multipliers: 1223C ------------------------------------------------------- 1224C 1225 KXBIAJK = KEND0 1226 KEBIAJK = KXBIAJK + NTAIKJ(ISYMTB) 1227 KZBIAJK = KEBIAJK + NTAIKJ(ISYMTB) 1228 KCBIAJK = KZBIAJK + NTAIKJ(ISYCTB) 1229 KEND1 = KCBIAJK + NTAIKJ(ISYCTB) 1230 1231 KXCIAJK = KEND1 1232 KECIAJK = KXCIAJK + NTAIKJ(ISYMTC) 1233 KZCIAJK = KECIAJK + NTAIKJ(ISYMTC) 1234 KCCIAJK = KZCIAJK + NTAIKJ(ISYCTC) 1235 KEND1 = KCCIAJK + NTAIKJ(ISYCTC) 1236 1237 LEND1 = LWORK - KEND1 1238 1239*=====================================================================* 1240* CCSD part for DELTA1 1241*=====================================================================* 1242 1243 DTIME = SECOND() 1244C 1245C ------------------------------------------ 1246C DELTBC = <Zeta_2| [[tH^B,T2C], tau_1] |HF> 1247C ------------------------------------------ 1248C 1249* read double excitation part of response vector TAmpC: 1250 KT2AMPC = KEND1 1251 KDELTBC = KT2AMPC + NT2AM(ISYMTC) 1252 KEND2 = KDELTBC + NT1AM(ISYRES) 1253 LEND2 = LWORK - KEND2 1254 1255 IF (LEND2 .LT. 0) THEN 1256 CALL QUIT('Insufficient work space in CC_GMAT.') 1257 END IF 1258 1259 IOPT = 2 1260 Call CC_RDRSP(LISTC,IDLSTC,ISYMTC,IOPT,MODEL, 1261 & WORK(KDUM),WORK(KT2AMPC)) 1262 1263 Call CCLR_DIASCL(WORK(KT2AMPC),TWO,ISYMTC) 1264 1265 Call DZERO(WORK(KDELTBC),NT1AM(ISYRES)) 1266 1267 Call CCG_21CCSD(WORK(KDELTBC),ISYRES,LISTA,IDLSTA, 1268 & WORK(KFOCKB(ITRAN)), ISYMFB, 1269 & WORK(KT1AMPB(ITRAN)),ISYMTB, 1270 & WORK(KT2AMPC), ISYMTC, 1271 & LISTC, IDLSTC, 1272 & WORK(KXBIAJK), WORK(KEBIAJK), 1273 & WORK(KZBIAJK), WORK(KCBIAJK), 1274 & WORK(KEND2), LEND2 ) 1275 1276 LSAME = (LISTC.EQ.LISTB .AND. IDLSTC.EQ.IDLSTB) 1277 1278 IF (.NOT. LSAME) THEN 1279 CALL DAXPY(NT1AM(ISYRES),ONE,WORK(KDELTBC),1, 1280 & WORK(KDELTA1(ITRAN)),1) 1281 ELSE 1282 CALL DAXPY(NT1AM(ISYRES),TWO,WORK(KDELTBC),1, 1283 & WORK(KDELTA1(ITRAN)),1) 1284 1285 CALL DCOPY(NTAIKJ(ISYMTB),WORK(KXBIAJK),1,WORK(KXCIAJK),1) 1286 CALL DCOPY(NTAIKJ(ISYMTB),WORK(KEBIAJK),1,WORK(KECIAJK),1) 1287 CALL DCOPY(NTAIKJ(ISYCTB),WORK(KZBIAJK),1,WORK(KZCIAJK),1) 1288 CALL DCOPY(NTAIKJ(ISYCTB),WORK(KCBIAJK),1,WORK(KCCIAJK),1) 1289 END IF 1290 1291 1292 IF (LOCDBG) THEN 1293 WRITE (LUPRI,*) MSGDBG, '1. CCG_21CCSD contrib. to Delta1:' 1294 WRITE (LUPRI,*) MSGDBG, 'Norm of DELTBC:', 1295 & DSQRT(DDOT(NT1AM(ISYRES),WORK(KDELTBC),1,WORK(KDELTBC),1)) 1296 XNORM = ZERO 1297 IF (ISYMTC.EQ.ISYRES) XNORM = 1298 & DDOT(NT1AM(ISYRES),WORK(KDELTBC),1,WORK(KT1AMPC(ITRAN)),1) 1299 WRITE(LUPRI,*) MSGDBG, 'DELTBC dotted with C amplitudes:', 1300 & XNORM 1301 XNORM = ZERO 1302 IF (ISYMTB.EQ.ISYRES) XNORM = 1303 & DDOT(NT1AM(ISYRES),WORK(KDELTBC),1,WORK(KT1AMPB(ITRAN)),1) 1304 WRITE (LUPRI,*) MSGDBG, 'DELTBC dotted with B amplitudes:', 1305 & XNORM 1306 WRITE (LUPRI,*) MSGDBG, 'DELTBC:' 1307 CALL CC_PRP(WORK(KDELTBC),WORK(KDUM),ISYRES,1,0) 1308 CALL FLSHFO(LUPRI) 1309 END IF 1310C 1311C ------------------------------------------ 1312C DELTCB: <Zeta_2| [[tH^C,T2B], tau_1] |HF> 1313C ------------------------------------------ 1314C 1315 IF (.NOT. LSAME) THEN 1316 1317* read double excitation part of response vector TAmpB: 1318 KT2AMPB = KEND1 1319 KDELTCB = KT2AMPB + NT2AM(ISYMTB) 1320 KEND2 = KDELTCB + NT1AM(ISYRES) 1321 LEND2 = LWORK - KEND2 1322 1323 IF (LEND2 .LT. 0) THEN 1324 CALL QUIT('Insufficient work space in CC_GMAT.') 1325 END IF 1326 1327 IOPT = 2 1328 Call CC_RDRSP(LISTB,IDLSTB,ISYMTB,IOPT,MODEL, 1329 & WORK(KDUM),WORK(KT2AMPB)) 1330 1331 Call CCLR_DIASCL(WORK(KT2AMPB),TWO,ISYMTB) 1332 1333 Call DZERO(WORK(KDELTCB),NT1AM(ISYRES)) 1334 1335 Call CCG_21CCSD (WORK(KDELTCB),ISYRES,LISTA,IDLSTA, 1336 & WORK(KFOCKC(ITRAN)), ISYMFC, 1337 & WORK(KT1AMPC(ITRAN)), ISYMTC, 1338 & WORK(KT2AMPB), ISYMTB, 1339 & LISTB, IDLSTB, 1340 & WORK(KXCIAJK), WORK(KECIAJK), 1341 & WORK(KZCIAJK), WORK(KCCIAJK), 1342 & WORK(KEND2), LEND2 ) 1343 1344 Call DAXPY(NT1AM(ISYRES),ONE,WORK(KDELTCB),1, 1345 & WORK(KDELTA1(ITRAN)),1) 1346 1347 END IF 1348 1349 IF (LOCDBG .AND. .NOT. LSAME) THEN 1350 WRITE (LUPRI,*) MSGDBG, '2. CCG_21CCSD contrib. to Delta1:' 1351 WRITE (LUPRI,*) MSGDBG, 'Norm of DELTCB:', 1352 & DSQRT(DDOT(NT1AM(ISYRES),WORK(KDELTCB),1,WORK(KDELTCB),1)) 1353 WRITE (LUPRI,*) MSGDBG, 'DELTCB dotted with C amplitudes:', 1354 & DDOT(NT1AM(ISYRES),WORK(KDELTCB),1,WORK(KT1AMPC(ITRAN)),1) 1355 WRITE (LUPRI,*) MSGDBG, 'DELTCB dotted with B amplitudes:', 1356 & DDOT(NT1AM(ISYRES),WORK(KDELTCB),1,WORK(KT1AMPB(ITRAN)),1) 1357 WRITE (LUPRI,*) MSGDBG, 'DELTCB:' 1358 Call CC_PRP(WORK(KDELTCB),WORK(KDUM),ISYRES,1,0) 1359 Call FLSHFO(LUPRI) 1360 END IF 1361 1362 TIM21SD = TIM21SD + SECOND() - DTIME 1363 1364 IF (LDOUBL) THEN 1365 1366*=====================================================================* 1367* CCSD part for DELTA2 1368*=====================================================================* 1369 1370 DTIME = SECOND() 1371C 1372C -------------------------------------------------------- 1373C allocate work space for double excitation result vector: 1374C -------------------------------------------------------- 1375C 1376 KDELTA2 = KEND1 1377 KEND1 = KDELTA2 + NT2AM(ISYRES) 1378 LEND1 = LWORK - KEND1 1379 1380 IF (LEND1 .LT. 0) THEN 1381 CALL QUIT('Insufficient work space in CC_GMAT. '// 1382 & '(doubles part)') 1383 END IF 1384 1385* initialize DELTA2 array: 1386 Call DZERO(WORK(KDELTA2), NT2AM(ISYRES)) 1387C 1388C ------------------------------------- 1389C B term: requires packed multipliers 1390C and squared (ia|jb) integrals 1391C ------------------------------------- 1392C 1393 ISYTCTB = MULD2H(ISYMTC,ISYMTB) 1394 ISYC4O = MULD2H(ISYTCTB,ISYCTR) 1395 1396* memory allocation: 1397* 1) squared integrals & scratch space for packed integrals 1398 KXIAJB = KEND1 1399 KSCR = KXIAJB + NT2SQ(ISYOVOV) 1400 KENDB2 = KSCR + NT2AM(ISYOVOV) 1401 LENDB2 = LWORK - KENDB2 1402 1403* 2) packed multipliers & double one-index transf. multipl.: 1404 KZETA2 = KSCR 1405 KC4O = KZETA2 + NT2AM(ISYCTR) 1406 KENDB3 = KC4O + NGAMMA(ISYC4O) 1407 LENDB3 = LWORK - KENDB3 1408 1409 IF ( LENDB2 .LT. 0 .OR. LENDB3 .LT. 0 ) THEN 1410 CALL QUIT('Insufficient work space in CC_GMAT. (B2/B3)') 1411 END IF 1412 1413* read packed (ia|jb) integrals and square them up 1414 Call CCG_RDIAJB(WORK(KSCR),NT2AM(ISYOVOV)) 1415 1416 Call CC_T2SQ (WORK(KSCR), WORK(KXIAJB), ISYOVOV) 1417 1418* read packed multipliers: 1419 KDUM = - 99 999 999 ! dummy address 1420 IOPT = 2 1421 Call CC_RDRSP(LISTA,IDLSTA,ISYCTR,IOPT,MODEL, 1422 & WORK(KDUM),WORK(KZETA2)) 1423 1424* calculate double one-index transformed lagr. multipliers: 1425 IOPT = 1 1426 Call CCG_4O(WORK(KC4O),ISYC4O, WORK(KZETA2),ISYCTR, 1427 & WORK(KT1AMPB(ITRAN)),ISYMTB, 1428 & WORK(KT1AMPC(ITRAN)),ISYMTC, 1429 & WORK(KENDB3),LENDB3,IOPT) 1430 1431* contract them with the (ia|jb) integrals: 1432 IOPT = 2 1433 Call CCRHS_A(WORK(KDELTA2),WORK(KXIAJB),WORK(KC4O), 1434 & WORK(KENDB3), LENDB3, ISYC4O, ISYOVOV, IOPT) 1435 1436C 1437C ----------------------------------------- 1438C A term: requires packed (ia|jb) integrals 1439C and squared multipliers 1440C ----------------------------------------- 1441C 1442* set symmetry for double transformed (ik|jl) integrals 1443 ISYTCTB = MULD2H(ISYMTC,ISYMTB) 1444 ISYX4O = MULD2H(ISYTCTB,ISYOVOV) 1445 1446* memory allocation: 1447* 1) squared multipliers & scratch space for packed multipliers 1448 KZETA2 = KEND1 1449 KSCR = KZETA2 + NT2SQ(ISYCTR) 1450 KENDA2 = KSCR + NT2AM(ISYCTR) 1451 LENDA2 = LWORK - KENDA2 1452 1453* 2) packed (ia|jb) integrals & (ik|jl) integrals 1454 KXIAJB = KSCR 1455 KX4O = KXIAJB + NT2AM(ISYOVOV) 1456 KENDA3 = KX4O + NGAMMA(ISYX4O) 1457 LENDA3 = LWORK - KENDA3 1458 1459 IF ( LENDA2 .LT. 0 .OR. LENDA3 .LT. 0 ) THEN 1460 CALL QUIT('Insufficient work space in CC_GMAT. (A2/A3)') 1461 END IF 1462 1463* read packed lagrangian multipliers and square them up 1464 IOPT = 2 1465 Call CC_RDRSP(LISTA,IDLSTA,ISYCTR,IOPT,MODEL, 1466 & WORK(KDUM),WORK(KSCR)) 1467 1468 Call CC_T2SQ (WORK(KSCR), WORK(KZETA2), ISYCTR) 1469 1470* read (ia|jb) integrals: 1471 Call CCG_RDIAJB(WORK(KXIAJB),NT2AM(ISYOVOV)) 1472 1473 1474* calculate double one-index transformed (ik|jl) integrals: 1475 IOPT = 1 1476 Call CCG_4O(WORK(KX4O),ISYX4O, WORK(KXIAJB),ISYOVOV, 1477 & WORK(KT1AMPB(ITRAN)),ISYMTB, 1478 & WORK(KT1AMPC(ITRAN)),ISYMTC, 1479 & WORK(KENDA3),LENDA3,IOPT) 1480 1481 1482* contract them with the multipliers: 1483 IOPT = 2 1484 Call CCRHS_A(WORK(KDELTA2),WORK(KZETA2),WORK(KX4O), 1485 & WORK(KENDA3),LENDA3,ISYX4O,ISYCTR,IOPT) 1486 1487 IF (LOCDBG) THEN 1488 WRITE (LUPRI,*) MSGDBG, 'finished A term.' 1489 Call CC_PRP(WORK(KDELTA1(ITRAN)),WORK(KDELTA2),ISYRES,1,1) 1490 Call FLSHFO(LUPRI) 1491 END IF 1492C 1493C ------------------------------------- 1494C E term: requires squared multipliers 1495C (still in memory form A term) 1496C ------------------------------------- 1497C 1498 KZETA2 = KEND1 1499 KEMAT1 = KZETA2 + NT2SQ(ISYCTR) 1500 KEMAT2 = KEMAT1 + NMATAB(ISYFBC) 1501 KENDE3 = KEMAT2 + NMATIJ(ISYFBC) 1502 LENDE3 = LWORK - KENDE3 1503 1504 IF ( LENDE3 .LT. 0 ) THEN 1505 CALL QUIT('Insufficient work space in CC_GMAT. (E3)') 1506 END IF 1507 1508* transpose ttF^BC(a b) --> EMAT1(b a) 1509 DO ISYMA = 1, NSYM 1510 ISYMB = MULD2H(ISYMA,ISYFBC) 1511 DO A = 1, NVIR(ISYMA) 1512 DO B = 1, NVIR(ISYMB) 1513 NAB = IMATAB(ISYMA,ISYMB) + NVIR(ISYMA)*(B-1) + A 1514 NBA = IMATAB(ISYMB,ISYMA) + NVIR(ISYMB)*(A-1) + B 1515 1516 WORK(KEMAT1 - 1 + NBA) = WORK(KFBCVV(ITRAN) - 1 + NAB) 1517 END DO 1518 END DO 1519 END DO 1520 1521 1522* transpose ttF^BC(i j) --> EMAT2(j i) 1523 DO ISYMI = 1, NSYM 1524 ISYMJ = MULD2H(ISYMI,ISYFBC) 1525 DO I = 1, NRHF(ISYMI) 1526 DO J = 1, NRHF(ISYMJ) 1527 NIJ = IMATIJ(ISYMI,ISYMJ) + NRHF(ISYMI)*(J-1) + I 1528 NJI = IMATIJ(ISYMJ,ISYMI) + NRHF(ISYMJ)*(I-1) + J 1529 1530 WORK(KEMAT2 - 1 + NJI) = WORK(KFBCOO(ITRAN) - 1 + NIJ) 1531 END DO 1532 END DO 1533 END DO 1534 1535* combine EMAT1/EMAT2 with lagrangian multipliers: 1536 Call CCRHS_E(WORK(KDELTA2),WORK(KZETA2),WORK(KEMAT1), 1537 & WORK(KEMAT2), WORK(KENDE3),LENDE3,ISYCTR,ISYFBC) 1538 1539 IF (LOCDBG) THEN 1540 WRITE (LUPRI,*) MSGDBG, 'finished E term.' 1541 Call CC_PRP(WORK(KDELTA1(ITRAN)),WORK(KDELTA2),ISYRES,1,1) 1542 Call FLSHFO(LUPRI) 1543 END IF 1544 1545*---------------------------------------------------------------------* 1546* D term & C term: 1547*---------------------------------------------------------------------* 1548 1549* recover work space from integrals and Lagrangian multipliers, 1550* and allocate work space for C & D term as square matrix: 1551 KCDTERM = KZETA2 1552 KEND2 = KCDTERM + NT2SQ(ISYRES) 1553 LEND2 = LWORK - KEND2 1554 1555 IF ( LEND2 .LT. 0 ) THEN 1556 CALL QUIT('Insufficient work space in CC_GMAT. (C & D)') 1557 END IF 1558 1559 Call CCG_22CD(WORK(KDELTA2),WORK(KCDTERM),ISYRES, 1560 & WORK(KXBIAJK),WORK(KEBIAJK),ISYMTB, 1561 & WORK(KZBIAJK),WORK(KCBIAJK),ISYCTB, 1562 & WORK(KXCIAJK),WORK(KECIAJK),ISYMTC, 1563 & WORK(KZCIAJK),WORK(KCCIAJK),ISYCTC, 1564 & LSAME, WORK(KEND2), LEND2) 1565 1566 1567 TIMDBL = TIMDBL + SECOND() - DTIME 1568 1569 END IF ! (LDOUBL) 1570 1571*=====================================================================* 1572* add CC3 corrections: 1573*=====================================================================* 1574 DTIME = SECOND() 1575 1576 IF (CCSDT) THEN 1577 IF (NODDY_GMAT) THEN 1578 CALL CCSDT_GMAT_NODDY(LISTA,IDLSTA,LISTB,IDLSTB, 1579 & LISTC,IDLSTC, 1580 & WORK(KDELTA1(ITRAN)),WORK(KDELTA2), 1581 & WORK(KEND1),LEND1) 1582 ELSE 1583 1584 LUCKJD = -1 1585 LUTOC = -1 1586 LU3VI = -1 1587 LUDKBC3 = -1 1588 LU3FOPX = -1 1589 LU3FOP2X = -1 1590C 1591 CALL WOPEN2(LUCKJD,FNCKJD,64,0) 1592 CALL WOPEN2(LUTOC,FNTOC,64,0) 1593 CALL WOPEN2(LU3VI,FN3VI,64,0) 1594 CALL WOPEN2(LUDKBC3,FNDKBC3,64,0) 1595 CALL WOPEN2(LU3FOPX,FN3FOPX,64,0) 1596 CALL WOPEN2(LU3FOP2X,FN3FOP2X,64,0) 1597 1598 CALL CC3_GMATNEW(LISTA,IDLSTA,LISTB,IDLSTB, 1599 * LISTC,IDLSTC, 1600 * WORK(KDELTA1(ITRAN)),WORK(KDELTA2), 1601 * ISYRES, 1602 * WORK(KEND1),LEND1, 1603 * LUTOC,FNTOC, 1604 * LU3VI,FN3VI,LUDKBC3,FNDKBC3, 1605 * LU3FOPX,FN3FOPX,LU3FOP2X,FN3FOP2X) 1606 1607 CALL WCLOSE2(LUCKJD,FNCKJD,'KEEP') 1608 CALL WCLOSE2(LUTOC,FNTOC,'KEEP') 1609 CALL WCLOSE2(LU3VI,FN3VI,'KEEP') 1610 CALL WCLOSE2(LUDKBC3,FNDKBC3,'KEEP') 1611 CALL WCLOSE2(LU3FOPX,FN3FOPX,'KEEP') 1612 CALL WCLOSE2(LU3FOP2X,FN3FOP2X,'KEEP') 1613 1614 1615 END IF 1616 END IF 1617 1618 TIMTRIP = TIMTRIP + SECOND() - DTIME 1619 1620C 1621C CCSLV part, JK+OC 1622C The result vectors are ind WORK(KDELTA1(ITRAN)) and WORK(KDELTA2) 1623C Futhermore, we need to know the list numbers, their ID's and symmetries 1624C For the zeta (A) vector these are: LISTA,IDLSTA,ISYCTR 1625C For the B trial vector these are : LISTB,IDLSTB,ISYMTB 1626C For the C trial vector these are : LISTC,IDLSTC,ISYMTC 1627C Also, we have as arguments the symmetry of B x C : ISYTCTB 1628C and the symmetry of A x B X C : ISYRES 1629C work arrays used are WORK(KEND2) and LEND2 1630C 1631 IF (CCSLV) THEN 1632 IF (.NOT. CCMM) CALL CCSL_GTR(WORK(KDELTA1(ITRAN)), 1633 * WORK(KDELTA2),ISYRES, 1634 * LISTA,IDLSTA,ISYCTR, 1635 * LISTB,IDLSTB,ISYMTB, 1636 * LISTC,IDLSTC,ISYMTC, 1637 * MODEL,WORK(KEND2),LEND2) 1638C 1639 IF (CCMM) THEN 1640 IF (.NOT. NYQMMM) THEN 1641 CALL CCMM_GTR(WORK(KDELTA1(ITRAN)), 1642 * WORK(KDELTA2),ISYRES, 1643 * LISTA,IDLSTA,ISYCTR, 1644 * LISTB,IDLSTB,ISYMTB, 1645 * LISTC,IDLSTC,ISYMTC, 1646 * MODEL,WORK(KEND2),LEND2) 1647C 1648 ELSE IF (NYQMMM) THEN 1649 RSPTYP='G' 1650 CALL CCMM_QRTRANSFORMER(WORK(KDELTA1(ITRAN)), 1651 * WORK(KDELTA2),ISYRES, 1652 * LISTB,IDLSTB,ISYMTB, 1653 * LISTC,IDLSTC,ISYMTC, 1654 * MODEL,RSPTYP,WORK(KEND2),LEND2) 1655 END IF 1656 END IF 1657 END IF 1658 IF (USE_PELIB()) THEN 1659 RSPTYP='G' 1660 CALL PELIB_IFC_QRTRANSFORMER(WORK(KDELTA1(ITRAN)), 1661 & WORK(KDELTA2),ISYRES,LISTB,IDLSTB,ISYMTB, 1662 & LISTC,IDLSTC,ISYMTC,MODEL,RSPTYP,WORK(KEND2),LEND2) 1663 END IF 1664 1665*=====================================================================* 1666* CCSD part for CC-R12 1667*=====================================================================* 1668 1669 IF (CCSDR12) THEN 1670C 1671C ------------------------------------------------------------ 1672C allocate work space for r12 double excitation result vector 1673C right after conv. doubles and reset KEND1: 1674C ------------------------------------------------------------ 1675C 1676 KDELTAR12 = KEND1 1677 KEND1 = KDELTAR12 + NTR12AM(ISYRES) 1678 LEND1 = LWORK - KEND1 1679 1680 IF (LEND1 .LT. 0) THEN 1681 CALL QUIT('Insufficient work space in CC_GMAT. '// 1682 & '(r12 doubles part)') 1683 END IF 1684 1685* initialize DELTAR12 array: 1686 Call DZERO(WORK(KDELTAR12), NTR12AM(ISYRES)) 1687C 1688C add B' term 1689C 1690 !read V_(alpha beta)^(kl) from disk 1691 KVABKL = KEND1 1692 KEND2 = KVABKL + NVABKL(1) 1693 LEND2 = LWORK - KEND2 1694 IF (LEND2 .LT. 0) THEN 1695 CALL QUIT('Insuff. work space for V^(albe)_kl in CC_GMAT') 1696 END IF 1697 lunit = -1 1698 call gpopen(lunit,fvabkl,'unknown',' ','unformatted', 1699 & idum,.false.) 1700 read (lunit)(work(kvabkl+i-1),i=1,nvabkl(1)) 1701 call gpclose(lunit,'KEEP') 1702C 1703 IOPT = 2 1704 CALL CC_R12MKVIRT(WORK(KVABKL),WORK(KLAMDPB(ITRAN)),ISYMTB, 1705 & WORK(KLAMDPC(ITRAN)),ISYMTC, 1706 & 'R12VCBDBKL',IOPT,WORK(KEND2),LEND2) 1707C 1708 !read doubles Lagrangian multipliers 1709 KZETA2 = KEND1 1710 KR12SCR= KZETA2 + NT2AM(ISYCTR) 1711 KEND2 = KR12SCR+ NTR12SQ(ISYRES) 1712 LEND2 = LWORK - KEND2 1713 IF (LEND2 .LT. 0) THEN 1714 CALL QUIT('Insuff. work space in CC_GMAT') 1715 END IF 1716C 1717 CALL DZERO(WORK(KR12SCR),NTR12SQ(ISYRES)) 1718C 1719 IOPT = 2 1720 CALL CC_RDRSP(LISTA,IDLSTA,ISYCTR,IOPT,MODEL, 1721 & WORK(KDUM),WORK(KZETA2)) 1722C 1723 !calculate contribution 1724 CALL CCRHS_BPP(WORK(KR12SCR),WORK(KZETA2),ISYCTR,.TRUE., 1725 & 'R12VCBDBKL',ISYTCTB,WORK(KEND2),LEND2) 1726C 1727C resort result into a symmetry packed triangular matrix 1728C 1729 IOPT = 1 1730 CALL CCR12PCK2(WORK(KDELTAR12),ISYRES,.FALSE., 1731 & WORK(KR12SCR),'T',IOPT) 1732 CALL CCLR_DIASCLR12(WORK(KDELTAR12),0.5D0*KETSCL,ISYRES) 1733C 1734 END IF !(CCSDR12) 1735 1736*=====================================================================* 1737* end CCSD part for DELTA2 / DELTAR12; save result vector on file or 1738* calculate the required dot products: 1739*=====================================================================* 1740C 1741C --------------------------------- 1742C here we continue for CCS and CC2: 1743C --------------------------------- 1744 888 CONTINUE 1745 IF (CC2) THEN 1746 KDELTA2 = KEND0 1747 KEND1 = KDELTA2 + NT2AM(ISYRES) 1748 LEND1 = LWORK - KEND1 1749 1750 IF (LEND1 .LT. 0) THEN 1751 CALL QUIT('Insufficient work space in CC_GMAT. '// 1752 & '(save section)') 1753 END IF 1754 1755 CALL DZERO(WORK(KDELTA2),NT2AM(ISYRES)) 1756C 1757 IF (CCSLV) THEN 1758 IF (.NOT. CCMM) CALL CCSL_GTR(WORK(KDELTA1(ITRAN)), 1759 * WORK(KDELTA2),ISYRES, 1760 * LISTA,IDLSTA,ISYCTR, 1761 * LISTB,IDLSTB,ISYMTB, 1762 * LISTC,IDLSTC,ISYMTC, 1763 * MODEL,WORK(KEND1),LEND1) 1764C 1765 IF (CCMM) THEN 1766 IF (.NOT.NYQMMM) THEN 1767 CALL CCMM_GTR(WORK(KDELTA1(ITRAN)), 1768 & WORK(KDELTA2),ISYRES, 1769 & LISTA,IDLSTA,ISYCTR, 1770 & LISTB,IDLSTB,ISYMTB, 1771 & LISTC,IDLSTC,ISYMTC, 1772 & MODEL,WORK(KEND1),LEND1) 1773 ELSE IF (NYQMMM) THEN 1774 RSPTYP = 'G' 1775 CALL CCMM_QRTRANSFORMER(WORK(KDELTA1(ITRAN)), 1776 & WORK(KDELTA2),ISYRES, 1777 & LISTB,IDLSTB,ISYMTB, 1778 & LISTC,IDLSTC,ISYMTC, 1779 & MODEL,RSPTYP,WORK(KEND1),LEND1) 1780 END IF 1781 END IF 1782 END IF 1783 IF (USE_PELIB()) THEN 1784 RSPTYP = 'G' 1785 CALL PELIB_IFC_QRTRANSFORMER(WORK(KDELTA1(ITRAN)), 1786 & WORK(KDELTA2),ISYRES,LISTB,IDLSTB,ISYMTB, 1787 & LISTC,IDLSTC,ISYMTC,MODEL,RSPTYP,WORK(KEND1),LEND1) 1788 END IF 1789 END IF 1790C 1791 IF (LOCDBG) THEN 1792 WRITE (LUPRI,*) 1793 & 'Final result for G matrix transformation ',ITRAN 1794 Call CC_PRP(WORK(KDELTA1(ITRAN)),WORK(KDELTA2),ISYRES,1,1) 1795 IF (CCSDR12) CALL CC_PRPR12(WORK(KDELTAR12),ISYRES,1,.TRUE.) 1796 END IF 1797 1798 IFILE = IGTRAN(4,ITRAN) 1799 IDLSTD = IGTRAN(4,ITRAN) 1800 LEN1 = NT1AM(ISYRES) 1801 LEN2 = NT2AM(ISYRES) 1802 LENR12 = NTR12AM(ISYRES) 1803 1804 IF (CCS) THEN 1805 IOPTW = 1 1806 MODELW = 'CCS ' 1807 ELSE IF (CC2) THEN 1808 IOPTW = 3 1809 MODELW = 'CC2 ' 1810 ELSE IF (CCSD) THEN 1811 IOPTW = 3 1812 MODELW = 'CCSD ' 1813 ELSE IF (CC3) THEN 1814 IOPTW = 3 1815 IOPTWE = 24 1816 MODELW = 'CC3 ' 1817 ELSE 1818 CALL QUIT('Unknown coupled cluster model in CC_GMAT.') 1819 END IF 1820 IF (CCR12) THEN 1821 APROXR12 = ' ' 1822 CALL CCSD_MODEL(MODELW,LENMOD,10,MODELW,10,APROXR12) 1823 IOPTWR12 = 32 1824 END IF 1825 1826 IF ( IOPTRES.EQ.0 .OR. IOPTRES.EQ.1 ) THEN 1827C 1828C ---------------------------------- 1829C save result on direct access file: 1830C ---------------------------------- 1831C 1832 IGTRAN(4,ITRAN) = IADRDEL 1833 1834 DTIME = SECOND() 1835 1836 CALL PUTWA2(LUGMAT,FILGMA,WORK(KDELTA1(ITRAN)),IADRDEL,LEN1) 1837 IADRDEL = IADRDEL + LEN1 1838 1839 IF (.NOT.CCS) THEN 1840 CALL PUTWA2(LUGMAT,FILGMA,WORK(KDELTA2),IADRDEL,LEN2) 1841 IADRDEL = IADRDEL + LEN2 1842 END IF 1843 1844 IF (CCSDR12) THEN 1845 CALL PUTWA2(LUGMAT,FILGMA,WORK(KDELTAR12),IADRDEL,LENR12) 1846 IADRDEL = IADRDEL + LENR12 1847 END IF 1848 1849 TIMIO = TIMIO + SECOND() - DTIME 1850 1851 ELSE IF (IOPTRES.EQ.3) THEN 1852C 1853C ----------------------------- 1854C write to file using CC_WRRSP: 1855C ----------------------------- 1856C 1857 DTIME = SECOND() 1858 CALL CC_WRRSP(LISTD,IFILE,ISYRES,IOPTW,MODELW,WORK(KDUM), 1859 & WORK(KDELTA1(ITRAN)),WORK(KDELTA2),WORK(KEND1),LEND1) 1860 IF (CCSDR12) THEN 1861 CALL CC_WRRSP(LISTD,IFILE,ISYRES,IOPTWR12,MODELW, 1862 & WORK(KDUM),WORK(KDUM),WORK(KDELTAR12), 1863 & WORK(KEND1),LEND1) 1864 END IF 1865 IF (CCSDT) THEN 1866 CALL CC_WRRSP(LISTD,IFILE,ISYRES,IOPTWE,MODELW,WORK(KDUM), 1867 & WORK(KDELTA1(ITRAN)),WORK(KDELTA2),WORK(KEND1),LEND1) 1868 END IF 1869 TIMIO = TIMIO + SECOND() - DTIME 1870 1871 ELSE IF (IOPTRES.EQ.4) THEN 1872C 1873C ----------------------------- 1874C write to file using CC_WARSP: 1875C ----------------------------- 1876C 1877 DTIME = SECOND() 1878 CALL CC_WARSP(LISTD,IFILE,ISYRES,IOPTW,MODELW,WORK(KDUM), 1879 & WORK(KDELTA1(ITRAN)),WORK(KDELTA2),WORK(KEND1),LEND1) 1880 IF (CCSDR12) THEN 1881 CALL CC_WARSP(LISTD,IFILE,ISYRES,IOPTWR12,MODELW, 1882 & WORK(KDUM),WORK(KDUM),WORK(KDELTAR12), 1883 & WORK(KEND1),LEND1) 1884 END IF 1885 IF (CCSDT) THEN 1886 CALL CC_WARSP(LISTD,IFILE,ISYRES,IOPTWE,MODELW,WORK(KDUM), 1887 & WORK(KDELTA1(ITRAN)),WORK(KDELTA2),WORK(KEND1),LEND1) 1888 END IF 1889 TIMIO = TIMIO + SECOND() - DTIME 1890 1891 ELSE IF (IOPTRES.EQ.5) THEN 1892C 1893C -------------------------------- 1894C calculate required dot products: 1895C -------------------------------- 1896C 1897 IF (.NOT.CCS) CALL CCLR_DIASCL(WORK(KDELTA2),TWO,ISYRES) 1898 CALL CCDOTRSP(IGDOTS,GCON,IOPTW,LISTD,ITRAN,NGTRAN,MXVEC, 1899 & WORK(KDELTA1(ITRAN)),WORK(KDELTA2),ISYRES, 1900 & WORK(KEND1),LEND1) 1901 IF (CCSDR12) THEN 1902 CALL CCDOTRSP(IGDOTS,GCON,IOPTWR12,LISTD,ITRAN,NGTRAN, 1903 & MXVEC,WORK(KDUM),WORK(KDELTAR12),ISYRES, 1904 & WORK(KEND1),LEND1) 1905 END IF 1906 1907 ELSE 1908 CALL QUIT('Illegal value for IOPTRES in CC_GMAT.') 1909 END IF 1910 1911 TIMTRN(ITRAN) = TIMTRN(ITRAN) + SECOND() 1912 1913 IF (IPRINT.GT.0) THEN 1914 IF (IOPTRES.EQ.5) THEN 1915 IVEC = 1 1916 DO WHILE (IGDOTS(IVEC,ITRAN).NE.0 .AND. IVEC.LE.MXVEC) 1917 IVEC = IVEC + 1 1918 END DO 1919 WRITE (LUPRI,'(1X,3(A,I5),A,I6,A,F10.2,A)') '| ',IDLSTA, 1920 & ' | ',IDLSTB,' | ',IDLSTC,' | ',IVEC-1, 1921 & ' | ',TIMTRN(ITRAN),' |' 1922 ELSE 1923 WRITE (LUPRI,'(1X,3(A,I5),A,I6,A,F10.2,A)') '| ',IDLSTA, 1924 & ' | ',IDLSTB,' | ',IDLSTC,' | ',IDLSTD, 1925 & ' | ',TIMTRN(ITRAN),' |' 1926 END IF 1927 END IF 1928 1929 END DO ! ITRAN 1930 1931*=====================================================================* 1932* finished with all G matrix transformations for this time: 1933* for IOPTRES = 1 and enough memory read result back into memory, 1934* in any case close G matrix file & return 1935*=====================================================================* 1936 1937* check size of work space: 1938 IF (IOPTRES .EQ. 1) THEN 1939 LENALL = IADRDEL - 1 1940 IF (LENALL .GT. LWORK) IOPTRES = 0 1941 END IF 1942 1943* read result vectors back into memory: 1944 IF (IOPTRES .EQ. 1) THEN 1945 DTIME = SECOND() 1946 CALL GETWA2(LUGMAT,FILGMA,WORK,1,LENALL) 1947 TIMIO = TIMIO + SECOND() - DTIME 1948 END IF 1949 1950* close G matrix file, if open: 1951 IF ( IOPTRES.EQ.0 .OR. IOPTRES.EQ.1 ) THEN 1952 1953 CALL WCLOSE2(LUGMAT,FILGMA,'KEEP') 1954 1955 END IF 1956 1957*---------------------------------------------------------------------* 1958* print timings and return: 1959*---------------------------------------------------------------------* 1960 TIMALL = SECOND() - TIMALL 1961 1962 IF (IPRINT.GE.0) THEN 1963 WRITE (LUPRI,'(1X,A1,50("-"),A1)') '+','+' 1964 WRITE (LUPRI,'(1X,A,I4,A,F10.2,A)') 1965 & '| total time for',NGTRAN,' G transforms.:',TIMALL,' secs.|' 1966 IF (TIMALL.GT.1.0D0) THEN 1967 WRITE (LUPRI,'(1X,A1,50("-"),A1)') '+','+' 1968 CONVRT = 100.0D0/TIMALL 1969 WRITE (LUPRI, 1970 & '(1X,"| % of time used in ",A,":",F10.2," |")') 1971 & 'ERI ', TIMINT*CONVRT, 1972 & 'CCRDAO ', TIMRDAO*CONVRT, 1973 & 'CCTRBT ', TIMTRBT*CONVRT, 1974 & 'CCG_21CD ', TIM21CD*CONVRT, 1975 & 'Fock interm. ', TIMFCK*CONVRT, 1976 & 'CCS part ', TIMCCS*CONVRT, 1977 & '21 CCSD part ', TIM21SD*CONVRT, 1978 & 'doubles part ', TIMDBL*CONVRT, 1979 & 'triples part ', TIMTRIP*CONVRT, 1980 & 'I/O ', TIMIO*CONVRT 1981 END IF 1982 1983 WRITE (LUPRI,'(1X,A1,50("="),A1,//)') '+','+' 1984 END IF 1985 1986 CALL QEXIT('CC_GMAT') 1987 1988 RETURN 1989 END 1990*---------------------------------------------------------------------* 1991c/* Deck CCG_1ITROO */ 1992*=====================================================================* 1993 SUBROUTINE CCG_1ITROO(TFMAT,ISYMTF,FMAT,ISYMFM,T1,ISYMT1) 1994*---------------------------------------------------------------------* 1995* 1996* Purpose: calculate occupied/occupied block of one-index 1997* transformed (Fock) matrix F: 1998* 1999* tF(ij) = sum(b) F(ib) * T1(bj) = [F,T1]_ij 2000* 2001* needed, e.g., for double one-index transformed Fock mat. 2002* 2003* Christof Haettig, August 1996 2004*=====================================================================* 2005#if defined (IMPLICIT_NONE) 2006 IMPLICIT NONE 2007#else 2008# include "implicit.h" 2009#endif 2010#include "priunit.h" 2011#include "ccorb.h" 2012#include "ccsdsym.h" 2013 2014 INTEGER ISYMTF, ISYMT1, ISYMFM 2015 2016#if defined (SYS_CRAY) 2017 REAL TFMAT(*) ! dimension (NMATIJ(ISYMTF)) 2018 REAL FMAT(*) ! dimension (NT1AM(ISYMFM)) 2019 REAL T1(*) ! dimension (NT1AM(ISYMT1)) 2020#else 2021 DOUBLE PRECISION TFMAT(*) ! dimension (NMATIJ(ISYMTF)) 2022 DOUBLE PRECISION FMAT(*) ! dimension (NT1AM(ISYMFM)) 2023 DOUBLE PRECISION T1(*) ! dimension (NT1AM(ISYMT1)) 2024#endif 2025 2026 INTEGER ISYMI, ISYMJ, ISYMB, NIJ, NBI, NBJ 2027 2028 CALL QENTER('CCG_1ITROO') 2029 2030* check symmetries: 2031 IF (ISYMTF .NE. MULD2H(ISYMFM,ISYMT1)) THEN 2032 CALL QUIT('Symmetry mismatch in CCG_1ITROO.') 2033 END IF 2034 2035* initialize resulting matrix: 2036 CALL DZERO(TFMAT, NMATIJ(ISYMTF)) 2037 2038 DO ISYMI = 1, NSYM 2039 ISYMJ = MULD2H(ISYMI,ISYMTF) 2040 ISYMB = MULD2H(ISYMI,ISYMFM) 2041 2042 DO I = 1, NRHF(ISYMI) 2043 DO J = 1, NRHF(ISYMJ) 2044 NIJ = IMATIJ(ISYMI,ISYMJ) + NRHF(ISYMI)*(J - 1) + I 2045 2046 DO B = 1, NVIR(ISYMB) 2047 NBI = IT1AM(ISYMB,ISYMI) + NVIR(ISYMB)*(I - 1) + B 2048 NBJ = IT1AM(ISYMB,ISYMJ) + NVIR(ISYMB)*(J - 1) + B 2049 TFMAT(NIJ) = TFMAT(NIJ) + FMAT(NBI) * T1(NBJ) 2050 END DO 2051 2052 END DO 2053 END DO 2054 2055 END DO 2056 2057 CALL QEXIT('CCG_1ITROO') 2058 2059 RETURN 2060 END 2061*---------------------------------------------------------------------* 2062c/* Deck CCG_1ITRVO */ 2063*=====================================================================* 2064 SUBROUTINE CCG_1ITRVO(tF,ISYMTF,FOO,FVV,ISYMFM,T1,ISYMT1,IOPT) 2065*---------------------------------------------------------------------* 2066* 2067* Purpose: calculate vir/occ or occ/vir block of one-index 2068* transformed (Fock) matrix F: 2069* 2070* iopt=1 : tF(ai) = + sum_b F(ab) * T1(bi) - sum_j T1(aj) * F(ji) 2071* 2072* iopt=2 : tF(ia) = + sum_b T1(ib) * F(ba) - sum_j F(ij) * T1(ja) 2073* 2074* needed, e.g., for Delta1 vector in CCQR_G 2075* 2076* iopt=1 referes to a one-index transformation of F with a vir/occ 2077* matrix (e.g. single excitation amplitudes) 2078* 2079* iopt=2 referes to a one-index transformation of F with a occ/vir 2080* matrix (e.g. single excitation part of zeta vector) 2081* 2082* Christof Haettig, August 1996 2083*=====================================================================* 2084#if defined (IMPLICIT_NONE) 2085 IMPLICIT NONE 2086#else 2087# include "implicit.h" 2088#endif 2089#include "priunit.h" 2090#include "ccorb.h" 2091#include "ccsdsym.h" 2092 2093 INTEGER ISYMTF, ISYMT1, ISYMFM, IOPT 2094 2095#if defined (SYS_CRAY) 2096 REAL TF(*) ! dimension (NT1AM(ISYMTF)) 2097 REAL FOO(*) ! dimension (NMATIJ(ISYMFM)) 2098 REAL FVV(*) ! dimension (NMATAB(ISYMFM)) 2099 REAL T1(*) ! dimension (NT1AM(ISYMT1)) 2100#else 2101 DOUBLE PRECISION TF(*) ! dimension (NT1AM(ISYMTF)) 2102 DOUBLE PRECISION FOO(*) ! dimension (NMATIJ(ISYMFM)) 2103 DOUBLE PRECISION FVV(*) ! dimension (NMATAB(ISYMFM)) 2104 DOUBLE PRECISION T1(*) ! dimension (NT1AM(ISYMT1)) 2105#endif 2106 2107 INTEGER ISYMI,ISYMJ,ISYMA,ISYMB,NAI,NAB,NBA,NBI,NAJ,NJI,NIJ 2108 2109 CALL QENTER('CCG_1ITRVO') 2110 2111* check symmetries: 2112 IF (ISYMTF .NE. MULD2H(ISYMFM,ISYMT1)) THEN 2113 CALL QUIT('Symmetry mismatch in CCG_1ITRVO.') 2114 END IF 2115 2116* initialize resulting matrix: 2117 CALL DZERO(tF, NT1AM(ISYMTF)) 2118 2119 DO ISYMI = 1, NSYM 2120 ISYMA = MULD2H(ISYMI,ISYMTF) 2121 ISYMB = MULD2H(ISYMI,ISYMT1) 2122 ISYMJ = MULD2H(ISYMA,ISYMT1) 2123 2124 IF (IOPT .EQ. 1) THEN 2125 2126 DO I = 1, NRHF(ISYMI) 2127 DO A = 1, NVIR(ISYMA) 2128 NAI = IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I - 1) + A 2129 2130 DO B = 1, NVIR(ISYMB) 2131 NBI = IT1AM(ISYMB,ISYMI) + NVIR(ISYMB)*(I - 1) + B 2132 NAB = IMATAB(ISYMA,ISYMB) + NVIR(ISYMA)*(B - 1) + A 2133 tF(NAI) = tF(NAI) + FVV(NAB) * T1(NBI) 2134 END DO 2135 2136 DO J = 1, NRHF(ISYMJ) 2137 NAJ = IT1AM(ISYMA,ISYMJ) + NVIR(ISYMA)*(J - 1) + A 2138 NJI = IMATIJ(ISYMJ,ISYMI) + NRHF(ISYMJ)*(I - 1) + J 2139 tF(NAI) = tF(NAI) - T1(NAJ) * FOO(NJI) 2140 END DO 2141 2142 END DO 2143 END DO 2144 2145 ELSE IF (IOPT .EQ. 2) THEN 2146 2147 DO I = 1, NRHF(ISYMI) 2148 DO A = 1, NVIR(ISYMA) 2149 NAI = IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I - 1) + A 2150 2151 DO B = 1, NVIR(ISYMB) 2152 NBI = IT1AM(ISYMB,ISYMI) + NVIR(ISYMB)*(I - 1) + B 2153 NBA = IMATAB(ISYMB,ISYMA) + NVIR(ISYMB)*(A - 1) + B 2154 tF(NAI) = tF(NAI) + T1(NBI) * FVV(NBA) 2155 END DO 2156 2157 DO J = 1, NRHF(ISYMJ) 2158 NAJ = IT1AM(ISYMA,ISYMJ) + NVIR(ISYMA)*(J - 1) + A 2159 NIJ = IMATIJ(ISYMI,ISYMJ) + NRHF(ISYMI)*(J - 1) + I 2160 tF(NAI) = tF(NAI) - FOO(NIJ) * T1(NAJ) 2161 END DO 2162 2163 END DO 2164 END DO 2165 2166 ELSE 2167 CALL QUIT('CCG_1ITRVO called with illegal value for IOPT.') 2168 ENDIF 2169 2170 END DO 2171 2172 CALL QEXIT('CCG_1ITRVO') 2173 2174 RETURN 2175 END 2176*---------------------------------------------------------------------* 2177c/* Deck CCG_1ITRVV */ 2178*=====================================================================* 2179 SUBROUTINE CCG_1ITRVV(TFMAT,ISYMTF,FMAT,ISYMFM,T1,ISYMT1) 2180*---------------------------------------------------------------------* 2181* 2182* Purpose: calculate virtual/virtual block of one-index 2183* transformed (Fock) matrix F: 2184* 2185* tF(ab) = - sum(i) T1(ai) * F(ib) = [F,T1]_ab 2186* 2187* needed, e.g., for double one-index transformed Fock mat. 2188* 2189* Christof Haettig, August 1996 2190*=====================================================================* 2191#if defined (IMPLICIT_NONE) 2192 IMPLICIT NONE 2193#else 2194# include "implicit.h" 2195#endif 2196#include "priunit.h" 2197#include "ccorb.h" 2198#include "ccsdsym.h" 2199 2200 INTEGER ISYMTF, ISYMT1, ISYMFM 2201 2202#if defined (SYS_CRAY) 2203 REAL TFMAT(*) 2204 REAL FMAT(*), T1(*) 2205#else 2206 DOUBLE PRECISION TFMAT(*) 2207 DOUBLE PRECISION FMAT(*), T1(*) 2208#endif 2209 2210 INTEGER ISYMA, ISYMI, ISYMB, NAI, NBI, NAB 2211 2212 CALL QENTER('CCG_1ITRVV') 2213 2214* check symmetries: 2215 IF (ISYMTF .NE. MULD2H(ISYMFM,ISYMT1)) THEN 2216 CALL QUIT('Symmetry mismatch in CCG_1ITRVV.') 2217 END IF 2218 2219* initialize resulting matrix: 2220 CALL DZERO(TFMAT, NMATAB(ISYMTF)) 2221 2222 DO ISYMA = 1, NSYM 2223 ISYMB = MULD2H(ISYMA,ISYMTF) 2224 ISYMI = MULD2H(ISYMA,ISYMT1) 2225 2226 DO A = 1, NVIR(ISYMA) 2227 DO B = 1, NVIR(ISYMB) 2228 NAB = IMATAB(ISYMA,ISYMB) + NVIR(ISYMA)*(B - 1) + A 2229 2230 DO I = 1, NRHF(ISYMI) 2231 NAI = IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I - 1) + A 2232 NBI = IT1AM(ISYMB,ISYMI) + NVIR(ISYMB)*(I - 1) + B 2233 TFMAT(NAB) = TFMAT(NAB) - T1(NAI) * FMAT(NBI) 2234 END DO 2235 2236 END DO 2237 END DO 2238 2239 END DO 2240 2241 CALL QEXIT('CCG_1ITRVV') 2242 2243 RETURN 2244 END 2245*---------------------------------------------------------------------* 2246c/* Deck CCG_21CCSD */ 2247*=====================================================================* 2248 SUBROUTINE CCG_21CCSD (DELTA1, ISYRES, ! out: Delta vector 2249 & LISTA, IZETVA, ! inp: A resp. Zeta vec. 2250 & FockB, ISYMFB, ! inp: B transf. Fock Ma 2251 & T1AMPB, ISYMTB, ! inp: resp. B singles 2252 & T2AMPC, ISYMTC, ! inp: resp. C doubles 2253 & LISTC, IDLSTC, ! inp: resp. C 2254 & XIAJK, EIAJK, ! out: transf. integrals 2255 & ZIAJK, CIAJK, ! out: transf. multipl. 2256 & WORK, LWORK )! work space 2257*---------------------------------------------------------------------* 2258* 2259* Purpose: calculate CCSD contribution to DELTA1: 2260* 2261* DELTA1 = DELTA1 + <Zeta_2| [[tH^B,T2C], tau_1] |HF> 2262* 2263* 2264* symmetries/variables: 2265* 2266* ISYRES : result vector DELTA1 2267* ISYCTR : lagrangian multipliers (zeta vector) CTR1, CTR2 2268* ISYMTB : response amplitudes T1AMPB, (T2AmpB) 2269* ISYMTC : response amplitudes T2AMPC, (T1AMPC) 2270* 2271* Note: the double excitation response amplitudes T2AMPC are 2272* expected in packed lower triangular storage 2273* 2274* the double excitation part of lagrangian multiplier vector 2275* CTR2 and the (ia|jb) MO integrals will be read from disk 2276* 2277* required vectors of the size (1/2) V^2 O^2 : 2278* CCSD: (ia|jb), CTR2, T2AMPC 2279* 2280* (ia|jb) integrals are assumed to have the symmetry ISYMOP, 2281* which is transferred via COMMON /CCSDSYM/ 2282* 2283* Written by Christof Haettig, August 1996. 2284* 2285*=====================================================================* 2286#if defined (IMPLICIT_NONE) 2287 IMPLICIT NONE 2288#else 2289# include "implicit.h" 2290#endif 2291#include "priunit.h" 2292#include "maxorb.h" 2293#include "ccorb.h" 2294#include "ccsdsym.h" 2295#include "ccsdinp.h" 2296#include "r12int.h" 2297#include "ccr12int.h" 2298 2299* local parameters: 2300 CHARACTER MSGDBG*(19) 2301 PARAMETER (MSGDBG='[debug] CC_21CCSD> ') 2302#if defined (SYS_CRAY) 2303 REAL ONE, ZERO, HALF, TWO 2304#else 2305 DOUBLE PRECISION ONE, ZERO, HALF, TWO 2306#endif 2307 PARAMETER (ONE = 1.0d0, ZERO = 0.0d0, HALF = 0.5d0, TWO = 2.0d0) 2308 LOGICAL LOCDBG 2309 PARAMETER (LOCDBG = .FALSE.) 2310 2311 2312* input variables: 2313 CHARACTER*(*) LISTA, LISTC 2314 CHARACTER*(10) MODEL 2315 INTEGER ISYRES, ISYCTR, ISYMFB, ISYMTB, ISYMTC, LWORK, IZETVA 2316 INTEGER IDLSTC 2317 2318#if defined (SYS_CRAY) 2319 REAL DELTA1(*), FOCKB(*), T1AMPB(*), T2AMPC(*) 2320 REAL XIAJK(*), EIAJK(*), ZIAJK(*), CIAJK(*) 2321 REAL WORK(LWORK) 2322 REAL SWAP 2323#else 2324 DOUBLE PRECISION DELTA1(*), FOCKB(*), T1AMPB(*), T2AMPC(*) 2325 DOUBLE PRECISION XIAJK(*), EIAJK(*), ZIAJK(*), CIAJK(*) 2326 DOUBLE PRECISION WORK(LWORK) 2327 DOUBLE PRECISION SWAP 2328#endif 2329 2330* local variables: 2331 LOGICAL LMOINT 2332 INTEGER ITAIKJ(8,8), NTAIKJ(8) 2333 INTEGER KEND0, KEND1, KEND1B, KEND2, KEND4 2334 INTEGER LEND0, LEND1, LEND1B, LEND2, LEND4 2335 INTEGER KXCMAT, KYCMAT, KDXYMAT, KXIAJB, KCTR2 2336 INTEGER KSCR, KDUM, KMCINT, KKCINT 2337 INTEGER ISYTCTB, ISYOVOV, ISYMKC, ISYMMC, ISYMXY, ISYDXY 2338 INTEGER ISYTCIN, ISYMC, ISYMJKL, NCI, KOFF, ISYMS, ISYMKLJ 2339 INTEGER ISYMI, ISYMJ, MAXJ, NIJ, NJI, IOPT, IPOS 2340 INTEGER ISYM, ICOUNT, ISYMAIK, ISYMB, ISYCTB, NTOTC 2341 INTEGER ISYMK, ISYIAJ, ISYMAI, ISYMBJ, NBJ, KOFF1, KOFF2, KOFF3 2342 INTEGER KT2AMP, ISYMDLK, KOFFAIK, KOFFDKL, ISYMAL, ISYML, ISYMA 2343 INTEGER KCINT, ISYMDKL, KEND1C, LEND1C, NTOTAI, NTOTB, NTOTJKL 2344 INTEGER KWTINT, KVTINT, KWHINT, KVHINT, KEXC, NTOTDKL, NTOTA 2345 INTEGER IOPTTCME, KTR12AM, KVINT, IAN, LUNIT 2346 2347#if defined (SYS_CRAY) 2348 REAL DDOT 2349#else 2350 DOUBLE PRECISION DDOT 2351#endif 2352 2353 INTEGER ILSTSYM 2354 2355 CALL QENTER('CCG_21CCSD ') 2356 2357 IF (LOCDBG) THEN 2358 WRITE (LUPRI,*) MSGDBG, 'just entered this routine...' 2359 WRITE (LUPRI,*) MSGDBG, 'LWORK:',LWORK 2360 Call FLSHFO(LUPRI) 2361 END IF 2362 2363 DO ISYM = 1, NSYM 2364 ICOUNT = 0 2365 DO ISYMAIK = 1, NSYM 2366 ISYMJ = MULD2H(ISYMAIK,ISYM) 2367 ITAIKJ(ISYMAIK,ISYMJ) = ICOUNT 2368 ICOUNT = ICOUNT + NT2BCD(ISYMAIK)*NRHF(ISYMJ) 2369 END DO 2370 NTAIKJ(ISYM) = ICOUNT 2371 END DO 2372 2373*---------------------------------------------------------------------* 2374* begin: set start of work space and set & check symmetries 2375*---------------------------------------------------------------------* 2376 ISYOVOV = ISYMOP ! 2-el integrals 2377 ISYCTR = ILSTSYM(LISTA,IZETVA) ! lagrangian multipliers CTR2 2378 ISYTCTB = MULD2H(ISYMTC,ISYMTB) ! Wtilde, Vtilde 2379 ISYCTB = MULD2H(ISYCTR,ISYMTB) ! one-index transformed CTR2 2380 ISYMKC = MULD2H(ISYOVOV,ISYMTC) ! K^C intermediate 2381 ISYMMC = MULD2H(ISYCTR,ISYMTC) ! M^C intermediate 2382 2383 IF (MULD2H(ISYRES,ISYOVOV) .NE. MULD2H(ISYCTR,ISYTCTB)) THEN 2384 CALL QUIT('Symmetry mismatch in CCG_21CCSD.') 2385 END IF 2386 2387 KWTINT = 1 2388 KVTINT = KWTINT + NTAIKJ(ISYTCTB) 2389 KWHINT = KVTINT + NTAIKJ(ISYTCTB) 2390 KVHINT = KWHINT + NTAIKJ(ISYRES) 2391 KEND0 = KVHINT + NTAIKJ(ISYRES) 2392 2393 KKCINT = KEND0 2394 KMCINT = KKCINT + N3ORHF(ISYMKC) 2395 KEND0 = KMCINT + N3ORHF(ISYMMC) 2396 2397 LEND0 = LWORK - KEND0 2398 IF (LEND0 .LT. 0) THEN 2399 CALL QUIT('Insufficient work space in CCG_21CCSD. (0)') 2400 END IF 2401 2402 2403*---------------------------------------------------------------------* 2404* read (ia|jb) integrals from disk and square up: 2405*---------------------------------------------------------------------* 2406 KXIAJB = KEND0 2407 KEND1 = KXIAJB + NT2SQ(ISYOVOV) 2408 LEND1 = LWORK - KEND1 2409 2410 If (LEND1 .LT. NT2AM(ISYOVOV)) THEN 2411 CALL QUIT('Insufficient work space in CCG_21CCSD. (1a)') 2412 END IF 2413 2414* read & unpack (ia|jb) integrals: 2415 Call CCG_RDIAJB(WORK(KEND1), NT2AM(ISYOVOV)) 2416 2417CCN 2418C WRITE(LUPRI,*) '(ia|jb) integrals:' 2419C CALL CC_PRP(WORK(KDUM),WORK(KEND1),ISYOVOV,0,1) 2420C WRITE(LUPRI,*) 'T^C(el,dj) response amplitudes:' 2421C CALL CC_PRP(WORK(KDUM),T2AMPC,ISYMTC,0,1) 2422CCN 2423 2424 Call CC_T2SQ(WORK(KEND1), WORK(KXIAJB), ISYOVOV) 2425 2426*---------------------------------------------------------------------* 2427* precalculate K^C intermediate: K(i,jkl) = sum_ed (ie|kd) T^C(el,dj) 2428* this and the (ia|jk) integrals are the only intermediates, which 2429* require the (ia|jb) integrals as full square matrix: 2430*---------------------------------------------------------------------* 2431 ISYMKC = MULD2H(ISYOVOV,ISYMTC) ! K^C intermediate 2432 2433 Call CC_MI(WORK(KKCINT),WORK(KXIAJB), ISYOVOV, 2434 & T2AMPC,ISYMTC,WORK(KEND1),LEND1) 2435 2436*---------------------------------------------------------------------* 2437* for CCSD(R12) (and higher) add Sum_mn (V^\dagger)_(im,kn) T^C(ml,nj) 2438*---------------------------------------------------------------------* 2439 IF (CCR12.AND..NOT.(CCS.OR.CC2)) THEN 2440 KTR12AM = KEND1 2441 KVINT = KTR12AM + NTR12AM(ISYMTC) 2442 KEND1 = KVINT + NTR12AM(1) 2443 LEND1 = LWORK - KEND1 2444 IF (LEND1.LT.0) THEN 2445 CALL QUIT('Insufficient work space in CCG_21CCSD. (1b)') 2446 END IF 2447 2448 KDUM = - 99 999 999 ! dummy address 2449 IOPT = 32 2450 CALL CC_RDRSP(LISTC,IDLSTC,ISYMTC,IOPT,MODEL,WORK(KDUM), 2451 & WORK(KTR12AM)) 2452 2453 LUNIT = -1 2454 CALL GPOPEN(LUNIT,FCCR12V,'OLD',' ','UNFORMATTED', 2455 & KDUM,.FALSE.) 2456 9999 READ(LUNIT) IAN 2457 READ(LUNIT) (WORK(KVINT-1+I), I=1, NTR12AM(1)) 2458 IF (IAN.NE.IANR12) GOTO 9999 2459 CALL GPCLOSE(LUNIT,'KEEP') 2460 2461 IOPT = 2 2462 CALL CC_R12CV(WORK(KKCINT),WORK(KTR12AM),ISYMTC,WORK(KVINT), 2463 & 1,IOPT,WORK(KEND1),LEND1) 2464 2465 END IF 2466 2467 IF (.FALSE.) THEN 2468 WRITE (LUPRI,*) MSGDBG,'response K intermediate:' 2469 WRITE(*,'(5f12.6)') (WORK(KKCINT-1+I),I=1,N3ORHF(ISYMKC)) 2470 END IF 2471 2472*---------------------------------------------------------------------* 2473* calculate (ia|jk) integrals, stored as I^j(ai,k) 2474* and (ja|ik) integrals, stored as I^j(ai,k) 2475* finally replace (ia|jk) by (ia|jk) - 0.5*(ja|ik) 2476*---------------------------------------------------------------------* 2477 CALL CCG_TRANS4(XIAJK,ISYMTB,WORK(KXIAJB),ISYOVOV,T1AMPB,ISYMTB) 2478 2479 CALL CCG_SORT1(XIAJK,EIAJK,ISYMTB,0) 2480 2481 CALL DAXPY(NTAIKJ(ISYMTB),-HALF,EIAJK,1,XIAJK,1) 2482 2483*---------------------------------------------------------------------* 2484* read & unpack lagrangian multipliers and re-read packed integrals: 2485*---------------------------------------------------------------------* 2486* 1) squared multipliers & scratch space for packed multipliers: 2487 KCTR2 = KEND0 2488 KSCR = KCTR2 + NT2SQ(ISYCTR) 2489 KEND1B = KSCR + NT2AM(ISYCTR) 2490 LEND1B = LWORK - KEND1B 2491 2492* 2) squared multipliers & packed integrals: 2493 KXIAJB = KCTR2 + NT2SQ(ISYCTR) 2494 KEND1 = KXIAJB + NT2AM(ISYOVOV) 2495 LEND1 = LWORK - KEND1 2496 2497 If (LEND1B .LT. 0 .OR. LEND1 .LT.0) THEN 2498 CALL QUIT('Insufficient work space in CCG_21CCSD. (1b/1)') 2499 END IF 2500 2501* read packed lagrangian multipliers and square them up 2502 KDUM = - 99 999 999 ! dummy address 2503 2504 IOPT = 2 2505 Call CC_RDRSP(LISTA,IZETVA,ISYCTR,IOPT,MODEL, 2506 & WORK(KDUM),WORK(KSCR)) 2507 2508 Call CC_T2SQ(WORK(KSCR), WORK(KCTR2), ISYCTR) 2509 2510* read packed integrals (ia|jb) and calculate in place L(iajb): 2511 CALL CCG_RDIAJB(WORK(KXIAJB), NT2AM(ISYOVOV)) 2512 IOPTTCME = 1 2513 CALL CCSD_TCMEPK(WORK(KXIAJB),ONE,ISYOVOV,IOPTTCME) 2514 2515*---------------------------------------------------------------------* 2516* calculate X^C and Y^C intermediates and evaluate A+B and E terms: 2517*---------------------------------------------------------------------* 2518 ISYMXY = MULD2H(ISYCTR,ISYMTC) 2519 ISYDXY = MULD2H(ISYCTR,ISYTCTB) 2520 2521 KXCMAT = KEND1 2522 KYCMAT = KXCMAT + NMATIJ(ISYMXY) 2523 KDXYMAT = KYCMAT + NMATAB(ISYMXY) 2524 KSCR = KDXYMAT + NT1AM(ISYDXY) 2525 KEND2 = KSCR + NT1AM(ISYRES) 2526 LEND2 = LWORK - KEND2 2527 2528 IF (LEND2 .LT. 0) THEN 2529 CALL QUIT('Insufficient work space in CCG_21CCSD. (A+B)') 2530 END IF 2531 2532* calculate X^C & Y^C intermediate: 2533 Call CC_XI(WORK(KXCMAT),WORK(KCTR2), ISYCTR, 2534 & T2AMPC,ISYMTC,WORK(KEND2),LEND2) 2535 2536 Call CC_YI(WORK(KYCMAT),WORK(KCTR2), ISYCTR, 2537 & T2AMPC,ISYMTC,WORK(KEND2),LEND2) 2538 2539 IF (LOCDBG) THEN 2540 WRITE (LUPRI,*) MSGDBG,'response X intermediate:' 2541 WRITE (LUPRI,'(5f12.6)') (WORK(KXCMAT-1+I),I=1,NMATIJ(ISYMXY)) 2542 WRITE (LUPRI,*) MSGDBG,'response Y intermediate:' 2543 WRITE (LUPRI,'(2f12.6)') (WORK(KYCMAT-1+I),I=1,NMATAB(ISYMXY)) 2544 END IF 2545 2546* calculate XY^C: XY^C_ij = X^C_ji, XY^C_bd = -Y^C_bd 2547* 1.) transpose X^C intermediate 2548 DO ISYMI = 1, NSYM 2549 ISYMJ = MULD2H(ISYMI,ISYMXY) 2550 IF (ISYMJ .LE. ISYMI) THEN 2551 DO I = 1, NRHF(ISYMI) 2552 MAXJ = NRHF(ISYMJ) 2553 IF (ISYMJ .EQ. ISYMI) MAXJ = I-1 2554 DO J = 1, MAXJ 2555 NIJ = IMATIJ(ISYMI,ISYMJ) + NRHF(ISYMI)*(J-1) + I 2556 NJI = IMATIJ(ISYMJ,ISYMI) + NRHF(ISYMJ)*(I-1) + J 2557 SWAP = WORK(KXCMAT-1+NIJ) 2558 WORK(KXCMAT-1+NIJ) = WORK(KXCMAT-1+NJI) 2559 WORK(KXCMAT-1+NJI) = SWAP 2560 END DO 2561 END DO 2562 END IF 2563 END DO 2564 2565* 2.) multiply Y^C intermediate with -1: 2566 Call DSCAL(NMATAB(ISYMXY), -ONE, WORK(KYCMAT), 1) 2567 2568*.....................................................................* 2569* A+B term: 2570*.....................................................................* 2571* calculate effective density matrix DXY = [XY^C,T1^B] 2572 IOPT = 1 2573 Call CCG_1ITRVO(WORK(KDXYMAT),ISYDXY,WORK(KXCMAT),WORK(KYCMAT), 2574 & ISYMXY,T1AMPB,ISYMTB,IOPT) 2575 2576* contract with L(ia|jb): 2577 Call CCG_LXD(WORK(KSCR), ISYRES, WORK(KDXYMAT), ISYDXY, 2578 & WORK(KXIAJB),ISYOVOV,0) 2579 2580* add contribution to DELTA1 2581 Call DAXPY (NT1AM(ISYRES), ONE, WORK(KSCR),1, DELTA1, 1) 2582 2583 IF (LOCDBG) THEN 2584 WRITE (LUPRI,*) MSGDBG,'T1AMPB:' 2585 WRITE (LUPRI,'(5f12.6)') (T1AMPB(I),I=1,NT1AM(ISYMTB)) 2586 WRITE (LUPRI,*) MSGDBG,'DXY intermediate:' 2587 WRITE (LUPRI,'(5f12.6)') (WORK(KDXYMAT-1+I),I=1,NT1AM(ISYRES)) 2588 WRITE (LUPRI,*) MSGDBG, 'A+B term:' 2589 CALL CC_PRP(WORK(KSCR),WORK,ISYRES,1,0) 2590 END IF 2591 2592*.....................................................................* 2593* E term: 2594*.....................................................................* 2595 2596* do one-index transformation of XY^C with F^B: 2597 IOPT = 2 2598 Call CCG_1ITRVO(WORK(KSCR),ISYRES, 2599 & WORK(KXCMAT),WORK(KYCMAT),ISYMXY, 2600 & FockB, ISYMFB, IOPT ) 2601 2602* add contribution to DELTA1 2603 Call DAXPY (NT1AM(ISYRES), ONE, WORK(KSCR),1, DELTA1, 1) 2604 2605 IF (LOCDBG) THEN 2606 WRITE (LUPRI,*) MSGDBG,'FockB:' 2607 WRITE (LUPRI,'(5f12.6)') (FockB(I),I=1,NT1AM(ISYMFB)) 2608 WRITE (LUPRI,*) MSGDBG, 'E term:' 2609 CALL CC_PRP(WORK(KSCR),WORK,ISYRES,1,0) 2610 WRITE (LUPRI,*) MSGDBG, 'Delta1 now:' 2611 CALL CC_PRP(DELTA1,WORK,ISYRES,1,0) 2612 END IF 2613 2614*---------------------------------------------------------------------* 2615 2616*---------------------------------------------------------------------* 2617* precalculate M^C intermediate: M(i,jkl) = sum_ed Z(ei,dk) T^C(el,dj) 2618* (overwrites X^C and Y^C intermediates) 2619*---------------------------------------------------------------------* 2620 ISYMMC = MULD2H(ISYCTR,ISYMTC) 2621 2622 Call CC_MI(WORK(KMCINT),WORK(KCTR2),ISYCTR, 2623 & T2AMPC,ISYMTC,WORK(KEND1),LEND1) 2624 2625 IF (.FALSE.) THEN 2626 WRITE (LUPRI,*) MSGDBG,'response M intermediate:' 2627 WRITE(LUPRI,'(5f12.6)') (WORK(KMCINT-1+I),I=1,N3ORHF(ISYMMC)) 2628 END IF 2629 2630*---------------------------------------------------------------------* 2631* calculate Zeta(ia|jk), stored as Z^j(ai,k) 2632* and Zeta(ja|ik), stored as C^j(ai,k) 2633* finally replace Zeta(ja|ik) by 2 Zeta(ja|ik) + Zeta(ia|jk) 2634*---------------------------------------------------------------------* 2635 CALL CCG_TRANS4(ZIAJK,ISYCTB,WORK(KCTR2),ISYCTR,T1AMPB,ISYMTB) 2636 2637 CALL CCG_SORT1(ZIAJK,CIAJK,ISYCTB,0) 2638 2639 CALL DAXPY(NTAIKJ(ISYCTB),HALF,ZIAJK,1,CIAJK,1) 2640 2641*---------------------------------------------------------------------* 2642* calculate Wtilde and What intermediates: 2643* Wtilde(dl,k;j) = sum_ia [2(ia|jk) - (ja|ik)] [2T(ai,dl)-T(al,di)] 2644* What(dl,k;j) = sum_ia Z(ia|jk) [2T(ai,dl)-T(al,di)] 2645* calculate Vtilde and Vhat intermediates: 2646* Vhat(dl,k;j) = sum_ia [2Z(ja|ik)+Z(ia|jk)] T(al,di) 2647* Vtilde(dl,k;j) = sum_ia (ja|ik) T(al,di) 2648*---------------------------------------------------------------------* 2649 KT2AMP = KEND0 2650 KEND1 = KT2AMP + NT2SQ(ISYMTC) 2651 LEND1 = LWORK - KEND1 2652 2653 IF (LEND1 .LT. 0) THEN 2654 CALL QUIT('Insufficient work space in CCG_21CCSD. (1b)') 2655 END IF 2656 2657* square up T2 amplitudes and calculate in place 2 Cou - Exc comb.: 2658 CALL CC_T2SQ(T2AMPC,WORK(KT2AMP),ISYMTC) 2659 CALL CCRHS_T2TR(WORK(KT2AMP),WORK(KEND1),LEND1,ISYMTC) 2660 2661* Wtilde intermediate: 2662 DO ISYMJ = 1, NSYM 2663 ISYMAIK = MULD2H(ISYMTB,ISYMJ) 2664 ISYMDLK = MULD2H(ISYMAIK,ISYMTC) 2665 2666 DO J = 1, NRHF(ISYMJ) 2667 KOFFAIK = 1+ITAIKJ(ISYMAIK,ISYMJ)+NT2BCD(ISYMAIK)*(J-1) 2668 KOFFDKL = ITAIKJ(ISYMDLK,ISYMJ) + NT2BCD(ISYMDLK)*(J-1) 2669 2670 IOPT = 1 2671 CALL CC_ZWVI(WORK(KWTINT+KOFFDKL),WORK(KT2AMP),ISYMTC, 2672 & XIAJK(KOFFAIK),ISYMAIK,WORK(KEND1),LEND1,IOPT) 2673 END DO 2674 END DO 2675 2676* What intermediate: 2677 DO ISYMJ = 1, NSYM 2678 ISYMAIK = MULD2H(ISYCTB,ISYMJ) 2679 ISYMDLK = MULD2H(ISYMAIK,ISYMTC) 2680 2681 DO J = 1, NRHF(ISYMJ) 2682 KOFFAIK = 1+ITAIKJ(ISYMAIK,ISYMJ)+NT2BCD(ISYMAIK)*(J-1) 2683 KOFFDKL = ITAIKJ(ISYMDLK,ISYMJ) + NT2BCD(ISYMDLK)*(J-1) 2684 2685 IOPT = 1 2686 CALL CC_ZWVI(WORK(KWHINT+KOFFDKL),WORK(KT2AMP),ISYMTC, 2687 & ZIAJK(KOFFAIK),ISYMAIK,WORK(KEND1),LEND1,IOPT) 2688 END DO 2689 END DO 2690 2691 2692* recover squared T2 amplitudes and transpose occupied indeces: 2693 CALL CC_T2SQ(T2AMPC,WORK(KT2AMP),ISYMTC) 2694 CALL CCSD_T2TP(WORK(KT2AMP),WORK(KEND1),LEND1,ISYMTC) 2695 2696* Vtilde intermediate: 2697 DO ISYMJ = 1, NSYM 2698 ISYMAIK = MULD2H(ISYMTB,ISYMJ) 2699 ISYMDLK = MULD2H(ISYMAIK,ISYMTC) 2700 2701 DO J = 1, NRHF(ISYMJ) 2702 KOFFAIK = 1+ITAIKJ(ISYMAIK,ISYMJ)+NT2BCD(ISYMAIK)*(J-1) 2703 KOFFDKL = ITAIKJ(ISYMDLK,ISYMJ) + NT2BCD(ISYMDLK)*(J-1) 2704 2705 IOPT = 1 2706 CALL CC_ZWVI(WORK(KVTINT+KOFFDKL), WORK(KT2AMP), ISYMTC, 2707 & EIAJK(KOFFAIK),ISYMAIK,WORK(KEND1),LEND1,IOPT) 2708 END DO 2709 END DO 2710 2711* Vhat intermediate: 2712 DO ISYMJ = 1, NSYM 2713 ISYMAIK = MULD2H(ISYCTB,ISYMJ) 2714 ISYMDLK = MULD2H(ISYMAIK,ISYMTC) 2715 2716 DO J = 1, NRHF(ISYMJ) 2717 KOFFAIK = 1+ITAIKJ(ISYMAIK,ISYMJ)+NT2BCD(ISYMAIK)*(J-1) 2718 KOFFDKL = ITAIKJ(ISYMDLK,ISYMJ) + NT2BCD(ISYMDLK)*(J-1) 2719 2720 IOPT = 1 2721 CALL CC_ZWVI(WORK(KVHINT+KOFFDKL), WORK(KT2AMP), ISYMTC, 2722 & CIAJK(KOFFAIK),ISYMAIK,WORK(KEND1),LEND1,IOPT) 2723 END DO 2724 END DO 2725 2726*---------------------------------------------------------------------* 2727* reorganization of work space: 2728* KKCINT : K^C intermediate 2729* KMCINT : M^C intermediate 2730* KCTR2 : zeta vector (packed) 2731* KXIAJB : (ia|jb) integrals (packed) 2732* KEND2 : free work space 2733*---------------------------------------------------------------------* 2734 KCTR2 = KEND0 2735 KXIAJB = KCTR2 + NT2AM(ISYCTR) 2736 KEND2 = KXIAJB + NT2AM(ISYOVOV) 2737 LEND2 = LWORK - KEND2 2738 2739 IF (LEND2 .LT. 0) THEN 2740 CALL QUIT('Insufficient work space in CCG_21CCSD.') 2741 END IF 2742 2743* read (ia|jb) integrals: 2744 Call CCG_RDIAJB(WORK(KXIAJB), NT2AM(ISYOVOV)) 2745 2746* read packed lagrangian multipliers and square them up 2747 IOPT = 2 2748 Call CC_RDRSP(LISTA,IZETVA,ISYCTR,IOPT,MODEL, 2749 & WORK(KDUM),WORK(KCTR2)) 2750 2751*---------------------------------------------------------------------* 2752* H term: contract Wtilde/Vtilde interm. with the Lagrangian multipl.: 2753*---------------------------------------------------------------------* 2754 DO ISYMA = 1, NSYM 2755 2756 ISYMDKL = MULD2H(ISYCTR,ISYMA) 2757 ISYMI = MULD2H(ISYTCTB,ISYMDKL) 2758 2759 KCINT = KEND2 2760 KEXC = KCINT + NT2BCD(ISYMDKL) 2761 KEND1C = KEXC + NT2BCD(ISYMDKL) 2762 LEND1C = LWORK - KEND1C 2763 2764 IF (LEND1C .LT. 0) THEN 2765 CALL QUIT('Insufficient work space in CCG_21CCSD. (1c)') 2766 END IF 2767 2768 DO A = 1, NVIR(ISYMA) 2769 2770 CALL CCG_TI(WORK(KCINT),ISYMDKL,WORK(KCTR2),ISYCTR,A,ISYMA) 2771 2772 CALL DCOPY(NT2BCD(ISYMDKL),WORK(KCINT),1,WORK(KEXC),1) 2773 CALL CCLT_P21I(WORK(KEXC),ISYMDKL,WORK(KEND1C),LEND1C, 2774 & IT2BCD, NT2BCD, IT1AM, NT1AM, NVIR) 2775 CALL DAXPY(NT2BCD(ISYMDKL),HALF,WORK(KCINT),1,WORK(KEXC),1) 2776 2777 KOFF1 = ITAIKJ(ISYMDKL,ISYMI) 2778 KOFF2 = IT1AM(ISYMA,ISYMI) + A 2779 NTOTDKL = MAX(NT2BCD(ISYMDKL),1) 2780 NTOTA = MAX(NVIR(ISYMA),1) 2781 2782 CALL DGEMV('T',NT2BCD(ISYMDKL),NRHF(ISYMI), 2783 & -ONE,WORK(KWTINT+KOFF1),NTOTDKL,WORK(KCINT),1, 2784 & ONE,DELTA1(KOFF2),NTOTA) 2785 2786 CALL DGEMV('T',NT2BCD(ISYMDKL),NRHF(ISYMI), 2787 & ONE,WORK(KVTINT+KOFF1),NTOTDKL,WORK(KEXC),1, 2788 & ONE,DELTA1(KOFF2),NTOTA) 2789 2790 END DO 2791 END DO 2792 2793 2794 IF (LOCDBG) THEN 2795 WRITE (LUPRI,*) MSGDBG, 'Delta1 after H term in CCG_21CCSD:' 2796 Call CC_PRP(DELTA1,WORK,ISYRES,1,0) 2797 END IF 2798 2799*---------------------------------------------------------------------* 2800* I term: contract What/Vhat intermediates with the integrals: 2801*---------------------------------------------------------------------* 2802 DO ISYMA = 1, NSYM 2803 2804 ISYMDKL = MULD2H(ISYOVOV,ISYMA) 2805 ISYMI = MULD2H(ISYRES,ISYMDKL) 2806 2807 KCINT = KEND2 2808 KEXC = KCINT + NT2BCD(ISYMDKL) 2809 KEND1C = KEXC + NT2BCD(ISYMDKL) 2810 LEND1C = LWORK - KEND1C 2811 2812 IF (LEND1C .LT. 0) THEN 2813 CALL QUIT('Insufficient work space in CCG_21CCSD. (1c)') 2814 END IF 2815 2816 DO A = 1, NVIR(ISYMA) 2817 2818 CALL CCG_TI(WORK(KCINT),ISYMDKL,WORK(KXIAJB),ISYOVOV,A,ISYMA) 2819 2820 CALL DCOPY(NT2BCD(ISYMDKL),WORK(KCINT),1,WORK(KEXC),1) 2821 CALL CCLT_P21I(WORK(KEXC),ISYMDKL,WORK(KEND1C),LEND1C, 2822 & IT2BCD, NT2BCD, IT1AM, NT1AM, NVIR) 2823 CALL DAXPY(NT2BCD(ISYMDKL),-HALF,WORK(KEXC),1,WORK(KCINT),1) 2824 2825 KOFF1 = ITAIKJ(ISYMDKL,ISYMI) 2826 KOFF2 = IT1AM(ISYMA,ISYMI) + A 2827 NTOTDKL = MAX(NT2BCD(ISYMDKL),1) 2828 NTOTA = MAX(NVIR(ISYMA),1) 2829 2830 CALL DGEMV('T',NT2BCD(ISYMDKL),NRHF(ISYMI), 2831 & -ONE,WORK(KWHINT+KOFF1),NTOTDKL,WORK(KCINT),1, 2832 & ONE,DELTA1(KOFF2),NTOTA) 2833 2834 CALL DGEMV('T',NT2BCD(ISYMDKL),NRHF(ISYMI), 2835 & ONE,WORK(KVHINT+KOFF1),NTOTDKL,WORK(KEXC),1, 2836 & ONE,DELTA1(KOFF2),NTOTA) 2837 2838 END DO 2839 END DO 2840 2841 IF (LOCDBG) THEN 2842 WRITE (LUPRI,*) MSGDBG, 'Delta1 after I term in CCG_21CCSD:' 2843 Call CC_PRP(DELTA1,WORK,ISYRES,1,0) 2844 END IF 2845 2846*---------------------------------------------------------------------* 2847* resort transformed integrals and multipliers for G and F term: 2848*---------------------------------------------------------------------* 2849 2850 CALL DAXPY(NTAIKJ(ISYMTB),HALF,EIAJK,1,XIAJK,1) 2851 CALL CCG_SORT1(XIAJK,EIAJK,ISYMTB,4) 2852 2853 CALL CCG_SORT1(ZIAJK,CIAJK,ISYCTB,4) 2854 2855*---------------------------------------------------------------------* 2856* loop over symmetry blocks and calculate G and F term contributions: 2857*---------------------------------------------------------------------* 2858 IF (LOCDBG) THEN 2859 WRITE (LUPRI,*) MSGDBG, 'Delta1 before G and F term:' 2860 Call CC_PRP(DELTA1,WORK,ISYRES,1,0) 2861 END IF 2862 DO ISYMC = 1, NSYM 2863 2864 ISYMI = MULD2H(ISYMC,ISYRES) 2865 NCI = IT1AM(ISYMC,ISYMI) + 1 2866 NTOTC = MAX(NVIR(ISYMC),1) 2867 2868C -------------------------------------- 2869C G term: contract (k l^B|j c) integrals 2870C with M^C(i,klj) intermediate: 2871C -------------------------------------- 2872 ISYMJKL = MULD2H(ISYMI,ISYMMC) 2873 2874 KOFF1 = KMCINT + I3ORHF(ISYMJKL,ISYMI) 2875 KOFF2 = ISJIKA(ISYMJKL,ISYMC) + 1 2876 NTOTJKL = MAX(NMAIJK(ISYMJKL),1) 2877 2878 CALL DGEMM('N','N',NVIR(ISYMC),NRHF(ISYMI),NMAIJK(ISYMJKL), 2879 & ONE,EIAJK(KOFF2),NTOTC,WORK(KOFF1),NTOTJKL, 2880 & ONE,DELTA1(NCI), NTOTC) 2881 2882C ---------------------------------------- 2883C F term: contract S^{B,c}(klj) multiplier 2884C with K^C(i,klj) intermediate 2885C ---------------------------------------- 2886 ISYMS = MULD2H(ISYCTR,ISYMTB) 2887 ISYMKLJ = MULD2H(ISYMC,ISYMS) 2888 2889 KOFF1 = KKCINT + I3ORHF(ISYMKLJ,ISYMI) 2890 KOFF2 = ISJIKA(ISYMKLJ,ISYMC) + 1 2891 NTOTJKL = MAX(NMAIJK(ISYMKLJ),1) 2892 2893 CALL DGEMM('N','N',NVIR(ISYMC),NRHF(ISYMI),NMAIJK(ISYMKLJ), 2894 & ONE,CIAJK(KOFF2),NTOTC,WORK(KOFF1),NTOTJKL, 2895 & ONE,DELTA1(NCI), NTOTC) 2896CCN 2897C WRITE(LUPRI,*) 'K^C(i,klj):' 2898C CALL OUTPUT(WORK(KOFF1),1,NTOTJKL,1,NRHF(ISYMI), 2899C & NTOTJKL,NRHF(ISYMI),1,LUPRI) 2900C WRITE(LUPRI,*) 'CIAJK:' 2901C CALL OUTPUT(CIAJK(KOFF2),1,NVIR(ISYMC), 2902C & 1,NMAIJK(ISYMKLJ),NVIR(ISYMC),NMAIJK(ISYMKLJ), 2903C & 1,LUPRI) 2904C WRITE (LUPRI,*) MSGDBG, 'after Symmetry block:',ISYMC 2905C Call CC_PRP(DELTA1,WORK,ISYRES,1,0) 2906CCN 2907 END DO ! ISYMC 2908 2909CCN 2910C WRITE(LUPRI,*) 'S^{B,c}(klj): ' 2911C DO ISYMC = 1, NSYM 2912C ISYMS = MULD2H(ISYCTR,ISYMTB) 2913C ISYMKLJ = MULD2H(ISYMC,ISYMS) 2914C WRITE(LUPRI,*) 'Symmetry block: ',ISYMKLJ,ISYMC 2915C CALL OUTPUT(CIAJK(1+ISJIKA(ISYMKLJ,ISYMC)),1,NVIR(ISYMC), 2916C & 1,NMAIJK(ISYMKLJ),NVIR(ISYMC),NMAIJK(ISYMKLJ), 2917C & 1,LUPRI) 2918C END DO 2919CCN 2920 2921*---------------------------------------------------------------------* 2922* print debug output and return: 2923*---------------------------------------------------------------------* 2924 IF (LOCDBG) THEN 2925 WRITE (LUPRI,*) MSGDBG, 'Delta1 at the end of CCG_21CCSD:' 2926 Call CC_PRP(DELTA1,WORK,ISYRES,1,0) 2927 END IF 2928 2929 CALL QEXIT('CCG_21CCSD ') 2930 2931 RETURN 2932 END 2933*---------------------------------------------------------------------* 2934c/* Deck CCG_21CD */ 2935*=====================================================================* 2936 SUBROUTINE CCG_21CD( RESVEC,ISYRES,CTR2,ISYCTR,XLAMDH0,XLAMDP0, 2937 & XLAMDHA,XLAMDPA,ISYMLA, 2938 & XLAMDHB,XLAMDPB,ISYMLB, 2939 & DSRHFC,ISYMXC,IDEL,ISYDEL, 2940 & LSAME, FACTOR,WORK,LWORK) 2941*---------------------------------------------------------------------* 2942* 2943* Purpose: to calculate G and K matrix contributions which are 2944* derivatives of the 21C and 21D contribution to RHO1 2945* (the left hand side transformed trial vector). 2946* for that this routine has to be called 3 times with 2947* appropriate permutations of the XLAMDH/XLAMDP matrices 2948* 2949* 2950* 2951* 21C contribution: 2952* 2953* RESVEC(a i) = RESVEC(a i) - FACTOR * XLAMDP0(del i) * sum_{k j c} 2954* [(del j^C|c^B k^A) + (del j^C|c^A k^B)] * CTR2(a j c k) 2955* 2956* 2957* 21D contribution: 2958* 2959* RESVEC(a i) = RESVEC(a i) + sum_{bet} FACTOR * XLAMDH0(bet a) * 2960* sum_{j alp} (del j^C|alp bet) * 2961* [ CTR2(alp^B i j del^A) + CTR2(alp^A i j delB)] 2962* 2963* 2964* For LSAME=.true. XLAMDHA & XLAMDHB and XLAMDPA & XLAMDPB are assumed 2965* to be the same --> some intermediates are multiplied by 2 instead 2966* of calculating them again with A and B exchanged... 2967* 2968* symmetries & variables: 2969* 2970* ISYRES : result vector RESVEC 2971* ISYCTR : lagrangian multipliers (zeta vector) CTR2 2972* ISYMXC : symmetry of DSRHFC 2973* ISYDEL : symmetry of SAO basis function IDEL 2974* ISYMLA : Lambda matrices XLAMDHA XLAMDPA 2975* ISYMLB : Lambda matrices XLAMDHB, XLAMDPB 2976* Lambda matrices XLAMDH0, XLAMDP0 assumed tot.sym. 2977* 2978* used for CC2 part of CCQR_G 2979* 2980* Christof Haettig, August 1996, revised in August 1998 2981*=====================================================================* 2982#if defined (IMPLICIT_NONE) 2983 IMPLICIT NONE 2984#else 2985# include "implicit.h" 2986#endif 2987#include "dummy.h" 2988#include "priunit.h" 2989#include "ccsdsym.h" 2990#include "ccorb.h" 2991 2992 2993#if defined (SYS_CRAY) 2994 REAL ZERO, ONE, TWO, FACT, FACTOR 2995#else 2996 DOUBLE PRECISION ZERO, ONE, TWO, FACT, FACTOR 2997#endif 2998 PARAMETER (ZERO = 0.0d0, ONE = 1.0d0, TWO = 2.0d0) 2999 3000 LOGICAL LSAME 3001 INTEGER ISYRES, ISYCTR, ISYDEL, ISYMLA, ISYMLB, ISYMXC 3002 INTEGER IDEL, LWORK 3003 3004#if defined (SYS_CRAY) 3005 REAL RESVEC(*) ! dimension (NT1AM(ISYRES)) 3006 REAL CTR2(*) ! dimension (NT2AM(ISYCTR)) 3007 REAL DSRHFC(*) 3008 REAL XLAMDH0(*) ! dimension (NLAMDT) 3009 REAL XLAMDP0(*) ! dimension (NLAMDT) 3010 REAL XLAMDHA(*) ! dimension (NGLMDT(ISYMLA)) 3011 REAL XLAMDPA(*) ! dimension (NGLMDT(ISYMLA)) 3012 REAL XLAMDHB(*) ! dimension (NGLMDT(ISYMLB)) 3013 REAL XLAMDPB(*) ! dimension (NGLMDT(ISYMLB)) 3014 REAL WORK(LWORK) 3015 REAL DDOT 3016#else 3017 DOUBLE PRECISION RESVEC(*) ! dimension (NT1AM(ISYRES)) 3018 DOUBLE PRECISION CTR2(*) ! dimension (NT2AM(ISYCTR)) 3019 DOUBLE PRECISION DSRHFC(*) 3020 DOUBLE PRECISION XLAMDH0(*) ! dimension (NLAMDT) 3021 DOUBLE PRECISION XLAMDP0(*) ! dimension (NLAMDT) 3022 DOUBLE PRECISION XLAMDHA(*) ! dimension (NGLMDT(ISYMLA)) 3023 DOUBLE PRECISION XLAMDPA(*) ! dimension (NGLMDT(ISYMLA)) 3024 DOUBLE PRECISION XLAMDHB(*) ! dimension (NGLMDT(ISYMLB)) 3025 DOUBLE PRECISION XLAMDPB(*) ! dimension (NGLMDT(ISYMLB)) 3026 DOUBLE PRECISION WORK(LWORK) 3027 DOUBLE PRECISION DDOT 3028#endif 3029 3030 INTEGER ISYLALB, ISYMAO, ISYMMO 3031 INTEGER KDRES, KCRES, KEND0, LEND0, KEND1, LEND1 3032 INTEGER ISYCTA, KCTR2A, IOPT, KLAM0, KOFFR, ISYMBET, NTOTA 3033 INTEGER ISYMJ, ISYALBE, KSQC, KEND2, LEND2 3034 INTEGER KOFFC, KOFFB, ISYMALP, ISYMC, ISYMK, KHALF, KMO 3035 INTEGER KOFFIC, KOFFLB, KOFFLA 3036 INTEGER ISYMCI, ISYMI, KOFFZ, KOFFSV, NTOTC, NTOTBE, NTOTAL 3037 INTEGER ISYMCK, ISYMAJ, ISYMA, NAJ, KCTR2, ISYMD, KRESV 3038 INTEGER NTOTB, ISYMAI, KOFFZC, KOFFZB, ISYCTAB 3039 INTEGER KCTR2AB, NHALF, ISYCTB 3040 INTEGER KMOINT 3041 3042 3043C INTEGER INDEX 3044C INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 3045 3046 CALL QENTER('CCG_21CD') 3047 3048* set & check symmetries: 3049 ISYLALB = MULD2H(ISYMLA,ISYMLB) ! Lambda C x Lambda B 3050 ISYMAO = MULD2H(ISYMXC,ISYDEL) ! AO integrals 3051 ISYMMO = MULD2H(ISYMAO,ISYLALB) ! MO integrals 3052 3053 IF (ISYRES .NE. MULD2H(ISYMMO,ISYCTR)) THEN 3054 CALL QUIT('Symmetry mismatch in CCG_21CD') 3055 END IF 3056 3057* allocate scratch arrays for half-transformed result vectors: 3058 KDRES = 1 3059 KCRES = KDRES + NT1AO(ISYRES) 3060 KEND1 = KCRES + NVIRT 3061 LEND1 = LWORK - KEND1 3062 3063 Call DZERO(WORK(KDRES),NT1AO(ISYRES)) 3064 Call DZERO(WORK(KCRES),NVIRT) 3065 3066 IF (LEND1 .LT. 0) THEN 3067 CALL QUIT('Insufficient work space in CCG_21CD (0).') 3068 END IF 3069 3070*---------------------------------------------------------------------* 3071* calculate one-index back-transformed CTR2 vector using XLAMDPA: 3072*---------------------------------------------------------------------* 3073 ISYCTA = MULD2H(ISYDEL,MULD2H(ISYCTR,ISYMLA)) 3074 ISYCTB = MULD2H(ISYDEL,MULD2H(ISYCTR,ISYMLB)) 3075 ISYCTAB = MULD2H(ISYCTA,ISYMLB) 3076 3077 KCTR2A = KEND1 3078 KCTR2AB = KCTR2A + MAX(NT2BCD(ISYCTA),NT2BCD(ISYCTB)) 3079 KEND1 = KCTR2AB + NT2BGD(ISYCTAB) 3080 LEND1 = LWORK - KEND1 3081 3082 IF (LEND1 .LT. 0) THEN 3083 CALL QUIT('Insufficient work space in CCG_21CD (1b).') 3084 END IF 3085 3086 D = IDEL - IBAS(ISYDEL) 3087 3088* back transform the virtual indeces of CTR2 with XLAMDPA/XLAMDPB: 3089 IOPT = 1 3090 Call CC_T2AO(CTR2, XLAMDPA, ISYMLA, WORK(KCTR2A), 3091 & WORK(KEND1), LEND1, IDEL, ISYDEL, ISYCTR, IOPT) 3092 3093 IOPT = 0 3094 CALL CC_BFAO(WORK(KCTR2AB),WORK(KCTR2A),ISYCTA,XLAMDPB,ISYMLB, 3095 & D,ISYDEL,DUMMY,IDUMMY,DUMMY,IDUMMY,ZERO,IOPT) 3096 3097 IF (LSAME) THEN 3098 CALL DSCAL(NT2BGD(ISYCTAB),TWO,WORK(KCTR2AB),1) 3099 ELSE 3100 IOPT = 1 3101 Call CC_T2AO(CTR2, XLAMDPB, ISYMLB, WORK(KCTR2A), 3102 & WORK(KEND1), LEND1, IDEL, ISYDEL, ISYCTR, IOPT) 3103 IOPT = 0 3104 CALL CC_BFAO(WORK(KCTR2AB),WORK(KCTR2A),ISYCTB,XLAMDPA,ISYMLA, 3105 & D,ISYDEL,DUMMY,IDUMMY,DUMMY,IDUMMY,ONE,IOPT) 3106 END IF 3107 3108*---------------------------------------------------------------------* 3109* start loop over occupied index j: 3110*---------------------------------------------------------------------* 3111 DO ISYMJ = 1, NSYM 3112 ISYALBE = MULD2H(ISYMJ,ISYMXC) 3113 ISYMCK = MULD2H(ISYALBE,MULD2H(ISYMLA,ISYMLB)) 3114 3115 KSQC = KEND1 3116 KMOINT = KSQC + N2BST(ISYALBE) 3117 KEND2 = KMOINT + NT1AM(ISYMCK) 3118 LEND2 = LWORK - KEND2 3119 3120 IF (LEND2 .LT. 0) THEN 3121 CALL QUIT('Insufficient work space in CCG_21CD (2).') 3122 END IF 3123 3124 3125 DO J = 1, NRHF(ISYMJ) 3126 3127 CALL DZERO(WORK(KMOINT),NT1AM(ISYMCK)) 3128 3129C --------------------------------------------- 3130C unpack integral distribution (al be | j del): 3131C --------------------------------------------- 3132 KOFFC = IDSRHF(ISYALBE,ISYMJ) + NNBST(ISYALBE)*(J-1) + 1 3133 3134 Call CCSD_SYMSQ (DSRHFC(KOFFC), ISYALBE, WORK(KSQC)) 3135 3136 3137C ---------------------------------------------------- 3138C loop over symmetry blocks of unpacked distributions: 3139C ---------------------------------------------------- 3140 DO ISYMBET = 1, NSYM 3141 ISYMALP = MULD2H(ISYMBET,ISYALBE) 3142 3143*---------------------------------------------------------------------* 3144* calculate (c^B k^A |j^C del) + (c^A k^B |j^c del) integrals: 3145*---------------------------------------------------------------------* 3146 ISYMC = MULD2H(ISYMALP,ISYMLB) 3147 ISYMK = MULD2H(ISYMBET,ISYMLA) 3148 3149 NHALF = NBAS(ISYMALP) * NRHF(ISYMK) 3150 3151 IF (LEND2 .LT. NHALF) THEN 3152 CALL QUIT('Insufficient work space in CCG_21CD (3).') 3153 END IF 3154 3155* first XLAMDPB(alp c) * XLAMDHA(bet k) * (alp bet|j^C del): 3156 KOFFIC = KSQC + IAODIS(ISYMALP,ISYMBET) 3157 KOFFLA = IGLMRH(ISYMBET,ISYMK) + 1 3158 KOFFLB = IGLMVI(ISYMALP,ISYMC) + 1 3159 KMO = KMOINT + IT1AM(ISYMC,ISYMK) 3160 KHALF = KEND2 3161 3162 NTOTAL = MAX(NBAS(ISYMALP),1) 3163 NTOTBE = MAX(NBAS(ISYMBET),1) 3164 NTOTC = MAX(NVIR(ISYMC),1) 3165 3166 Call DGEMM('N','N',NBAS(ISYMALP),NRHF(ISYMK),NBAS(ISYMBET), 3167 & ONE, WORK(KOFFIC),NTOTAL, XLAMDHA(KOFFLA),NTOTBE, 3168 & ZERO, WORK(KHALF),NTOTAL) 3169 3170 Call DGEMM('T','N',NVIR(ISYMC), NRHF(ISYMK), NBAS(ISYMALP), 3171 & ONE, XLAMDPB(KOFFLB),NTOTAL, WORK(KHALF),NTOTAL, 3172 & ONE, WORK(KMO),NTOTC ) 3173 3174 3175* add XLAMDPA(alp c) * XLAMDHB(bet k) * (alp bet|j^C del): 3176 IF (.NOT. LSAME) THEN 3177 3178 ISYMC = MULD2H(ISYMALP,ISYMLA) 3179 ISYMK = MULD2H(ISYMBET,ISYMLB) 3180 3181 NHALF = NBAS(ISYMALP) * NRHF(ISYMK) 3182 3183 IF (LEND2 .LT. NHALF) THEN 3184 CALL QUIT('Insufficient work space in CCG_21CD (3).') 3185 END IF 3186 3187 KOFFIC = KSQC + IAODIS(ISYMALP,ISYMBET) 3188 KOFFLB = IGLMRH(ISYMBET,ISYMK) + 1 3189 KOFFLA = IGLMVI(ISYMALP,ISYMC) + 1 3190 KMO = KMOINT + IT1AM(ISYMC,ISYMK) 3191 KHALF = KEND2 3192 3193 NTOTAL = MAX(NBAS(ISYMALP),1) 3194 NTOTBE = MAX(NBAS(ISYMBET),1) 3195 NTOTC = MAX(NVIR(ISYMC),1) 3196 3197 Call DGEMM('N','N',NBAS(ISYMALP),NRHF(ISYMK),NBAS(ISYMBET), 3198 & ONE, WORK(KOFFIC),NTOTAL,XLAMDHB(KOFFLB),NTOTBE, 3199 & ZERO, WORK(KHALF),NTOTAL) 3200 3201 Call DGEMM('T','N',NVIR(ISYMC), NRHF(ISYMK), NBAS(ISYMALP), 3202 & ONE, XLAMDPA(KOFFLA),NTOTAL,WORK(KHALF),NTOTAL, 3203 & ONE, WORK(KMO),NTOTC ) 3204 END IF 3205 3206*---------------------------------------------------------------------* 3207* calculate the D contribution: 3208* FACTOR * sum_alpha I(alpha beta;j^C del) * CTR2(alpha i, j; del)^AB 3209*---------------------------------------------------------------------* 3210 ISYMAI = MULD2H(ISYMJ,ISYCTAB) 3211 ISYMI = MULD2H(ISYMALP,ISYMAI) 3212 3213 KOFFSV = KDRES + IT1AO(ISYMBET,ISYMI) 3214 KOFFZB = KCTR2AB + IT2BGD(ISYMAI,ISYMJ) 3215 & + NT1AO(ISYMAI)*(J-1) + IT1AO(ISYMALP,ISYMI) 3216 KOFFIC = KSQC + IAODIS(ISYMALP,ISYMBET) 3217 3218 NTOTBE = MAX(NBAS(ISYMBET),1) 3219 3220 NTOTAL = MAX(NBAS(ISYMALP),1) 3221 Call DGEMM('T','N',NBAS(ISYMBET),NRHF(ISYMI),NBAS(ISYMALP), 3222 & FACTOR, WORK(KOFFIC),NTOTAL,WORK(KOFFZB),NTOTAL, 3223 & ONE, WORK(KOFFSV),NTOTBE ) 3224 3225*---------------------------------------------------------------------* 3226* end loop over blocks of unpacked AO integrals 3227*---------------------------------------------------------------------* 3228 END DO ! ISYMBET 3229 3230*---------------------------------------------------------------------* 3231* calculate the C contribution: -sum_ck I(c k;j del) * CTR2(c k;a j) 3232*---------------------------------------------------------------------* 3233 ISYMAJ = MULD2H(ISYMCK,ISYCTR) 3234 ISYMA = MULD2H(ISYMJ,ISYMAJ) 3235 3236 IF (MULD2H(ISYMA,ISYDEL) .NE. ISYRES) THEN 3237 CALL QUIT('Symmetry mismatch in CCG_21CD.') 3238 END IF 3239 3240 FACT = ONE * FACTOR 3241 IF (LSAME) FACT = TWO * FACTOR 3242 3243 DO A = 1, NVIR(ISYMA) 3244 NAJ = IT1AM(ISYMA,ISYMJ) + NVIR(ISYMA)*(J-1) + A 3245 KCTR2 = IT2SQ(ISYMCK,ISYMAJ) + NT1AM(ISYMCK)*(NAJ-1) + 1 3246 KOFFR = KCRES-1 + IVIR(ISYMA)-NRHFT + A 3247 3248 WORK(KOFFR) = WORK(KOFFR) - 3249 & FACT * DDOT(NT1AM(ISYMCK),WORK(KMOINT),1,CTR2(KCTR2),1) 3250 END DO 3251 3252*---------------------------------------------------------------------* 3253* end loop over symmetry and index of J 3254*---------------------------------------------------------------------* 3255 END DO ! J 3256 END DO ! ISYMJ 3257 3258*---------------------------------------------------------------------* 3259* scale C contribution with XLAMDP0 3260*---------------------------------------------------------------------* 3261 DO ISYMI = 1, NSYM 3262 ISYMA = MULD2H(ISYMI,ISYRES) 3263 ISYMD = ISYMI 3264 3265 DO I = 1, NRHF(ISYMI) 3266 KOFFR = KCRES + IVIR(ISYMA)-NRHFT 3267 KRESV = IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I-1) + 1 3268 KLAM0 = IGLMRH(ISYMD,ISYMI) + NBAS(ISYMD)*(I-1) 3269 & + IDEL-IBAS(ISYDEL) 3270 3271 Call DAXPY(NVIR(ISYMA),XLAMDP0(KLAM0),WORK(KOFFR),1, 3272 & RESVEC(KRESV),1) 3273 3274 END DO 3275 END DO 3276 3277*---------------------------------------------------------------------* 3278* transform D contribution to MO respresentation using XLAMDH0 3279*---------------------------------------------------------------------* 3280 DO ISYMA = 1, NSYM 3281 ISYMI = MULD2H(ISYMA,ISYRES) 3282 ISYMBET = ISYMA 3283 3284 KOFFSV = KDRES + IT1AO(ISYMBET,ISYMI) 3285 KLAM0 = IGLMVI(ISYMBET,ISYMA) + 1 3286 KRESV = IT1AM(ISYMA,ISYMI) + 1 3287 3288 NTOTBE = MAX(NBAS(ISYMBET),1) 3289 NTOTA = MAX(NVIR(ISYMA),1) 3290 3291 Call DGEMM('T','N',NVIR(ISYMA),NRHF(ISYMI),NBAS(ISYMBET), 3292 & ONE, XLAMDH0(KLAM0),NTOTBE, WORK(KOFFSV),NTOTBE, 3293 & ONE, RESVEC(KRESV),NTOTA) 3294 3295 END DO 3296 3297*---------------------------------------------------------------------* 3298* that's it! return. 3299*---------------------------------------------------------------------* 3300 3301 CALL QEXIT('CCG_21CD') 3302 3303 RETURN 3304 END 3305*---------------------------------------------------------------------* 3306c/* Deck CCG_22CD */ 3307*=====================================================================* 3308 SUBROUTINE CCG_22CD(DELTA2,CDTERM,ISYRES, 3309 & XBIAJK,EBIAJK,ISYMXB, 3310 & ZBIAJK,CBIAJK,ISYMZB, 3311 & XCIAJK,ECIAJK,ISYMXC, 3312 & ZCIAJK,CCIAJK,ISYMZC, 3313 & LSAME, WORK, LWORK) 3314*---------------------------------------------------------------------* 3315* 3316* Purpose: to calculate the contribution to the G matrix which 3317* are analog to the 22C/D contribution to the left transf. 3318* 3319* symmetries & variables: 3320* 3321* ISYRES : result vector DELTA2, scratch vector CDTERM 3322* ISYMXB : one-index transf. integrals XBIAJK, EBIAJK 3323* ISYMZB : one-index transf. multipliers ZBIAJK, CBIAJK 3324* ISYMXC : one-index transf. integrals XCIAJK, ECIAJK 3325* ISYMZC : one-index transf. multipliers ZCIAJK, CCIAJK 3326* 3327* if LSAME = .TRUE. C and B perturbations are assumed to be identical 3328* 3329* Christof Haettig, September 1998 3330*=====================================================================* 3331#if defined (IMPLICIT_NONE) 3332 IMPLICIT NONE 3333#else 3334# include "implicit.h" 3335#endif 3336#include "priunit.h" 3337#include "ccsdsym.h" 3338#include "ccorb.h" 3339 3340 LOGICAL LSAME 3341 INTEGER ISYRES, ISYMXB, ISYMZB, ISYMXC, ISYMZC, LWORK 3342 3343#if defined (SYS_CRAY) 3344 REAL ZERO, HALF, ONE, TWO, THREE, FACT, DDOT 3345 REAL DELTA2(*), CDTERM(*), WORK(LWORK) 3346 REAL XBIAJK(*),EBIAJK(*),ZBIAJK(*),CBIAJK(*) 3347 REAL XCIAJK(*),ECIAJK(*),ZCIAJK(*),CCIAJK(*) 3348#else 3349 DOUBLE PRECISION DELTA2(*), CDTERM(*), WORK(LWORK) 3350 DOUBLE PRECISION XBIAJK(*),EBIAJK(*),ZBIAJK(*),CBIAJK(*) 3351 DOUBLE PRECISION XCIAJK(*),ECIAJK(*),ZCIAJK(*),CCIAJK(*) 3352 DOUBLE PRECISION ZERO, HALF, ONE, TWO, THREE, FACT, DDOT 3353#endif 3354 PARAMETER (ZERO = 0.0d0, ONE = 1.0d0, TWO = 2.0d0) 3355 PARAMETER (THREE = 3.0d0, HALF = 0.5d0) 3356 3357 INTEGER ITAIKJ(8,8), NTAIKJ(8), ISYMAIK, ISYM, ISYMJ, ICOUNT 3358 INTEGER ISYMAJ, ISYMBI, ISYMAI, ISYMBJ, ISYMKL 3359 INTEGER NTOTAJ, NTOTBI, NBJ, NJ, NAI, NAIBJ 3360 INTEGER KSCR, KOFF1, KOFF2, KOFF3 3361 3362 INTEGER INDEX 3363 INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 3364 3365 CALL QENTER('CCG_22CD') 3366 3367* check symmetries: 3368 IF (MULD2H(ISYMXB,ISYMZC).NE.ISYRES .OR. 3369 & MULD2H(ISYMXC,ISYMZB).NE.ISYRES ) THEN 3370 CALL QUIT('Symmetry mismatch in CCG_22CD.') 3371 END IF 3372 3373 DO ISYM = 1, NSYM 3374 ICOUNT = 0 3375 DO ISYMAIK = 1, NSYM 3376 ISYMJ = MULD2H(ISYMAIK,ISYM) 3377 ITAIKJ(ISYMAIK,ISYMJ) = ICOUNT 3378 ICOUNT = ICOUNT + NT2BCD(ISYMAIK)*NRHF(ISYMJ) 3379 END DO 3380 NTAIKJ(ISYM) = ICOUNT 3381 END DO 3382 3383 KSCR = 1 3384 3385 IF (LWORK .LT. MAX(NTAIKJ(ISYMXB),NTAIKJ(ISYMXC)) ) THEN 3386 CALL QUIT('Insufficient memory in CCG_22CD.') 3387 END IF 3388 3389C CALL DZERO(CDTERM,NT2SQ(ISYRES)) 3390 3391 FACT = HALF 3392 IF (LSAME) FACT = ONE 3393 3394C ---------------------------------------------------------- 3395C calculate 2 x Cou - Exc of one-index transformed integrals 3396C and resort to L(ib,kl) storage for DGEMM 3397C ---------------------------------------------------------- 3398 CALL CCG_SORT1(XBIAJK,WORK(KSCR),ISYMXB,0) 3399 CALL DSCAL(NTAIKJ(ISYMXB),-ONE,WORK(KSCR),1) 3400 CALL DAXPY(NTAIKJ(ISYMXB),TWO,XBIAJK,1,WORK(KSCR),1) 3401 CALL CCG_SORT1(WORK(KSCR),EBIAJK,ISYMXB,1) 3402 3403C ---------------------------------------------------- 3404C resort one-index transformed Lagrangian multipliers: 3405C ---------------------------------------------------- 3406 CALL CCG_SORT1(ZCIAJK,CCIAJK,ISYMZC,2) 3407 3408C --------------------------------------------- 3409C contract L^C with I^B to D term contribution: 3410C --------------------------------------------- 3411 DO ISYMBI = 1, NSYM 3412 ISYMAJ = MULD2H(ISYRES,ISYMBI) 3413 ISYMKL = MULD2H(ISYMXB,ISYMBI) 3414 3415 KOFF1 = ISAIKL(ISYMAJ,ISYMKL) + 1 3416 KOFF2 = ISAIKL(ISYMBI,ISYMKL) + 1 3417 KOFF3 = IT2SQ(ISYMAJ,ISYMBI) + 1 3418 3419 NTOTAJ = MAX(NT1AM(ISYMAJ),1) 3420 NTOTBI = MAX(NT1AM(ISYMBI),1) 3421 3422 CALL DGEMM('N','T',NT1AM(ISYMAJ),NT1AM(ISYMBI),NMATIJ(ISYMKL), 3423 & -FACT,CCIAJK(KOFF1),NTOTAJ,EBIAJK(KOFF2),NTOTBI, 3424 & ZERO,CDTERM(KOFF3),NTOTAJ) 3425 END DO 3426 3427C WRITE (LUPRI,*) 'NORM^2(EBIAJK):',DDOT(NTAIKJ(ISYMXB),EBIAJK,1,EBIAJK,1) 3428C WRITE (LUPRI,*) 'NORM^2(CCIAJK):',DDOT(NTAIKJ(ISYMZC),CCIAJK,1,CCIAJK,1) 3429 3430C WRITE (LUPRI,*) 'CDTERM after 1. contribution:' 3431C CALL CC_PRSQ(WORK,CDTERM,ISYRES,0,1) 3432 3433 IF (.NOT. LSAME) THEN 3434 3435C ---------------------------------------------------------- 3436C calculate 2 x Cou - Exc of one-index transformed integrals 3437C and resort to L(ib,kl) storage for DGEMM 3438C ---------------------------------------------------------- 3439 CALL CCG_SORT1(XCIAJK,WORK(KSCR),ISYMXC,0) 3440 CALL DSCAL(NTAIKJ(ISYMXC),-ONE,WORK(KSCR),1) 3441 CALL DAXPY(NTAIKJ(ISYMXC),TWO,XCIAJK,1,WORK(KSCR),1) 3442 CALL CCG_SORT1(WORK(KSCR),ECIAJK,ISYMXC,1) 3443 3444C ---------------------------------------------------- 3445C resort one-index transformed Lagrangian multipliers: 3446C ---------------------------------------------------- 3447 CALL CCG_SORT1(ZBIAJK,CBIAJK,ISYMZB,2) 3448 3449C WRITE (LUPRI,*) 'NORM^2(ECIAJK):',DDOT(NTAIKJ(ISYMXC),ECIAJK,1,ECIAJK,1) 3450C WRITE (LUPRI,*) 'NORM^2(CBIAJK):',DDOT(NTAIKJ(ISYMZB),CBIAJK,1,CBIAJK,1) 3451 3452C --------------------------------------------- 3453C contract L^B with I^C to D term contribution: 3454C --------------------------------------------- 3455 DO ISYMBI = 1, NSYM 3456 ISYMAJ = MULD2H(ISYRES,ISYMBI) 3457 ISYMKL = MULD2H(ISYMXC,ISYMBI) 3458 3459 KOFF1 = ISAIKL(ISYMAJ,ISYMKL) + 1 3460 KOFF2 = ISAIKL(ISYMBI,ISYMKL) + 1 3461 KOFF3 = IT2SQ(ISYMAJ,ISYMBI) + 1 3462 3463 NTOTAJ = MAX(NT1AM(ISYMAJ),1) 3464 NTOTBI = MAX(NT1AM(ISYMBI),1) 3465 3466 CALL DGEMM('N','T',NT1AM(ISYMAJ),NT1AM(ISYMBI),NMATIJ(ISYMKL), 3467 & -FACT,CBIAJK(KOFF1),NTOTAJ,ECIAJK(KOFF2),NTOTBI, 3468 & ONE,CDTERM(KOFF3),NTOTAJ) 3469 END DO 3470 3471 END IF 3472 3473C WRITE (LUPRI,*) 'CDTERM after 2. contribution:' 3474C CALL CC_PRSQ(WORK,CDTERM,ISYRES,0,1) 3475 3476C -------------------------------------- 3477C apply (2 - Pij) to the D contribution: 3478C -------------------------------------- 3479 CALL CCRHS_T2TR(CDTERM,WORK,LWORK,ISYRES) 3480 3481C ------------------------------------------------------------ 3482C transpose occupied indeces before adding the C contribution: 3483C ------------------------------------------------------------ 3484 CALL CCSD_T2TP(CDTERM,WORK,LWORK,ISYRES) 3485 3486C WRITE (LUPRI,*) 'CDTERM before 3. contribution:' 3487C CALL CC_PRSQ(WORK,CDTERM,ISYRES,0,1) 3488 3489C WRITE (LUPRI,*) 'ISYRES:',ISYRES 3490C WRITE (LUPRI,*) 'ISYMXB,ISYMXC:',ISYMXB,ISYMXC 3491C WRITE (LUPRI,*) 'ISYMZB,ISYMZC:',ISYMZB,ISYMZC 3492 3493C WRITE (LUPRI,*) 'NORM^2(XBIAJK):',DDOT(NTAIKJ(ISYMXB),XBIAJK,1,XBIAJK,1) 3494C WRITE (LUPRI,*) 'NORM^2(ZCIAJK):',DDOT(NTAIKJ(ISYMZC),ZCIAJK,1,ZCIAJK,1) 3495 3496C -------------------------------------------------- 3497C resort one-index transformed integrals to (lb,ik): 3498C -------------------------------------------------- 3499 CALL CCG_SORT1(XBIAJK,EBIAJK,ISYMXB,0) 3500C WRITE (LUPRI,*) 'NORM^2(EBIAJK):',DDOT(NTAIKJ(ISYMXB),EBIAJK,1,EBIAJK,1) 3501 CALL CCG_SORT1(EBIAJK,XBIAJK,ISYMXB,1) 3502 3503C ------------------------------------------------------------- 3504C calculate 2 x Exc + Cou of one-index transformed multipliers: 3505C and resort to Zeta(ja,kl) storage for DGEMM 3506C ------------------------------------------------------------- 3507 CALL CCG_SORT1(ZCIAJK,CCIAJK,ISYMZC,0) 3508C WRITE (LUPRI,*) 'NORM^2(CCIAJK):',DDOT(NTAIKJ(ISYMZC),CCIAJK,1,CCIAJK,1) 3509 CALL DAXPY(NTAIKJ(ISYMZC),TWO,CCIAJK,1,ZCIAJK,1) 3510C WRITE (LUPRI,*) 'NORM^2(ZCIAJK):',DDOT(NTAIKJ(ISYMZC),ZCIAJK,1,ZCIAJK,1) 3511 CALL CCG_SORT1(ZCIAJK,CCIAJK,ISYMZC,2) 3512 3513C WRITE (LUPRI,*) 'NORM^2(XBIAJK):',DDOT(NTAIKJ(ISYMXB),XBIAJK,1,XBIAJK,1) 3514C WRITE (LUPRI,*) 'NORM^2(CCIAJK):',DDOT(NTAIKJ(ISYMZC),CCIAJK,1,CCIAJK,1) 3515 3516C --------------------------------------------- 3517C contract L^C with I^B to C term contribution: 3518C --------------------------------------------- 3519 DO ISYMBI = 1, NSYM 3520 ISYMAJ = MULD2H(ISYRES,ISYMBI) 3521 ISYMKL = MULD2H(ISYMXB,ISYMBI) 3522 3523 KOFF1 = ISAIKL(ISYMAJ,ISYMKL) + 1 3524 KOFF2 = ISAIKL(ISYMBI,ISYMKL) + 1 3525 KOFF3 = IT2SQ(ISYMAJ,ISYMBI) + 1 3526 3527 NTOTAJ = MAX(NT1AM(ISYMAJ),1) 3528 NTOTBI = MAX(NT1AM(ISYMBI),1) 3529 3530 CALL DGEMM('N','T',NT1AM(ISYMAJ),NT1AM(ISYMBI),NMATIJ(ISYMKL), 3531 & FACT,CCIAJK(KOFF1),NTOTAJ,XBIAJK(KOFF2),NTOTBI, 3532 & ONE,CDTERM(KOFF3),NTOTAJ) 3533 END DO 3534 3535C WRITE (LUPRI,*) 'CDTERM after 3. contribution:' 3536C CALL CC_PRSQ(WORK,CDTERM,ISYRES,0,1) 3537 3538 IF (.NOT. LSAME) THEN 3539 3540C -------------------------------------------------- 3541C resort one-index transformed integrals to (lb,ik): 3542C -------------------------------------------------- 3543 CALL CCG_SORT1(XCIAJK,ECIAJK,ISYMXC,0) 3544 CALL CCG_SORT1(ECIAJK,XCIAJK,ISYMXC,1) 3545 3546C ------------------------------------------------------------- 3547C calculate 2 x Exc + Cou of one-index transformed multipliers: 3548C and resort to Zeta(ja,kl) storage for DGEMM 3549C ------------------------------------------------------------- 3550 CALL CCG_SORT1(ZBIAJK,CBIAJK,ISYMZB,0) 3551 CALL DAXPY(NTAIKJ(ISYMZB),TWO,CBIAJK,1,ZBIAJK,1) 3552 CALL CCG_SORT1(ZBIAJK,CBIAJK,ISYMZB,2) 3553 3554C --------------------------------------------- 3555C contract L^B with I^C to C term contribution: 3556C --------------------------------------------- 3557 DO ISYMBI = 1, NSYM 3558 ISYMAJ = MULD2H(ISYRES,ISYMBI) 3559 ISYMKL = MULD2H(ISYMXC,ISYMBI) 3560 3561 KOFF1 = ISAIKL(ISYMAJ,ISYMKL) + 1 3562 KOFF2 = ISAIKL(ISYMBI,ISYMKL) + 1 3563 KOFF3 = IT2SQ(ISYMAJ,ISYMBI) + 1 3564 3565 NTOTAJ = MAX(NT1AM(ISYMAJ),1) 3566 NTOTBI = MAX(NT1AM(ISYMBI),1) 3567 3568 CALL DGEMM('N','T',NT1AM(ISYMAJ),NT1AM(ISYMBI),NMATIJ(ISYMKL), 3569 & FACT,CBIAJK(KOFF1),NTOTAJ,XCIAJK(KOFF2),NTOTBI, 3570 & ONE,CDTERM(KOFF3),NTOTAJ) 3571 END DO 3572 3573 END IF 3574 3575C WRITE (LUPRI,*) 'CDTERM after 4. contribution:' 3576C CALL CC_PRSQ(WORK,CDTERM,ISYRES,0,1) 3577 3578C --------------------------------------------- 3579C 'back'-transposition of the occupied indeces: 3580C --------------------------------------------- 3581 CALL CCSD_T2TP(CDTERM,WORK,LWORK,ISYRES) 3582 3583*---------------------------------------------------------------------* 3584* storage in result vector: 3585*---------------------------------------------------------------------* 3586 DO ISYMBJ = 1, NSYM 3587 ISYMAI = MULD2H(ISYMBJ,ISYRES) 3588 3589 IF (ISYMAI .EQ. ISYMBJ) THEN 3590 3591 DO NBJ = 1, NT1AM(ISYMBJ) 3592 NJ = IT2SQ(ISYMAI,ISYMBJ)+NT1AM(ISYMAI)*(NBJ-1) 3593 3594 DO NAI = 1, NT1AM(ISYMAI) 3595 NAIBJ = IT2AM(ISYMAI,ISYMBJ) + INDEX(NAI,NBJ) 3596 DELTA2(NAIBJ) = DELTA2(NAIBJ) + CDTERM(NJ + NAI) 3597 END DO 3598 3599 NAIBJ = IT2AM(ISYMBJ,ISYMBJ) + INDEX(NBJ,NBJ) 3600 DELTA2(NAIBJ) = DELTA2(NAIBJ) + CDTERM(NJ + NBJ) 3601 3602 END DO 3603 3604 ELSE IF (ISYMAI .LT. ISYMBJ) THEN 3605 3606 NJ = IT2SQ(ISYMAI,ISYMBJ) + 1 3607 NAIBJ = IT2AM(ISYMAI,ISYMBJ) + 1 3608 CALL DAXPY(NT1AM(ISYMAI)*NT1AM(ISYMBJ),ONE,CDTERM(NJ),1, 3609 & DELTA2(NAIBJ),1) 3610 3611 ELSE IF (ISYMAI .GT. ISYMBJ) THEN 3612 3613 DO NBJ = 1, NT1AM(ISYMBJ) 3614 NJ = IT2SQ(ISYMAI,ISYMBJ) + NT1AM(ISYMAI)*(NBJ-1) + 1 3615 NAIBJ = IT2AM(ISYMBJ,ISYMAI) + NBJ 3616 3617 Call DAXPY(NT1AM(ISYMAI),ONE,CDTERM(NJ),1, 3618 & DELTA2(NAIBJ),NT1AM(ISYMBJ)) 3619 END DO 3620 3621 END IF 3622 3623 END DO 3624 3625 CALL QEXIT('CCG_22CD') 3626 3627 RETURN 3628 END 3629*---------------------------------------------------------------------* 3630c/* Deck CCG_4O */ 3631*=====================================================================* 3632 SUBROUTINE CCG_4O(X4O, ISYM4O, XOVOV, ISYOVOV, T1AMPB, ISYMTB, 3633 & T1AMPC, ISYMTC, WORK, LWORK, IOPT) 3634*---------------------------------------------------------------------* 3635* 3636* Purpose: transform (ia|jb) integrals with B1 and C1 amplitudes to 3637* 3638* X4O(ik,jl) = (i k^B|j l^C) + (i k^C|j l^B) 3639* = sum(ab) (ia|jb) * T1AMPB(a k) * T1AMPC(b l) 3640* +sum(ab) (ia|jb) * T1AMPC(a k) * T1AMPB(b l) 3641* 3642* IOPT = 1 : packed integrals (jb|id), output in NGAMMA format 3643* IOPT = 2 : full square of integrals, output in NGAMMA format 3644* IOPT = 3 : packed integrals (jb|id), output in N3ORHF format 3645* IOPT = 4 : full square of integrals, output in N3ORHF format 3646* 3647* needed for CCSD part of CCQR_G (22A/22B contributions) 3648* and for H matrix (CC_HMAT) 3649* 3650* Christof Haettig, August 1996, changes for H matrix February 1998 3651*=====================================================================* 3652#if defined (IMPLICIT_NONE) 3653 IMPLICIT NONE 3654#else 3655# include "implicit.h" 3656#endif 3657#include "priunit.h" 3658#include "ccsdinp.h" 3659#include "ccsdsym.h" 3660#include "ccorb.h" 3661 3662 CHARACTER MSGDBG*(16) 3663 PARAMETER (MSGDBG='[debug] CCG_4O> ') 3664 LOGICAL LOCDBG 3665 PARAMETER (LOCDBG = .FALSE.) 3666 3667 INTEGER ISYM4O, ISYOVOV, ISYMTB, ISYMTC, IOPT 3668 INTEGER LWORK 3669 3670#if defined (SYS_CRAY) 3671 REAL X4O(*) 3672 REAL XOVOV(*) 3673 REAL T1AMPB(*) ! dimension (NT1AM(ISYMTB)) 3674 REAL T1AMPC(*) ! dimension (NT1AM(ISYMTC)) 3675 REAL WORK(LWORK) 3676#else 3677 DOUBLE PRECISION X4O(*) 3678 DOUBLE PRECISION XOVOV(*) 3679 DOUBLE PRECISION T1AMPB(*) ! dimension (NT1AM(ISYMTB)) 3680 DOUBLE PRECISION T1AMPC(*) ! dimension (NT1AM(ISYMTC)) 3681 DOUBLE PRECISION WORK(LWORK) 3682#endif 3683 3684 INTEGER KEND0, ISYMB, ISYMIAJ, ISYCIKJ, ISYBIKJ, KXCIKJ, KXBIKJ 3685 INTEGER KEND1, LEND1, ISYTCTB, IOPT2 3686 3687C INTEGER INDEX 3688C INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 3689 3690 CALL QENTER('CCG_4O') 3691 3692* begin: 3693 IF (LOCDBG) THEN 3694 WRITE (LUPRI,*) MSGDBG, 'just entered this routine...' 3695 END IF 3696 3697* check symmetries: 3698 ISYTCTB = MULD2H(ISYMTC,ISYMTB) 3699 IF (ISYM4O .NE. MULD2H(ISYTCTB,ISYOVOV)) THEN 3700 CALL QUIT('Symmetry mismatch in CCG_4O.') 3701 END IF 3702 3703* check IOPT and initialize result array: 3704 If (IOPT.EQ.1 .OR. IOPT.EQ.2) THEN 3705 CALL DZERO( X4O, NGAMMA(ISYM4O) ) 3706 IOPT2 = IOPT 3707 ELSE IF (IOPT.EQ.3 .OR. IOPT.EQ.4) THEN 3708 CALL DZERO( X4O, N3ORHF(ISYM4O) ) 3709 IOPT2 = IOPT - 2 3710 ELSE 3711 CALL QUIT('Illegal value for IOPT in CCG_4O.') 3712 ENDIF 3713 3714* set start of work space: 3715 KEND0 = 1 3716 3717 3718* start outer loop over symmetry blocks: 3719 DO ISYMB = 1, NSYM 3720 3721*---------------------------------------------------------------------* 3722* memory allocation: 3723*---------------------------------------------------------------------* 3724 ISYMIAJ = MULD2H(ISYMB,ISYOVOV) 3725 ISYCIKJ = MULD2H(ISYMIAJ,ISYMTC) 3726 ISYBIKJ = MULD2H(ISYMIAJ,ISYMTB) 3727 3728 KXCIKJ = KEND0 3729 KXBIKJ = KEND0 3730 KEND1 = KEND0 + MAX( NMAIJK(ISYCIKJ) , NMAIJK(ISYBIKJ) ) 3731 LEND1 = LWORK - KEND1 3732 3733 IF (LEND1 .LT. 0) THEN 3734 CALL QUIT('Insufficient work space in CCG_4O.') 3735 END IF 3736 3737 DO B = 1, NVIR(ISYMB) 3738 3739*---------------------------------------------------------------------* 3740* calculate one-index transformed (i k^C|j b) distributions 3741* and update X4O(ik,jl) by (i k^C|j b) * B1(b l) for all i,k,j,l 3742*---------------------------------------------------------------------* 3743 Call CCG_TRANS2(WORK(KXCIKJ), ISYCIKJ, XOVOV, ISYOVOV, 3744 & T1AMPC, ISYMTC, B, ISYMB, IOPT2) 3745 3746 Call CCG_4OS(X4O, ISYM4O, WORK(KXCIKJ), ISYCIKJ, 3747 & T1AMPB, ISYMTB, B, ISYMB, IOPT) 3748*---------------------------------------------------------------------* 3749* calculate one-index transformed (i k^B|j b) distributions 3750* and update X4O(ik,jl) by (i k^B|j b) * C1(b l) for all i,k,j,l 3751*---------------------------------------------------------------------* 3752 Call CCG_TRANS2(WORK(KXBIKJ), ISYBIKJ, XOVOV, ISYOVOV, 3753 & T1AMPB, ISYMTB, B, ISYMB, IOPT2) 3754 3755 Call CCG_4OS(X4O, ISYM4O, WORK(KXBIKJ), ISYBIKJ, 3756 & T1AMPC, ISYMTC, B, ISYMB, IOPT) 3757 3758*---------------------------------------------------------------------* 3759 END DO ! B 3760 END DO ! ISYMB 3761 3762 IF (LOCDBG) THEN 3763 WRITE (LUPRI,*) MSGDBG, 'leaving this routine...' 3764 END IF 3765 3766 CALL QEXIT('CCG_4O') 3767 3768 RETURN 3769 END 3770*---------------------------------------------------------------------* 3771c/* Deck CCG_4OS */ 3772*=====================================================================* 3773 SUBROUTINE CCG_4OS(X4O,ISYM4O, XIKJB,ISYIKJ, 3774 & T1,ISYMT1, B,ISYMB, IOPT) 3775*---------------------------------------------------------------------* 3776* 3777* Purpose: update I(ik,jl) integrals by: 3778* 3779* X4O(ik,jl) = X4O(ik,jl) + (i k|j b) * T1(b l) for all i,k,j,l 3780* 3781* called by CCG_4O 3782* 3783* Christof Haettig, August 1996, changes for H matrix: Februar 1998 3784*=====================================================================* 3785#if defined (IMPLICIT_NONE) 3786 IMPLICIT NONE 3787#else 3788# include "implicit.h" 3789#endif 3790#include "priunit.h" 3791#include "ccsdinp.h" 3792#include "ccsdsym.h" 3793#include "ccorb.h" 3794 3795 CHARACTER MSGDBG*(17) 3796 PARAMETER (MSGDBG='[debug] CCG_4OS> ') 3797 LOGICAL LOCDBG 3798 PARAMETER (LOCDBG = .FALSE.) 3799 3800 INTEGER ISYM4O, ISYIKJ, ISYMT1, ISYMB, IOPT 3801 3802#if defined (SYS_CRAY) 3803 REAL X4O(*) 3804 REAL XIKJB(*) ! dimension (NMAIJK(ISYIKJ)) 3805 REAL T1(*) ! dimension (NT1AM(ISYMT1)) 3806#else 3807 DOUBLE PRECISION X4O(*) 3808 DOUBLE PRECISION XIKJB(*) ! dimension (NMAIJK(ISYIKJ)) 3809 DOUBLE PRECISION T1(*) ! dimension (NT1AM(ISYMT1)) 3810#endif 3811 3812 INTEGER ISYML, NBL, ISYMJ, ISYMJL, ISYMIK, NJL, NIK, NIKJ, NIKJL 3813 INTEGER NLJ, NLJKI, OFFIK, OFFI, ISYLJK, ISYMI, ISYMLJ, ISYMK 3814 3815 INTEGER INDEX 3816 INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 3817 3818 CALL QENTER('CCG_4OS') 3819 3820* begin: 3821 IF (LOCDBG) THEN 3822 WRITE (LUPRI,*) MSGDBG, 'just entered this routine...' 3823 END IF 3824 3825* check symmetries: 3826 IF (MULD2H(ISYM4O,ISYMT1) .NE. MULD2H(ISYIKJ,ISYMB)) THEN 3827 CALL QUIT('Symmetry mismatch in CCG_4OS.') 3828 END IF 3829 3830*=====================================================================* 3831* store result in NGAMMA format: 3832*=====================================================================* 3833 IF (IOPT.EQ.1 .OR. IOPT.EQ.2) THEN 3834 3835*---------------------------------------------------------------------* 3836* start loop over occupied index l: 3837*---------------------------------------------------------------------* 3838 ISYML = MULD2H(ISYMB,ISYMT1) 3839 3840 DO L = 1, NRHF(ISYML) 3841 3842 NBL = IT1AM(ISYMB,ISYML) + NVIR(ISYMB)*(L-1) + B 3843 3844*---------------------------------------------------------------------* 3845* add contribution to X4O: 3846*---------------------------------------------------------------------* 3847 DO ISYMJ = 1, NSYM 3848 3849 ISYMJL = MULD2H(ISYMJ,ISYML) 3850 ISYMIK = MULD2H(ISYMJL,ISYM4O) 3851 3852 IF (ISYMIK .LT. ISYMJL) THEN 3853 DO J = 1, NRHF(ISYMJ) 3854 NJL = IMATIJ(ISYMJ,ISYML) + NRHF(ISYMJ)*(L-1) + J 3855 3856 DO NIK = 1, NMATIJ(ISYMIK) 3857 NIKJ = IMAIJK(ISYMIK,ISYMJ) 3858 & + NMATIJ(ISYMIK)*(J-1) + NIK 3859 NIKJL = IGAMMA(ISYMIK,ISYMJL) 3860 & + NMATIJ(ISYMIK)*(NJL-1) + NIK 3861 3862 X4O(NIKJL) = X4O(NIKJL) + XIKJB(NIKJ) * T1(NBL) 3863 END DO 3864 END DO 3865 3866 ELSE IF (ISYMIK .EQ. ISYMJL) THEN 3867 3868 DO J = 1, NRHF(ISYMJ) 3869 NJL = IMATIJ(ISYMJ,ISYML) + NRHF(ISYMJ)*(L-1) + J 3870 3871 DO NIK = 1, NJL 3872 NIKJ = IMAIJK(ISYMIK,ISYMJ) 3873 & + NMATIJ(ISYMIK)*(J-1) + NIK 3874 NIKJL = IGAMMA(ISYMIK,ISYMJL) + INDEX(NIK,NJL) 3875 3876 X4O(NIKJL) = X4O(NIKJL) + XIKJB(NIKJ) * T1(NBL) 3877 END DO 3878 END DO 3879 3880 END IF 3881 3882 END DO ! ISYMJ 3883 3884 END DO ! L 3885 3886*=====================================================================* 3887* store result in N3ORHF format: 3888*=====================================================================* 3889 ELSE IF (IOPT.EQ.3 .OR. IOPT.EQ.4) THEN 3890 3891 ISYML = MULD2H(ISYMB,ISYMT1) 3892 3893 DO ISYMJ = 1, NSYM 3894 DO ISYMI = 1, NSYM 3895 3896 ISYMLJ = MULD2H(ISYML,ISYMJ) 3897 ISYMIK = MULD2H(ISYMLJ,ISYM4O) 3898 ISYMK = MULD2H(ISYMIK,ISYMI) 3899 ISYLJK = MULD2H(ISYMLJ,ISYMK) 3900 3901 DO I = 1, NRHF(ISYMI) 3902 OFFI = I3ORHF(ISYLJK,ISYMI) + NMAIJK(ISYLJK)*(I-1) 3903 3904 DO K = 1, NRHF(ISYMK) 3905 3906 OFFIK = OFFI + IMAIJK(ISYMLJ,ISYMK) + NMATIJ(ISYMLJ)*(K-1) 3907 NIK = IMATIJ(ISYMI,ISYMK) + NRHF(ISYMI)*(K-1) + I 3908 3909 DO J = 1, NRHF(ISYMJ) 3910 3911 NIKJ = IMAIJK(ISYMIK,ISYMJ) + NMATIJ(ISYMIK)*(J-1) + NIK 3912 3913 DO L = 1, NRHF(ISYML) 3914 3915 NBL = IT1AM(ISYMB,ISYML) + NVIR(ISYMB)*(L-1) + B 3916 NLJ = IMATIJ(ISYML,ISYMJ) + NRHF(ISYML)*(J-1) + L 3917 NLJKI = OFFIK + NLJ 3918 3919 X4O(NLJKI) = X4O(NLJKI) + XIKJB(NIKJ) * T1(NBL) 3920 3921 END DO ! L 3922 3923 END DO ! J 3924 3925 END DO ! K 3926 3927 END DO ! I 3928 3929 END DO ! ISYMI 3930 END DO ! ISYMJ 3931 3932*=====================================================================* 3933* return: 3934*=====================================================================* 3935 ELSE 3936 CALL QUIT('Illegal value for IOPT in CCG_4OS.') 3937 END IF 3938 3939 IF (LOCDBG) THEN 3940 WRITE (LUPRI,*) MSGDBG, 'leaving this routine...' 3941 END IF 3942 3943 CALL QEXIT('CCG_4OS') 3944 3945 RETURN 3946 END 3947*---------------------------------------------------------------------* 3948c/* Deck CCG_OOVV */ 3949*=====================================================================* 3950 SUBROUTINE CCG_OOVV(Jint, Kint, ISYJKINT, XOVOV, ISYOVOV, 3951 & B1AMP, ISYMTB, C1AMP, ISYMTC, 3952 & WORK, LWORK, D, ISYMD, IOPT,IPCK) 3953*---------------------------------------------------------------------* 3954* 3955* Purpose: transform a three-index batch of (jb|id) integrals 3956* (with fixed index d) with B1 and C1 amplitudes to 3957* 3958* Jint(el,j) = (e^C l^B|j d) + (e^B l^C|j d) 3959* and 3960* Kint(el,j) = (j l^B|e^C d) + (j l^C|e^B d) 3961* 3962* 3963* IOPT = 1 : Jint only 3964* IOPT = 2 : Kint only 3965* IOPT = 3 : Jint & Kint 3966* 3967* IPCK = 1 : packed integrals (jb|id) 3968* IPCK = 2 : full square of integrals 3969* 3970* needed for CCSD part of CCQR_G (CCG_22C contribiution) 3971* 3972* Christof Haettig, August 1996 3973*=====================================================================* 3974#if defined (IMPLICIT_NONE) 3975 IMPLICIT NONE 3976#else 3977# include "implicit.h" 3978#endif 3979#include "priunit.h" 3980#include "ccsdsym.h" 3981#include "ccorb.h" 3982 3983#if defined (SYS_CRAY) 3984 REAL ZERO, ONE 3985#else 3986 DOUBLE PRECISION ZERO, ONE 3987#endif 3988 PARAMETER (ZERO = 0.0d0, ONE = 1.0d0) 3989 3990 INTEGER ISYJKINT, ISYOVOV, ISYMTB, ISYMTC, ISYMD 3991 INTEGER IPCK, LWORK, IOPT 3992 3993#if defined (SYS_CRAY) 3994 REAL Kint(*) ! dimension (NT2BCD(ISYJKINT)) 3995 REAL Jint(*) ! dimension (NT2BCD(ISYJKINT)) 3996 REAL XOVOV(*) 3997 REAL B1AMP(*) ! dimension (NT1AM(ISYMTB)) 3998 REAL C1AMP(*) ! dimension (NT1AM(ISYMTC)) 3999 REAL WORK(LWORK) 4000#else 4001 DOUBLE PRECISION Kint(*) ! dimension (NT2BCD(ISYJKINT)) 4002 DOUBLE PRECISION Jint(*) ! dimension (NT2BCD(ISYJKINT)) 4003 DOUBLE PRECISION XOVOV(*) 4004 DOUBLE PRECISION B1AMP(*) ! dimension (NT1AM(ISYMTB)) 4005 DOUBLE PRECISION C1AMP(*) ! dimension (NT1AM(ISYMTC)) 4006 DOUBLE PRECISION WORK(LWORK) 4007#endif 4008 4009 INTEGER ISYTCTB, ISYMIAJ, ISYCJLI, ISYBJLI, ISYJLE, ISYMEL 4010 INTEGER KXCJLI, KXBJLI, JXCJLI, JXBJLI, KEND1, LEND1, ISYML 4011 INTEGER ISYMJL, ISYME, KKSCR, KJSCR, KEND2, LEND2, ISYMI, NELJ 4012 INTEGER KOFFK, KOFFJ, KT1, NTOTJL, NTOTE, ISYMJ, NJL, NJLE, NEL 4013 4014C INTEGER INDEX 4015C INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 4016 4017 CALL QENTER('CCG_OOVV') 4018 4019* check symmetries: 4020 ISYTCTB = MULD2H(ISYMTC,ISYMTB) 4021 4022 IF (MULD2H(ISYOVOV,ISYTCTB) .NE. MULD2H(ISYJKINT,ISYMD)) THEN 4023 CALL QUIT('Symmetry mismatch in CCG_OOVV.') 4024 END IF 4025 4026*---------------------------------------------------------------------* 4027* memory allocation 4028*---------------------------------------------------------------------* 4029 ISYMIAJ = MULD2H(ISYMD,ISYOVOV) 4030 ISYCJLI = MULD2H(ISYMIAJ,ISYMTC) 4031 ISYBJLI = MULD2H(ISYMIAJ,ISYMTB) 4032 ISYJLE = MULD2H(ISYMIAJ,ISYTCTB) 4033 4034 KXCjli = 1 4035 KXBjli = KXCjli + NMAIJK(ISYCJLI) 4036 JXCjli = KXBjli + NMAIJK(ISYBJLI) 4037 JXBjli = JXCjli + NMAIJK(ISYCJLI) 4038 KEND1 = JXBjli + NMAIJK(ISYBJLI) 4039 LEND1 = LWORK - KEND1 4040 4041 IF (LEND1 .LT. 0) THEN 4042 CALL QUIT('Insufficient work space in CCG_OOVV.') 4043 END IF 4044 4045*---------------------------------------------------------------------* 4046* calculate one-index transformed (j l^C|i d) distribution 4047* and contruct distribution with j and i indeces exchanged 4048*---------------------------------------------------------------------* 4049 Call CCG_TRANS2(WORK(KXCjli),ISYCJLI,XOVOV,ISYOVOV,C1AMP,ISYMTC, 4050 & D,ISYMD,IPCK) 4051 4052 IF (IOPT .EQ. 1 .OR. IOPT .EQ. 3) THEN 4053 Call CCLT_PI3O(WORK(JXCjli),WORK(KXCjli),ISYCJLI) 4054 END IF 4055 4056*---------------------------------------------------------------------* 4057* calculate one-index transformed (j l^B|i d) distribution 4058* and contruct distribution with j and i indeces exchanged 4059*---------------------------------------------------------------------* 4060 Call CCG_TRANS2(WORK(KXBjli),ISYBJLI,XOVOV,ISYOVOV,B1AMP,ISYMTB, 4061 & D,ISYMD,IPCK) 4062 4063 IF (IOPT .EQ. 1 .OR. IOPT .EQ. 3) THEN 4064 Call CCLT_PI3O(WORK(JXBjli),WORK(KXBjli),ISYBJLI) 4065 END IF 4066 4067*---------------------------------------------------------------------* 4068* calculate I^d(el,j) = (j l^C|e^B d) + (j l^B|e^C d): 4069*---------------------------------------------------------------------* 4070 DO ISYMJL = 1, NSYM 4071 ISYME = MULD2H(ISYMJL,ISYJLE) 4072 4073 KKSCR = KEND1 4074 KJSCR = KKSCR + NMATIJ(ISYMJL) * NVIR(ISYME) 4075 KEND2 = KJSCR + NMATIJ(ISYMJL) * NVIR(ISYME) 4076 LEND2 = LWORK - KEND2 4077 4078 IF (LEND2 .LT. 0) THEN 4079 CALL QUIT('Insufficient work space in CCG_22C.') 4080 END IF 4081 4082*.....................................................................* 4083* first contribution Kint: -sum_i (j l^C|i d) * T1B(e i) 4084*.....................................................................* 4085 ISYMI = MULD2H(ISYMJL,ISYCJLI) 4086 4087 KOFFK = KXCjli + IMAIJK(ISYMJL,ISYMI) 4088 KOFFJ = JXCjli + IMAIJK(ISYMJL,ISYMI) 4089 KT1 = IT1AM(ISYME,ISYMI) + 1 4090 4091 NTOTJL = MAX(NMATIJ(ISYMJL),1) 4092 NTOTE = MAX(NVIR(ISYME),1) 4093 4094 IF (IOPT .EQ. 2 .OR. IOPT .EQ. 3) THEN 4095 Call DGEMM('N','T',NMATIJ(ISYMJL),NVIR(ISYME),NRHF(ISYMI), 4096 & -ONE, WORK(KOFFK), NTOTJL, B1AMP(KT1), NTOTE, 4097 & ZERO, WORK(KKSCR), NTOTJL ) 4098 END IF 4099 4100 IF (IOPT .EQ. 1 .OR. IOPT .EQ. 3) THEN 4101*.....................................................................* 4102* analogously for Jint: -sum_i (i l^C|j d) * T1B(e i) 4103*.....................................................................* 4104 Call DGEMM('N','T',NMATIJ(ISYMJL),NVIR(ISYME),NRHF(ISYMI), 4105 & -ONE, WORK(KOFFJ), NTOTJL, B1AMP(KT1), NTOTE, 4106 & ZERO, WORK(KJSCR), NTOTJL ) 4107 END IF 4108 4109 4110*.....................................................................* 4111* add second contribution: -sum_i (j l^C|i d) * T1C(e i) 4112*.....................................................................* 4113 ISYMI = MULD2H(ISYMJL,ISYBJLI) 4114 4115 KOFFK = KXBjli + IMAIJK(ISYMJL,ISYMI) 4116 KOFFJ = JXBjli + IMAIJK(ISYMJL,ISYMI) 4117 KT1 = IT1AM(ISYME,ISYMI) + 1 4118 4119 NTOTJL = MAX(NMATIJ(ISYMJL),1) 4120 NTOTE = MAX(NVIR(ISYME),1) 4121 4122 IF (IOPT .EQ. 2 .OR. IOPT .EQ. 3) THEN 4123 Call DGEMM('N','T',NMATIJ(ISYMJL),NVIR(ISYME),NRHF(ISYMI), 4124 & -ONE, WORK(KOFFK), NTOTJL, C1AMP(KT1), NTOTE, 4125 & +ONE, WORK(KKSCR), NTOTJL ) 4126 END IF 4127 4128 IF (IOPT .EQ. 1 .OR. IOPT .EQ. 3) THEN 4129*.....................................................................* 4130* analogously for Jint: -sum_i (i l^C|j d) * T1C(e i) 4131*.....................................................................* 4132 Call DGEMM('N','T',NMATIJ(ISYMJL),NVIR(ISYME),NRHF(ISYMI), 4133 & -ONE, WORK(KOFFJ), NTOTJL, C1AMP(KT1), NTOTE, 4134 & +ONE, WORK(KJSCR), NTOTJL ) 4135 END IF 4136 4137 4138*---------------------------------------------------------------------* 4139* resort result from I(j l;e d) to I(e l;j d): 4140*---------------------------------------------------------------------* 4141 DO ISYMJ = 1, NSYM 4142 ISYML = MULD2H(ISYMJ,ISYMJL) 4143 ISYMEL = MULD2H(ISYME,ISYML) 4144 4145 IF (IOPT .EQ. 2 .OR. IOPT .EQ. 3) THEN 4146 DO J = 1, NRHF(ISYMJ) 4147 DO L = 1, NRHF(ISYML) 4148 DO E = 1, NVIR(ISYME) 4149 NJL = IMATIJ(ISYMJ,ISYML) + NRHF(ISYMJ)*(L-1) + J 4150 NJLE = NMATIJ(ISYMJL)*(E-1) + NJL 4151 4152 NEL = IT1AM(ISYME,ISYML) + NVIR(ISYME)*(L-1) + E 4153 NELJ = IT2BCD(ISYMEL,ISYMJ) + NT1AM(ISYMEL)*(J-1) + NEL 4154 4155 Kint(NELJ) = WORK(KKSCR-1 + NJLE) 4156 4157 END DO 4158 END DO 4159 END DO 4160 END IF 4161 4162 IF (IOPT .EQ. 1 .OR. IOPT .EQ. 3) THEN 4163 DO J = 1, NRHF(ISYMJ) 4164 DO L = 1, NRHF(ISYML) 4165 DO E = 1, NVIR(ISYME) 4166 NJL = IMATIJ(ISYMJ,ISYML) + NRHF(ISYMJ)*(L-1) + J 4167 NJLE = NMATIJ(ISYMJL)*(E-1) + NJL 4168 4169 NEL = IT1AM(ISYME,ISYML) + NVIR(ISYME)*(L-1) + E 4170 NELJ = IT2BCD(ISYMEL,ISYMJ) + NT1AM(ISYMEL)*(J-1) + NEL 4171 4172 Jint(NELJ) = WORK(KJSCR-1 + NJLE) 4173 4174 END DO 4175 END DO 4176 END DO 4177 END IF 4178 4179 END DO ! ISYMJ 4180 4181 END DO ! ISYMJL 4182 4183 CALL QEXIT('CCG_OOVV') 4184 4185 RETURN 4186 END 4187*---------------------------------------------------------------------* 4188c/* Deck CCG_RDIAJB */ 4189*=====================================================================* 4190 SUBROUTINE CCG_RDIAJB(XIAJB,LIAJB) 4191*---------------------------------------------------------------------* 4192* Purpose: read packed (ia|jb) integrals into memory 4193* integrals read from unit 52, which is assumed open 4194* 4195* Christof Haettig, August 1996 4196*=====================================================================* 4197#if defined (IMPLICIT_NONE) 4198 IMPLICIT NONE 4199#else 4200# include "implicit.h" 4201#endif 4202#include "priunit.h" 4203#include "ccsdinp.h" 4204#include "ccsdsym.h" 4205#include "iratdef.h" 4206#include "ccinftap.h" 4207 4208 LOGICAL LOCDBG 4209 PARAMETER (LOCDBG = .FALSE.) 4210 4211 INTEGER LIAJB 4212 4213#if defined (SYS_CRAY) 4214 REAL XIAJB(LIAJB) 4215#else 4216 DOUBLE PRECISION XIAJB(LIAJB) 4217#endif 4218 4219 CALL QENTER('CCG_RDIAJB') 4220 4221* check dimension: 4222 IF (LIAJB .LT. NT2AM(ISYMOP)) THEN 4223 CALL QUIT( 4224 * 'Insufficient memory for (ia|jb) integrals in CCG_RDIAJB.') 4225 END IF 4226 4227* rewind integral file: 4228 REWIND(LUIAJB) 4229 4230* read (ia|jb) integrals: 4231 CALL READI(LUIAJB,IRAT*NT2AM(1),XIAJB) 4232 4233 IF (LOCDBG) THEN 4234 CALL AROUND( 'CCG_RDIAJB: Integrals (ia|jb) ' ) 4235 CALL CC_PRP(0.00D0,XIAJB,ISYMOP,0,1) 4236 ENDIF 4237 4238 CALL QEXIT('CCG_RDIAJB') 4239 4240 RETURN 4241 END 4242*---------------------------------------------------------------------* 4243c/* Deck CCG_TI */ 4244*=====================================================================* 4245 SUBROUTINE CCG_TI(Tint,ISYMTI,T2amp,ISYMTAM,C,ISYMC) 4246*---------------------------------------------------------------------* 4247* 4248* Purpose: to get a three-index submatrix T^{c}(bk,j) out of the 4249* packed matrix of T(dk,cj) amplitudes 4250* 4251* T2amp -- packed matrix of T2 amplitudes 4252* Tint -- three-index submatrix of T2 with fixed c 4253* 4254* needed for CCSD part of CCQR_G 4255* 4256* Christof Haettig, August 1996 4257*=====================================================================* 4258#if defined (IMPLICIT_NONE) 4259 IMPLICIT NONE 4260#else 4261# include "implicit.h" 4262#endif 4263#include "priunit.h" 4264#include "ccorb.h" 4265#include "ccsdsym.h" 4266 4267#if defined (SYS_CRAY) 4268 REAL T2amp(*) ! dimension (NT2AM(ISYMTAM)) 4269 REAL Tint(*) ! dimension (NT2BCD(ISYMTI)) 4270#else 4271 DOUBLE PRECISION T2amp(*) ! dimension (NT2AM(ISYMTAM)) 4272 DOUBLE PRECISION Tint(*) ! dimension (NT2BCD(ISYMTI)) 4273#endif 4274 4275 INTEGER ISYMTI, ISYMTAM, ISYMC, ISYMDK, ISYMCJ, ISYMJ 4276 INTEGER NCJ, NDK, NDKCJ, NDKJ 4277 4278 INTEGER INDEX 4279 INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 4280 4281 CALL QENTER('CCG_TI') 4282 4283* check symmetries: 4284 IF (ISYMTI .NE. MULD2H(ISYMTAM,ISYMC)) THEN 4285 CALL QUIT('SYMMETRY MISMATCH IN CCG_TI.') 4286 END IF 4287 4288 4289 DO ISYMDK = 1, NSYM 4290 ISYMCJ = MULD2H(ISYMDK,ISYMTAM) 4291 ISYMJ = MULD2H(ISYMCJ,ISYMC) 4292 4293 DO J = 1, NRHF(ISYMJ) 4294 4295 NCJ = IT1AM(ISYMC,ISYMJ) + NVIR(ISYMC)*(J-1) + C 4296 4297 IF (ISYMDK. EQ. ISYMCJ) THEN 4298 4299 DO NDK = 1, NT1AM(ISYMDK) 4300 NDKCJ = IT2AM(ISYMDK,ISYMCJ) + INDEX(NDK,NCJ) 4301 NDKJ = IT2BCD(ISYMDK,ISYMJ) + NT1AM(ISYMDK)*(J-1) + NDK 4302 Tint(NDKJ) = T2amp(NDKCJ) 4303 END DO 4304 4305 ELSE IF (ISYMDK. LT. ISYMCJ) THEN 4306 4307 NDKCJ = IT2AM(ISYMDK,ISYMCJ) + NT1AM(ISYMDK)*(NCJ-1) + 1 4308 NDKJ = IT2BCD(ISYMDK,ISYMJ) + NT1AM(ISYMDK)*(J-1) + 1 4309 Call DCOPY(NT1AM(ISYMDK),T2amp(NDKCJ),1,Tint(NDKJ),1) 4310 4311 ELSE IF (ISYMDK. GT. ISYMCJ) THEN 4312 4313 NDKCJ = IT2AM(ISYMCJ,ISYMDK) + NCJ 4314 NDKJ = IT2BCD(ISYMDK,ISYMJ) + NT1AM(ISYMDK)*(J-1) + 1 4315 Call DCOPY(NT1AM(ISYMDK),T2amp(NDKCJ),NT1AM(ISYMCJ), 4316 & Tint(NDKJ),1) 4317 4318 END IF 4319 4320 END DO ! J 4321 END DO ! ISYMDK 4322 4323 CALL QEXIT('CCG_TI') 4324C 4325 RETURN 4326 END 4327*---------------------------------------------------------------------* 4328c/* Deck CCG_TRANS13 */ 4329*=====================================================================* 4330 SUBROUTINE CCG_TRANS13(Xiaec,ISYMIAE,XOVOV,ISYOVOV, 4331 & B1AMP,ISYMB1,C,ISYMC,IPOS,WORK,LWORK) 4332*---------------------------------------------------------------------* 4333* 4334* Purpose: transform the first (IPOS=1) or the third (IPOS=3) index 4335* of a three-index batch of (i a|k c) integrals (with 4336* fixed index c) with the single excitation amplitudes 4337* stored in B1AMP. 4338* 4339* NOTE the minus sign in the transformation! 4340* 4341* 4342* IPOS = 1 : (e^B a|i c) = sum_k -B1(ek) * (k a|i c) 4343* 4344* IPOS = 3 : (i a|e^B c) = sum_k -B1(ek) * (i a|k c) 4345* 4346* the result is in both cases returned in Xiaec stored as 4347* (a i|e^B;c) using the dimension/adress arrays NCKATR & ICKATR 4348* 4349* (needed for CCSD part of CCQR_G) 4350* 4351* Christof Haettig, August 1996 4352*=====================================================================* 4353#if defined (IMPLICIT_NONE) 4354 IMPLICIT NONE 4355#else 4356# include "implicit.h" 4357#endif 4358#include "priunit.h" 4359#include "ccsdinp.h" 4360#include "ccsdsym.h" 4361#include "ccorb.h" 4362 4363 CHARACTER MSGDBG*(21) 4364 PARAMETER (MSGDBG='[debug] CCG_TRANS13> ') 4365 LOGICAL LOCDBG 4366 PARAMETER (LOCDBG = .FALSE.) 4367 4368 INTEGER ISYMIAE, ISYOVOV, ISYMB1, ISYMC, LWORK, IPOS 4369 4370#if defined (SYS_CRAY) 4371 REAL XIAEC(*) ! dimension (NCKATR(ISYMIAE)) 4372 REAL XOVOV(*) ! dimension (NT2AM(ISYOVOV)) 4373 REAL B1AMP(*) ! dimension (NT1AM(ISYMB1)) 4374 REAL WORK(LWORK) 4375#else 4376 DOUBLE PRECISION XIAEC(*) ! dimension (NCKATR(ISYMIAE)) 4377 DOUBLE PRECISION XOVOV(*) ! dimension (NT2AM(ISYOVOV)) 4378 DOUBLE PRECISION B1AMP(*) ! dimension (NT1AM(ISYMB1)) 4379 DOUBLE PRECISION WORK(LWORK) 4380#endif 4381 4382 INTEGER ISYMAIK, KAIK, KFREE, LFREE, ISYMAI, ISYMK, ISYME 4383 INTEGER NEK, NAI, NAIE, NAIK 4384 4385C INTEGER INDEX 4386C INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 4387 4388 CALL QENTER('CCG_TRANS13') 4389 4390* check symmetries and IPOS option: 4391 IF (MULD2H(ISYMIAE,ISYOVOV) .NE. MULD2H(ISYMB1,ISYMC)) THEN 4392 CALL QUIT('Symmetry mismatch in CCG_TRANS13.') 4393 END IF 4394 IF (IPOS .NE. 1 .AND. IPOS .NE. 3) THEN 4395 CALL QUIT('CCG_TRANS13 called with illegal value for IPOS.') 4396 END IF 4397 4398* initialize result: 4399 Call DZERO (Xiaec, NCKATR(ISYMIAE)) 4400 4401* get submatrix I(ai,k) of integrals: 4402 ISYMAIK = MULD2H(ISYMC,ISYOVOV) 4403 KAIK = 1 4404 KFREE = KAIK + NT2BCD(ISYMAIK) 4405 LFREE = LWORK - NT2BCD(ISYMAIK) 4406 4407 IF (LFREE .lt. 0) THEN 4408 CALL QUIT('Insufficient work space in CCG_TRANS13.') 4409 END IF 4410 4411 Call CCG_TI (WORK(KAIK),ISYMAIK,XOVOV,ISYOVOV,C,ISYMC) 4412 4413* if the first index should be transformed, exchange i and k: 4414 IF (IPOS .eq. 1) THEN 4415 Call CCLT_P21I(WORK(KAIK),ISYMAIK,WORK(KFREE),LFREE, 4416 & IT2BCD, NT2BCD, IT1AM, NT1AM, NVIR) 4417 END IF 4418 4419* do the transformation: 4420 DO ISYMAI = 1, NSYM 4421 ISYMK = MULD2H(ISYMAI,ISYMAIK) 4422 ISYME = MULD2H(ISYMK,ISYMB1) 4423 4424 DO K = 1, NRHF(ISYMK) 4425 DO E = 1, NVIR(ISYME) 4426 NEK = IT1AM(ISYME,ISYMK) + NVIR(ISYME)*(K-1) + E 4427 DO NAI = 1, NT1AM(ISYMAI) 4428 NAIE = ICKATR(ISYMAI,ISYME) + NT1AM(ISYMAI)*(E-1) + NAI 4429 NAIK = IT2BCD(ISYMAI,ISYMK) + NT1AM(ISYMAI)*(K-1) + NAI 4430 Xiaec(NAIE) = Xiaec(NAIE) - B1AMP(NEK) * WORK(KAIK-1+NAIK) 4431 END DO 4432 END DO 4433 END DO 4434 4435 IF (LOCDBG) THEN 4436 DO E = 1, NVIR(ISYME) 4437 WRITE (LUPRI,*) MSGDBG, 'ISYMAI,E:',ISYMAI,E 4438 NAIE = ICKATR(ISYMAI,ISYME) + NT1AM(ISYMAI)*(E-1) + 1 4439 Call CC_PRP(Xiaec(NAIE),WORK,ISYMAI,1,0) 4440 END DO 4441 END IF 4442 4443 END DO 4444 4445 4446 CALL QEXIT('CCG_TRANS13') 4447C 4448 RETURN 4449 END 4450*---------------------------------------------------------------------* 4451c/* Deck CCG_TRANS2 */ 4452*=====================================================================* 4453 SUBROUTINE CCG_TRANS2(Xjlic,ISYMJLI,XOVOV,ISYOVOV,B1AMP,ISYMB1, 4454 & C, ISYMC, IOPT) 4455*---------------------------------------------------------------------* 4456* 4457* Purpose: transform a three-index batch of (jd|ic) integrals 4458* (with fixed index c) with B1 amplitudes to (j l^B|i c) 4459* (one-index transformed integrals with second index 4460* transformed) 4461* 4462* I^c(jl;i) = sum_a (j a|i c) * B1AMP(a l) 4463* 4464* IOPT = 1 : XOVOV contains (ia|jb) integrals packed 4465* 2 : full square matrix of (ia|jb) in XOVOV 4466* --> transform virtual index of the leading pair 4467* 4468* needed for CCSD part of CCQR_G 4469* 4470* Christof Haettig, August 1996 4471*=====================================================================* 4472#if defined (IMPLICIT_NONE) 4473 IMPLICIT NONE 4474#else 4475# include "implicit.h" 4476#endif 4477#include "priunit.h" 4478#include "ccsdsym.h" 4479#include "ccorb.h" 4480 4481 INTEGER ISYMJLI, ISYMB1, ISYOVOV, ISYMC, IOPT 4482 4483#if defined (SYS_CRAY) 4484 REAL XJLIC(*) ! dimension (NMAIJK(ISYMJLI)) 4485 REAL XOVOV(*) 4486 REAL B1AMP(*) ! dimension (NT1AM(ISYMB1)) 4487#else 4488 DOUBLE PRECISION XJLIC(*) ! dimension (NMAIJK(ISYMJLI)) 4489 DOUBLE PRECISION XOVOV(*) 4490 DOUBLE PRECISION B1AMP(*) ! dimension (NT1AM(ISYMB1)) 4491#endif 4492 4493 INTEGER ISYMDJ, ISYMCI, ISYMJL, ISYMI, ISYMJ, ISYML, ISYMD 4494 INTEGER NCI, NDJ, NDL, NJL, NJLI, NDJCI 4495 4496 INTEGER INDEX 4497 INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 4498 4499 CALL QENTER('CCG_TRANS2') 4500 4501* check symmetries: 4502 IF (MULD2H(ISYMJLI,ISYMC) .NE. MULD2H(ISYMB1,ISYOVOV)) 4503 & CALL QUIT('SYMMETRY MISMATCH IN CCG_TRANS2.') 4504 4505* initialize result: 4506 Call DZERO (Xjlic, NMAIJK(ISYMJLI)) 4507 4508 DO ISYMDJ = 1, NSYM 4509 ISYMCI = MULD2H(ISYMDJ,ISYOVOV) 4510 ISYMI = MULD2H(ISYMCI,ISYMC) 4511 4512 DO I = 1, NRHF(ISYMI) 4513 4514 NCI = IT1AM(ISYMC,ISYMI) + NVIR(ISYMC)*(I-1) + C 4515 4516 DO ISYMD = 1, NSYM 4517 ISYMJ = MULD2H(ISYMD,ISYMDJ) 4518 ISYML = MULD2H(ISYMD,ISYMB1) 4519 ISYMJL = MULD2H(ISYMJ,ISYML) 4520 4521 IF (IOPT .EQ. 1) THEN 4522 4523 IF (ISYMDJ .EQ. ISYMCI) THEN 4524 4525 DO L = 1, NRHF(ISYML) 4526 DO J = 1, NRHF(ISYMJ) 4527 NJL = IMATIJ(ISYMJ,ISYML) + NRHF(ISYMJ)*(L-1) + J 4528 NJLI = IMAIJK(ISYMJL,ISYMI) + NMATIJ(ISYMJL)*(I-1) + NJL 4529 DO D = 1, NVIR(ISYMD) 4530 NDJ = IT1AM(ISYMD,ISYMJ) + NVIR(ISYMD)*(J-1) + D 4531 NDL = IT1AM(ISYMD,ISYML) + NVIR(ISYMD)*(L-1) + D 4532 NDJCI = IT2AM(ISYMDJ,ISYMCI) + INDEX(NDJ,NCI) 4533 Xjlic(NJLI) = Xjlic(NJLI) + XOVOV(NDJCI) * B1AMP(NDL) 4534 END DO 4535 END DO 4536 END DO 4537 4538 ELSE IF (ISYMDJ .LT. ISYMCI) THEN 4539 4540 DO L = 1, NRHF(ISYML) 4541 DO J = 1, NRHF(ISYMJ) 4542 NJL = IMATIJ(ISYMJ,ISYML) + NRHF(ISYMJ)*(L-1) + J 4543 NJLI = IMAIJK(ISYMJL,ISYMI) + NMATIJ(ISYMJL)*(I-1) + NJL 4544 DO D = 1, NVIR(ISYMD) 4545 NDJ = IT1AM(ISYMD,ISYMJ) + NVIR(ISYMD)*(J-1) + D 4546 NDL = IT1AM(ISYMD,ISYML) + NVIR(ISYMD)*(L-1) + D 4547 NDJCI = IT2AM(ISYMDJ,ISYMCI) + NT1AM(ISYMDJ)*(NCI-1)+NDJ 4548 Xjlic(NJLI) = Xjlic(NJLI) + XOVOV(NDJCI) * B1AMP(NDL) 4549 END DO 4550 END DO 4551 END DO 4552 4553 ELSE IF (ISYMDJ .GT. ISYMCI) THEN 4554 4555 DO L = 1, NRHF(ISYML) 4556 DO J = 1, NRHF(ISYMJ) 4557 NJL = IMATIJ(ISYMJ,ISYML) + NRHF(ISYMJ)*(L-1) + J 4558 NJLI = IMAIJK(ISYMJL,ISYMI) + NMATIJ(ISYMJL)*(I-1) + NJL 4559 DO D = 1, NVIR(ISYMD) 4560 NDJ = IT1AM(ISYMD,ISYMJ) + NVIR(ISYMD)*(J-1) + D 4561 NDL = IT1AM(ISYMD,ISYML) + NVIR(ISYMD)*(L-1) + D 4562 NDJCI = IT2AM(ISYMCI,ISYMDJ) + NT1AM(ISYMCI)*(NDJ-1)+NCI 4563 Xjlic(NJLI) = Xjlic(NJLI) + XOVOV(NDJCI) * B1AMP(NDL) 4564 END DO 4565 END DO 4566 END DO 4567 4568 END IF 4569 4570 ELSE IF (IOPT .EQ. 2) THEN 4571 4572 DO L = 1, NRHF(ISYML) 4573 DO J = 1, NRHF(ISYMJ) 4574 NJL = IMATIJ(ISYMJ,ISYML) + NRHF(ISYMJ)*(L-1) + J 4575 NJLI = IMAIJK(ISYMJL,ISYMI) + NMATIJ(ISYMJL)*(I-1) + NJL 4576 DO D = 1, NVIR(ISYMD) 4577 NDJ = IT1AM(ISYMD,ISYMJ) + NVIR(ISYMD)*(J-1) + D 4578 NDL = IT1AM(ISYMD,ISYML) + NVIR(ISYMD)*(L-1) + D 4579 NDJCI = IT2SQ(ISYMDJ,ISYMCI) + NT1AM(ISYMDJ)*(NCI-1)+ NDJ 4580 Xjlic(NJLI) = Xjlic(NJLI) + XOVOV(NDJCI) * B1AMP(NDL) 4581 END DO 4582 END DO 4583 END DO 4584 4585 ELSE 4586 CALL QUIT('Illegal value for parameter IOPT in CCG_TRANS2.') 4587 END IF 4588 4589 END DO ! ISYMD 4590 4591 END DO ! J 4592 END DO ! ISYMDK 4593 4594 CALL QEXIT('CCG_TRANS2') 4595 4596 RETURN 4597 END 4598*---------------------------------------------------------------------* 4599c/* Deck CCG_TRANS4 */ 4600*=====================================================================* 4601 SUBROUTINE CCG_TRANS4(XIAJK,ISYMIAJK,XOVOV,ISYOVOV,T1AMP,ISYMT1) 4602*---------------------------------------------------------------------* 4603* 4604* Purpose: transform (ia|jb) integrals to (ia|jk) using T1(bk) 4605* 4606* I^j(ai;k) = sum_b (i a|j b) * T1AMP(b k) 4607* 4608* XOVOV contains (ia|jb) integrals as full matrix 4609* 4610* needed for CCSD part of CCQR_G 4611* 4612* Christof Haettig, August 1998 4613*=====================================================================* 4614#if defined (IMPLICIT_NONE) 4615 IMPLICIT NONE 4616#else 4617# include "implicit.h" 4618#endif 4619#include "priunit.h" 4620#include "ccsdsym.h" 4621#include "ccorb.h" 4622 4623 INTEGER ISYMIAJK, ISYMT1, ISYOVOV 4624 4625#if defined (SYS_CRAY) 4626 REAL XIAJK(*), XOVOV(*), T1AMP(*) 4627 REAL ONE, ZERO 4628#else 4629 DOUBLE PRECISION XIAJK(*), XOVOV(*), T1AMP(*) 4630 DOUBLE PRECISION ONE, ZERO 4631#endif 4632 PARAMETER (ZERO = 0.0d0, ONE = 1.0d0) 4633 4634 INTEGER ITAIKJ(8,8), NTAIKJ(8) 4635 INTEGER ISYMB, ISYMK, ISYMAI, ISYMBJ, ISYMAIK, ISYMJ 4636 INTEGER NBJ, KOFF1, KOFF2, KOFF3, ICOUNT, ISYM, NTOTB, NTOTAI 4637 4638C INTEGER INDEX 4639C INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 4640 4641 CALL QENTER('CCG_TRANS4') 4642 4643*---------------------------------------------------------------------* 4644* check symmetries & and set index arrays: 4645*---------------------------------------------------------------------* 4646 IF (ISYMIAJK .NE. MULD2H(ISYMT1,ISYOVOV)) THEN 4647 CALL QUIT('SYMMETRY MISMATCH IN CCG_TRANS2.') 4648 END IF 4649 4650 DO ISYM = 1, NSYM 4651 ICOUNT = 0 4652 DO ISYMAIK = 1, NSYM 4653 ISYMJ = MULD2H(ISYMAIK,ISYM) 4654 ITAIKJ(ISYMAIK,ISYMJ) = ICOUNT 4655 ICOUNT = ICOUNT + NT2BCD(ISYMAIK)*NRHF(ISYMJ) 4656 END DO 4657 NTAIKJ(ISYM) = ICOUNT 4658 END DO 4659 4660*---------------------------------------------------------------------* 4661* calculate (ia|jk) integrals and store as I^j(ai,k) 4662*---------------------------------------------------------------------* 4663 DO ISYMBJ = 1, NSYM 4664 ISYMAI = MULD2H(ISYOVOV,ISYMBJ) 4665 4666 DO ISYMB = 1, NSYM 4667 ISYMK = MULD2H(ISYMB,ISYMT1) 4668 ISYMAIK = MULD2H(ISYMAI,ISYMK) 4669 ISYMJ = MULD2H(ISYMBJ,ISYMB) 4670 4671 DO J = 1, NRHF(ISYMJ) 4672 4673 NBJ = IT1AM(ISYMB,ISYMJ) + NVIR(ISYMB)*(J-1) + 1 4674 KOFF1 = IT1AM(ISYMB,ISYMK) + 1 4675 KOFF2 = IT2SQ(ISYMAI,ISYMBJ) + NT1AM(ISYMAI)*(NBJ-1) + 1 4676 KOFF3 = ITAIKJ(ISYMAIK,ISYMJ) + NT2BCD(ISYMAIK)*(J-1) 4677 & + IT2BCD(ISYMAI,ISYMK) + 1 4678 4679 NTOTB = MAX(NVIR(ISYMB),1) 4680 NTOTAI = MAX(NT1AM(ISYMAI),1) 4681 4682 CALL DGEMM('N','N',NT1AM(ISYMAI),NRHF(ISYMK),NVIR(ISYMB), 4683 & ONE,XOVOV(KOFF2),NTOTAI,T1AMP(KOFF1),NTOTB, 4684 & ZERO,XIAJK(KOFF3),NTOTAI) 4685 4686 END DO 4687 4688 END DO 4689 END DO 4690 4691 CALL QEXIT('CCG_TRANS4') 4692 4693 RETURN 4694 END 4695*---------------------------------------------------------------------* 4696c/* Deck CCG_SORT1 */ 4697*=====================================================================* 4698 SUBROUTINE CCG_SORT1(X,Y,ISYMX,IOPT) 4699*---------------------------------------------------------------------* 4700* 4701* Purpose: resort integrals stored as X(ai,k,j) to 4702* 4703* IOPT = 0 -- Y(ai,j,k) 4704* 1 -- Y(ai,jk) 4705* 2 -- Y(ai,kj) 4706* 4 -- Y(a,jki) 4707* 5 -- Y(jki,a) 4708* 4709* Christof Haettig, August 1998 4710*=====================================================================* 4711#if defined (IMPLICIT_NONE) 4712 IMPLICIT NONE 4713#else 4714# include "implicit.h" 4715#endif 4716#include "priunit.h" 4717#include "ccsdsym.h" 4718#include "ccorb.h" 4719 4720 INTEGER ISYMX, IOPT 4721 4722#if defined (SYS_CRAY) 4723 REAL X(*), Y(*) 4724#else 4725 DOUBLE PRECISION X(*), Y(*) 4726#endif 4727 INTEGER ITAIKJ(8,8), NTAIKJ(8) 4728 INTEGER ISYM,ICOUNT,ISYMJ,ISYMAIK,ISYMAI,ISYMK,KOFF1,KOFF2,NJK 4729 INTEGER ISYMJK,ISYMAJ,ISYMI,ISYMA,ISYMAJK,ISYJIK,NAI,NJKI 4730 4731 CALL QENTER('CCG_SORT1') 4732 4733*---------------------------------------------------------------------* 4734* set index arrays for X integrals: 4735*---------------------------------------------------------------------* 4736 DO ISYM = 1, NSYM 4737 ICOUNT = 0 4738 DO ISYMAIK = 1, NSYM 4739 ISYMJ = MULD2H(ISYMAIK,ISYM) 4740 ITAIKJ(ISYMAIK,ISYMJ) = ICOUNT 4741 ICOUNT = ICOUNT + NT2BCD(ISYMAIK)*NRHF(ISYMJ) 4742 END DO 4743 NTAIKJ(ISYM) = ICOUNT 4744 END DO 4745 4746 IF (IOPT.EQ.0) THEN 4747*---------------------------------------------------------------------* 4748* resort X(ai,k,j) to Y(aj,k,i) 4749*---------------------------------------------------------------------* 4750 DO ISYMJ = 1, NSYM 4751 DO ISYMK = 1, NSYM 4752 4753 ISYMJK = MULD2H(ISYMJ,ISYMK) 4754 ISYMAI = MULD2H(ISYMJK,ISYMX) 4755 ISYMAIK = MULD2H(ISYMAI,ISYMK) 4756 4757 DO J = 1, NRHF(ISYMJ) 4758 DO K = 1, NRHF(ISYMK) 4759 4760 DO ISYMI = 1, NSYM 4761 4762 ISYMA = MULD2H(ISYMI,ISYMAI) 4763 ISYMAJ = MULD2H(ISYMA,ISYMJ) 4764 ISYMAJK = MULD2H(ISYMK,ISYMAJ) 4765 4766 KOFF1 = ITAIKJ(ISYMAIK,ISYMJ) + NT2BCD(ISYMAIK)*(J-1) 4767 & + IT2BCD(ISYMAI,ISYMK) + NT1AM(ISYMAI)*(K-1) 4768 & + IT1AM(ISYMA,ISYMI) + 1 4769 4770 KOFF2 = ITAIKJ(ISYMAJK,ISYMI) 4771 & + IT2BCD(ISYMAJ,ISYMK) + NT1AM(ISYMAJ)*(K-1) 4772 & + IT1AM(ISYMA,ISYMJ) + NVIR(ISYMA)*(J-1) + 1 4773 4774 DO I = 1, NRHF(ISYMI) 4775 4776 CALL DCOPY(NVIR(ISYMA),X(KOFF1),1,Y(KOFF2),1) 4777 4778 KOFF2 = KOFF2 + NT2BCD(ISYMAJK) 4779 KOFF1 = KOFF1 + NVIR(ISYMA) 4780 4781 END DO 4782 END DO 4783 END DO 4784 END DO 4785 4786 END DO 4787 END DO 4788 4789 ELSE IF (IOPT.EQ.1 .OR. IOPT.EQ.2) THEN 4790*---------------------------------------------------------------------* 4791* resort X(ai,k,j) to Y(ai,jk) / Y(ai,kj) 4792*---------------------------------------------------------------------* 4793 DO ISYMJ = 1, NSYM 4794 DO ISYMK = 1, NSYM 4795 4796 ISYMJK = MULD2H(ISYMJ,ISYMK) 4797 ISYMAI = MULD2H(ISYMJK,ISYMX) 4798 ISYMAIK = MULD2H(ISYMAI,ISYMK) 4799 4800 DO J = 1, NRHF(ISYMJ) 4801 DO K = 1, NRHF(ISYMK) 4802 4803 IF (IOPT.EQ.1) THEN 4804 NJK = IMATIJ(ISYMJ,ISYMK) + NRHF(ISYMJ)*(K-1) + J 4805 ELSE 4806 NJK = IMATIJ(ISYMK,ISYMJ) + NRHF(ISYMK)*(J-1) + K 4807 END IF 4808 4809 KOFF1 = ITAIKJ(ISYMAIK,ISYMJ) + NT2BCD(ISYMAIK)*(J-1) 4810 & + IT2BCD(ISYMAI,ISYMK) + NT1AM(ISYMAI)*(K-1) + 1 4811 4812 KOFF2 = ISAIKL(ISYMAI,ISYMJK) + NT1AM(ISYMAI)*(NJK-1) + 1 4813 4814 CALL DCOPY(NT1AM(ISYMAI),X(KOFF1),1,Y(KOFF2),1) 4815 4816 END DO 4817 END DO 4818 4819 END DO 4820 END DO 4821 4822 ELSE IF (IOPT.EQ.4) THEN 4823*---------------------------------------------------------------------* 4824* resort X(ai,k,j) to Y(a,jki) 4825*---------------------------------------------------------------------* 4826 DO ISYMJ = 1, NSYM 4827 DO ISYMK = 1, NSYM 4828 DO ISYMI = 1, NSYM 4829 4830 ISYMJK = MULD2H(ISYMJ,ISYMK) 4831 ISYJIK = MULD2H(ISYMJK,ISYMI) 4832 ISYMAI = MULD2H(ISYMJK,ISYMX) 4833 ISYMAIK = MULD2H(ISYMAI,ISYMK) 4834 ISYMA = MULD2H(ISYMI,ISYMAI) 4835 4836 DO J = 1, NRHF(ISYMJ) 4837 DO K = 1, NRHF(ISYMK) 4838 DO I = 1, NRHF(ISYMI) 4839 4840 NJK = IMATIJ(ISYMJ,ISYMK) + NRHF(ISYMJ)*(K-1) + J 4841 NJKI = IMAIJK(ISYMJK,ISYMI) + NMATIJ(ISYMJK)*(I-1) + NJK 4842 NAI = IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I-1) + 1 4843 4844 KOFF1 = ITAIKJ(ISYMAIK,ISYMJ) + NT2BCD(ISYMAIK)*(J-1) 4845 & + IT2BCD(ISYMAI,ISYMK) + NT1AM(ISYMAI)*(K-1) + NAI 4846 4847 KOFF2 = ISJIKA(ISYJIK,ISYMA) + NVIR(ISYMA)*(NJKI-1) + 1 4848 4849 CALL DCOPY(NVIR(ISYMA),X(KOFF1),1,Y(KOFF2),1) 4850 4851 END DO 4852 END DO 4853 END DO 4854 4855 END DO 4856 END DO 4857 END DO 4858 4859 ELSE IF (IOPT.EQ.5) THEN 4860*---------------------------------------------------------------------* 4861* resort X(ai,k,j) to Y(jki,a) 4862*---------------------------------------------------------------------* 4863 DO ISYMJ = 1, NSYM 4864 DO ISYMK = 1, NSYM 4865 DO ISYMI = 1, NSYM 4866 4867 ISYMJK = MULD2H(ISYMJ,ISYMK) 4868 ISYJIK = MULD2H(ISYMJK,ISYMI) 4869 ISYMAI = MULD2H(ISYMJK,ISYMX) 4870 ISYMAIK = MULD2H(ISYMAI,ISYMK) 4871 ISYMA = MULD2H(ISYMI,ISYMAI) 4872 4873 DO J = 1, NRHF(ISYMJ) 4874 DO K = 1, NRHF(ISYMK) 4875 DO I = 1, NRHF(ISYMI) 4876 4877 NJK = IMATIJ(ISYMJ,ISYMK) + NRHF(ISYMJ)*(K-1) + J 4878 NJKI = IMAIJK(ISYMJK,ISYMI) + NMATIJ(ISYMJK)*(I-1) + NJK 4879 NAI = IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I-1) + 1 4880 4881 KOFF1 = ITAIKJ(ISYMAIK,ISYMJ) + NT2BCD(ISYMAIK)*(J-1) 4882 & + IT2BCD(ISYMAI,ISYMK) + NT1AM(ISYMAI)*(K-1) + NAI 4883 4884 KOFF2 = ISJIKA(ISYJIK,ISYMA) + NJKI 4885 4886 CALL DCOPY(NVIR(ISYMA),X(KOFF1),1,Y(KOFF2),NMAIJK(ISYJIK)) 4887 4888 END DO 4889 END DO 4890 END DO 4891 4892 END DO 4893 END DO 4894 END DO 4895 4896*---------------------------------------------------------------------* 4897 ELSE 4898 CALL QUIT('Illegal option in CCG_SORT1.') 4899 END IF 4900 4901 CALL QEXIT('CCG_SORT1') 4902 4903 RETURN 4904 END 4905*---------------------------------------------------------------------* 4906c/* Deck CCG_LXD */ 4907*=====================================================================* 4908 SUBROUTINE CCG_LXD(FTWO,ISYFCK,DENS,ISYDNS,XLIAJB,ISYOVOV,IOPT) 4909*---------------------------------------------------------------------* 4910* 4911* Purpose: calculate contraction of packed (squared) 4912* L(ia,jb) integrals with a density matrix D(jb) 4913* 4914* FTWO(ai) = sum(bj) L(ia,jb) * DENS(jb) 4915* 4916* IOPT = 0 : packed integrals 4917* IOPT = 1 : squared integrals 4918* 4919* needed, e.g., for the two electron contribution to 4920* the first order responses of the Fock matrix 4921* 4922* Christof Haettig, August 1996 4923*=====================================================================* 4924#if defined (IMPLICIT_NONE) 4925 IMPLICIT NONE 4926#else 4927# include "implicit.h" 4928#endif 4929#include "priunit.h" 4930#include "ccorb.h" 4931#include "ccsdsym.h" 4932 4933 INTEGER ISYFCK, ISYDNS, ISYOVOV, IOPT, KOFF1 4934 4935#if defined (SYS_CRAY) 4936 REAL FTWO(*), DENS(*), XLIAJB(*) 4937 REAL ZERO, ONE 4938#else 4939 DOUBLE PRECISION FTWO(*), DENS(*), XLIAJB(*) 4940 DOUBLE PRECISION ZERO, ONE 4941#endif 4942 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0) 4943 4944 INTEGER ISYMAI, ISYMBJ, NAI, NBJ, NAIBJ, NTOTAI, NTOTBJ 4945 4946 INTEGER INDEX 4947 INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 4948 4949 CALL QENTER('CCG_LXD') 4950 4951* check symmetries: 4952 IF (ISYFCK .NE. MULD2H(ISYDNS,ISYOVOV)) THEN 4953 CALL QUIT('Symmetry mismatch in CCG_LXD.') 4954 END IF 4955 4956* initialize result vector: 4957 CALL DZERO(FTWO, NT1AM(ISYFCK)) 4958 4959* calculate contraction FTWO(ia) = sum_bj L(ia,jb) * D(jb): 4960 ISYMAI = ISYFCK 4961 ISYMBJ = ISYDNS 4962 4963 IF (IOPT.EQ.0 .AND. ISYMAI. EQ. ISYMBJ) THEN 4964 4965 DO NAI = 1, NT1AM(ISYMAI) 4966 DO NBJ = 1, NT1AM(ISYMBJ) 4967 4968 NAIBJ = IT2AM(ISYMAI,ISYMBJ) + INDEX(NAI,NBJ) 4969 FTWO(NAI) = FTWO(NAI) + XLIAJB(NAIBJ) * DENS(NBJ) 4970 4971 END DO 4972 END DO 4973 4974 ELSE IF (IOPT.EQ.0 .AND. ISYMAI .GT. ISYMBJ) THEN 4975 4976 KOFF1 = IT2AM(ISYMBJ,ISYMAI) + 1 4977 NTOTBJ = MAX(1,NT1AM(ISYMBJ)) 4978 4979 CALL DGEMV('T',NT1AM(ISYMBJ),NT1AM(ISYMAI),ONE,XLIAJB(KOFF1), 4980 & NTOTBJ,DENS,1,ZERO,FTWO,1) 4981 4982 ELSE IF (IOPT.EQ.0 .AND. ISYMAI .LT. ISYMBJ) THEN 4983 4984 KOFF1 = IT2AM(ISYMAI,ISYMBJ) + 1 4985 NTOTAI = MAX(1,NT1AM(ISYMAI)) 4986 4987 CALL DGEMV('N',NT1AM(ISYMAI),NT1AM(ISYMBJ),ONE,XLIAJB(KOFF1), 4988 & NTOTAI,DENS,1,ZERO,FTWO,1) 4989 4990 ELSE IF (IOPT.EQ.1) THEN 4991 4992 KOFF1 = IT2SQ(ISYMAI,ISYMBJ) + 1 4993 NTOTAI = MAX(1,NT1AM(ISYMAI)) 4994 4995 CALL DGEMV('N',NT1AM(ISYMAI),NT1AM(ISYMBJ),ONE,XLIAJB(KOFF1), 4996 & NTOTAI,DENS,1,ZERO,FTWO,1) 4997 4998 END IF 4999 5000 CALL QEXIT('CCG_LXD') 5001 5002 RETURN 5003 END 5004*---------------------------------------------------------------------* 5005*=====================================================================* 5006 SUBROUTINE CCSDT_GMAT_NODDY(LISTL,IDLSTL,LISTB,IDLSTB, 5007 & LISTC,IDLSTC, 5008 & OMEGA1,OMEGA2,WORK,LWORK) 5009*---------------------------------------------------------------------* 5010* 5011* Purpose: compute triples contribution to G matrix transformation 5012* 5013* (G T^B T^C)^eff_1,2 = (G T^B T^C)_1,2(CCSD) 5014* + (G T^B T^C)_1,2(L_3) 5015* 5016* Written by Christof Haettig, April 2002 5017* based on CCSDT_FMAT_NODDY 5018* 5019*=====================================================================* 5020#if defined (IMPLICIT_NONE) 5021 IMPLICIT NONE 5022#else 5023# include "implicit.h" 5024#endif 5025#include "dummy.h" 5026#include "priunit.h" 5027#include "ccsdinp.h" 5028#include "maxorb.h" 5029#include "ccsdsym.h" 5030#include "ccfield.h" 5031#include "ccorb.h" 5032 5033 LOGICAL LOCDBG 5034 PARAMETER (LOCDBG=.FALSE.) 5035 5036 INTEGER ISYM0 5037 PARAMETER (ISYM0 = 1) 5038 5039 CHARACTER*3 LISTL, LISTB, LISTC 5040 INTEGER LWORK, IDLSTL, IDLSTB, IDLSTC 5041 5042#if defined (SYS_CRAY) 5043 REAL WORK(LWORK) 5044 REAL OMEGA1(NT1AMX), OMEGA2(NT2AMX) 5045 REAL DDOT, TWO, XNORMVAL, FREQL 5046#else 5047 DOUBLE PRECISION WORK(LWORK) 5048 DOUBLE PRECISION OMEGA1(NT1AMX), OMEGA2(NT2AMX) 5049 DOUBLE PRECISION DDOT, TWO, XNORMVAL, FREQL 5050#endif 5051 PARAMETER (TWO = 2.0D0) 5052 5053 CHARACTER*10 MODEL 5054 INTEGER KXINT, KXIAJB, KYIAJB, KT1AMP0, KLAMP0, KLAMH0, KEND1, 5055 & KFOCK0, LWRK1, KLAMPB, KLAMHB, KLAMPC, KLAMHC, KT1AMPB, 5056 & KT1AMPC, KINT1T0, KINT2T0, KINT1SBC, KINT2SBC, 5057 & KBIOVVO, KBIOOVV, KBIOOOO, KBIVVVV, 5058 & KCIOVVO, KCIOOVV, KCIOOOO, KCIVVVV, 5059 & KBCIOVVO, KBCIOOVV, KBCIOOOO, KBCIVVVV, 5060 & KOME1, KOME2, KL1AM, KL2AM, KL3AM, KT2AM, KFOCKD, 5061 & KSCR1, KDUM, KFIELD, KFIELDAO 5062 INTEGER ISYMD, ILLL, IDEL, ISYDIS, IJ, NIJ, LUSIFC, INDEX, 5063 & ISYMC, ISYMB, LUFOCK, KEND2, LWRK2, IOPT, ILSTSYM, ISYML 5064 5065 INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 5066 5067 KDUM = 1 5068 CALL QENTER('CCSDT_GMAT_NODDY') 5069 5070 IF (DIRECT) CALL QUIT('CCSDT_GMAT_NODDY: DIRECT NOT IMPLEMENTED') 5071 5072*---------------------------------------------------------------------* 5073* Memory allocation: 5074*---------------------------------------------------------------------* 5075 KSCR1 = 1 5076 KFOCKD = KSCR1 + NT1AMX 5077 KEND1 = KFOCKD + NORBT 5078 5079 KFOCK0 = KEND1 5080 KEND1 = KFOCK0 + NORBT*NORBT 5081 5082 IF (NONHF) THEN 5083 KFIELD = KEND1 5084 KFIELDAO = KFIELD + NORBT*NORBT 5085 KEND1 = KFIELDAO + NORBT*NORBT 5086 END IF 5087 5088 KT1AMP0 = KEND1 5089 KLAMP0 = KT1AMP0 + NT1AMX 5090 KLAMH0 = KLAMP0 + NLAMDT 5091 KEND1 = KLAMH0 + NLAMDT 5092 5093 KT1AMPB = KEND1 5094 KLAMPB = KT1AMPB + NT1AMX 5095 KLAMHB = KLAMPB + NLAMDT 5096 KEND1 = KLAMHB + NLAMDT 5097 5098 KT1AMPC = KEND1 5099 KLAMPC = KT1AMPC + NT1AMX 5100 KLAMHC = KLAMPC + NLAMDT 5101 KEND1 = KLAMHC + NLAMDT 5102 5103 KXIAJB = KEND1 5104 KYIAJB = KXIAJB + NT1AMX*NT1AMX 5105 KEND1 = KYIAJB + NT1AMX*NT1AMX 5106 5107 KINT1T0 = KEND1 5108 KINT2T0 = KINT1T0 + NT1AMX*NVIRT*NVIRT 5109 KEND1 = KINT2T0 + NRHFT*NRHFT*NT1AMX 5110 5111 KBIOVVO = KEND1 5112 KBIOOVV = KBIOVVO + NRHFT*NVIRT*NVIRT*NRHFT 5113 KBIOOOO = KBIOOVV + NRHFT*NVIRT*NVIRT*NRHFT 5114 KBIVVVV = KBIOOOO + NRHFT*NRHFT*NRHFT*NRHFT 5115 KEND1 = KBIVVVV + NVIRT*NVIRT*NVIRT*NVIRT 5116 5117 KCIOVVO = KEND1 5118 KCIOOVV = KCIOVVO + NRHFT*NVIRT*NVIRT*NRHFT 5119 KCIOOOO = KCIOOVV + NRHFT*NVIRT*NVIRT*NRHFT 5120 KCIVVVV = KCIOOOO + NRHFT*NRHFT*NRHFT*NRHFT 5121 KEND1 = KCIVVVV + NVIRT*NVIRT*NVIRT*NVIRT 5122 5123 KBCIOVVO = KEND1 5124 KBCIOOVV = KBCIOVVO + NRHFT*NVIRT*NVIRT*NRHFT 5125 KBCIOOOO = KBCIOOVV + NRHFT*NVIRT*NVIRT*NRHFT 5126 KBCIVVVV = KBCIOOOO + NRHFT*NRHFT*NRHFT*NRHFT 5127 KEND1 = KBCIVVVV + NVIRT*NVIRT*NVIRT*NVIRT 5128 5129 KINT1SBC = KEND1 5130 KINT2SBC = KINT1SBC + NT1AMX*NVIRT*NVIRT 5131 KEND1 = KINT2SBC + NRHFT*NRHFT*NT1AMX 5132 5133 KOME1 = KEND1 5134 KOME2 = KOME1 + NT1AMX 5135 KEND1 = KOME2 + NT1AMX*NT1AMX 5136 5137 KL1AM = KEND1 5138 KL2AM = KL1AM + NT1AMX 5139 KL3AM = KL2AM + NT1AMX*NT1AMX 5140 KEND1 = KL3AM + NT1AMX*NT1AMX*NT1AMX 5141 5142 LWRK1 = LWORK - KEND1 5143 IF (LWRK1 .LT. 0) THEN 5144 CALL QUIT('Insufficient space in CCSDT_GMAT_NODDY') 5145 ENDIF 5146 5147*---------------------------------------------------------------------* 5148* Read SCF orbital energies from file: 5149*---------------------------------------------------------------------* 5150 LUSIFC = -1 5151 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY, 5152 & .FALSE.) 5153 REWIND LUSIFC 5154 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 5155 READ (LUSIFC) 5156 READ (LUSIFC) (WORK(KFOCKD+I-1), I=1,NORBT) 5157 CALL GPCLOSE(LUSIFC,'KEEP') 5158 5159*---------------------------------------------------------------------* 5160* Get zeroth-order Lambda matrices: 5161*---------------------------------------------------------------------* 5162 IOPT = 1 5163 Call CC_RDRSP('R0',0,ISYM0,IOPT,MODEL,WORK(KT1AMP0),WORK(KDUM)) 5164 5165 Call LAMMAT(WORK(KLAMP0),WORK(KLAMH0),WORK(KT1AMP0), 5166 & WORK(KEND1),LWRK1) 5167 5168*---------------------------------------------------------------------* 5169* Get response Lambda matrices: 5170*---------------------------------------------------------------------* 5171 ISYMB = ILSTSYM(LISTB,IDLSTB) 5172 IOPT = 1 5173 Call CC_RDRSP(LISTB,IDLSTB,ISYMB,IOPT,MODEL, 5174 & WORK(KT1AMPB),WORK(KDUM)) 5175 5176 CALL CCLR_LAMTRA(WORK(KLAMP0),WORK(KLAMPB), 5177 & WORK(KLAMH0),WORK(KLAMHB),WORK(KT1AMPB),ISYMB) 5178 5179 ISYMC = ILSTSYM(LISTC,IDLSTC) 5180 IOPT = 1 5181 Call CC_RDRSP(LISTC,IDLSTC,ISYMC,IOPT,MODEL, 5182 & WORK(KT1AMPC),WORK(KDUM)) 5183 5184 CALL CCLR_LAMTRA(WORK(KLAMP0),WORK(KLAMPC), 5185 & WORK(KLAMH0),WORK(KLAMHC),WORK(KT1AMPC),ISYMC) 5186 5187*---------------------------------------------------------------------* 5188* read zeroth-order AO Fock matrix from file: 5189*---------------------------------------------------------------------* 5190 LUFOCK = -1 5191 CALL GPOPEN(LUFOCK,'CC_FCKH','OLD',' ','UNFORMATTED',IDUMMY, 5192 & .FALSE.) 5193 REWIND(LUFOCK) 5194 READ(LUFOCK) (WORK(KFOCK0-1+I),I=1,N2BST(ISYM0)) 5195 CALL GPCLOSE(LUFOCK,'KEEP') 5196 5197 CALL CC_FCKMO(WORK(KFOCK0),WORK(KLAMP0),WORK(KLAMH0), 5198 & WORK(KEND1),LWRK1,ISYM0,ISYM0,ISYM0) 5199 5200*---------------------------------------------------------------------* 5201* Get the one electron integrals from the fields 5202*---------------------------------------------------------------------* 5203 IF (NONHF) THEN 5204 CALL DZERO(WORK(KFIELDAO),NORBT*NORBT) 5205 DO I = 1, NFIELD 5206 CALL CC_ONEP(WORK(KFIELDAO),WORK(KEND1),LWRK1, 5207 * EFIELD(I),1,LFIELD(I)) 5208 ENDDO 5209 CALL DCOPY(NORBT*NORBT,WORK(KFIELDAO),1,WORK(KFIELD),1) 5210 5211 ! calculate external field in zero-order lambda basis 5212 CALL CC_FCKMO(WORK(KFIELD),WORK(KLAMP0),WORK(KLAMH0), 5213 * WORK(KEND1),LWRK1,ISYM0,ISYM0,ISYM0) 5214 ENDIF 5215 5216*---------------------------------------------------------------------* 5217* Compute integrals needed for the following contributions: 5218*---------------------------------------------------------------------* 5219 5220 CALL DZERO(WORK(KBIOVVO),NRHFT*NVIRT*NVIRT*NRHFT) 5221 CALL DZERO(WORK(KBIOOVV),NRHFT*NVIRT*NVIRT*NRHFT) 5222 CALL DZERO(WORK(KBIOOOO),NRHFT*NRHFT*NRHFT*NRHFT) 5223 CALL DZERO(WORK(KBIVVVV),NVIRT*NVIRT*NVIRT*NVIRT) 5224 5225 CALL DZERO(WORK(KCIOVVO),NRHFT*NVIRT*NVIRT*NRHFT) 5226 CALL DZERO(WORK(KCIOOVV),NRHFT*NVIRT*NVIRT*NRHFT) 5227 CALL DZERO(WORK(KCIOOOO),NRHFT*NRHFT*NRHFT*NRHFT) 5228 CALL DZERO(WORK(KCIVVVV),NVIRT*NVIRT*NVIRT*NVIRT) 5229 5230 CALL DZERO(WORK(KBCIOVVO),NRHFT*NVIRT*NVIRT*NRHFT) 5231 CALL DZERO(WORK(KBCIOOVV),NRHFT*NVIRT*NVIRT*NRHFT) 5232 CALL DZERO(WORK(KBCIOOOO),NRHFT*NRHFT*NRHFT*NRHFT) 5233 CALL DZERO(WORK(KBCIVVVV),NVIRT*NVIRT*NVIRT*NVIRT) 5234 5235 CALL DZERO(WORK(KXIAJB),NT1AMX*NT1AMX) 5236 CALL DZERO(WORK(KYIAJB),NT1AMX*NT1AMX) 5237 5238 CALL DZERO(WORK(KINT1T0),NT1AMX*NVIRT*NVIRT) 5239 CALL DZERO(WORK(KINT2T0),NRHFT*NRHFT*NT1AMX) 5240 5241 CALL DZERO(WORK(KINT1SBC),NT1AMX*NVIRT*NVIRT) 5242 CALL DZERO(WORK(KINT2SBC),NRHFT*NRHFT*NT1AMX) 5243 5244 DO ISYMD = 1, NSYM 5245 DO ILLL = 1,NBAS(ISYMD) 5246 IDEL = IBAS(ISYMD) + ILLL 5247 ISYDIS = MULD2H(ISYMD,ISYMOP) 5248 5249C ---------------------------- 5250C Work space allocation no. 2. 5251C ---------------------------- 5252 KXINT = KEND1 5253 KEND2 = KXINT + NDISAO(ISYDIS) 5254 LWRK2 = LWORK - KEND2 5255 IF (LWRK2 .LT. 0) THEN 5256 WRITE(LUPRI,*) 'Need : ',KEND2,'Available : ',LWORK 5257 CALL QUIT('Insufficient space in CCSDT_GMAT_NODDY') 5258 ENDIF 5259 5260C --------------------------- 5261C Read in batch of integrals. 5262C --------------------------- 5263 CALL CCRDAO(WORK(KXINT),IDEL,1,WORK(KEND2),LWRK2, 5264 * WORK(KDUM),DIRECT) 5265 5266C ---------------------------------- 5267C Calculate integrals needed in CC3: 5268C ---------------------------------- 5269 CALL CCSDT_TRAN1(WORK(KINT1T0),WORK(KINT2T0), 5270 & WORK(KLAMP0),WORK(KLAMH0), 5271 & WORK(KXINT),IDEL) 5272 5273 CALL CC3_TRAN2(WORK(KXIAJB),WORK(KYIAJB), 5274 & WORK(KLAMP0),WORK(KLAMH0), 5275 & WORK(KXINT),IDEL) 5276 5277 ! XINT1S = XINT1S + (C-barB K-barC|B D) 5278 ! XINT2S = XINT2S + (C-barB K-barC|L J) 5279 CALL CCSDT_TRAN3_R(WORK(KINT1SBC),WORK(KINT2SBC), 5280 & WORK(KLAMP0),WORK(KLAMH0), 5281 & WORK(KLAMPB),WORK(KLAMHC), 5282 & WORK(KLAMP0),WORK(KLAMH0), 5283 & WORK(KXINT),IDEL) 5284 5285 ! XINT1S = XINT1S + (C-barB K|B-barC D) 5286 ! XINT2S = XINT2S + (C-barB K|L J-barC) 5287 CALL CCSDT_TRAN3_R(WORK(KINT1SBC),WORK(KINT2SBC), 5288 & WORK(KLAMP0),WORK(KLAMH0), 5289 & WORK(KLAMPB),WORK(KLAMH0), 5290 & WORK(KLAMPC),WORK(KLAMHC), 5291 & WORK(KXINT),IDEL) 5292 5293 ! XINT1S = XINT1S + (C-barC K-barB|B D) 5294 ! XINT2S = XINT2S + (C-barC K-barB|L J) 5295 CALL CCSDT_TRAN3_R(WORK(KINT1SBC),WORK(KINT2SBC), 5296 & WORK(KLAMP0),WORK(KLAMH0), 5297 & WORK(KLAMPC),WORK(KLAMHB), 5298 & WORK(KLAMP0),WORK(KLAMH0), 5299 & WORK(KXINT),IDEL) 5300 5301 ! XINT1S = XINT1S + (C K-barB|B-barC D) 5302 ! XINT2S = XINT2S + (C K-barB|L J-barC) 5303 CALL CCSDT_TRAN3_R(WORK(KINT1SBC),WORK(KINT2SBC), 5304 & WORK(KLAMP0),WORK(KLAMH0), 5305 & WORK(KLAMP0),WORK(KLAMHB), 5306 & WORK(KLAMPC),WORK(KLAMHC), 5307 & WORK(KXINT),IDEL) 5308 5309 ! XINT1S = XINT1S + (C-barC K|B-barB D) 5310 ! XINT2S = XINT2S + (C-barC K|L J-barB) 5311 CALL CCSDT_TRAN3_R(WORK(KINT1SBC),WORK(KINT2SBC), 5312 & WORK(KLAMP0),WORK(KLAMH0), 5313 & WORK(KLAMPC),WORK(KLAMH0), 5314 & WORK(KLAMPB),WORK(KLAMHB), 5315 & WORK(KXINT),IDEL) 5316 5317 ! XINT1S = XINT1S + (C K-barC|B-barB D) 5318 ! XINT2S = XINT2S + (C K-barC|L J-barB) 5319 CALL CCSDT_TRAN3_R(WORK(KINT1SBC),WORK(KINT2SBC), 5320 & WORK(KLAMP0),WORK(KLAMH0), 5321 & WORK(KLAMP0),WORK(KLAMHC), 5322 & WORK(KLAMPB),WORK(KLAMHB), 5323 & WORK(KXINT),IDEL) 5324 5325C ---------------------------------------------- 5326C (OV|VO)-B (OO|VV)-B (OO|OO)-B (VV|VV)-B 5327C ---------------------------------------------- 5328 CALL CCFOP_TRAN1_R(WORK(KBIOVVO),WORK(KBIOOVV), 5329 & WORK(KBIOOOO),WORK(KBIVVVV), 5330 & WORK(KLAMP0),WORK(KLAMH0), 5331 & WORK(KLAMPB),WORK(KLAMHB), 5332 & WORK(KLAMP0),WORK(KLAMH0), 5333 & WORK(KXINT),IDEL) 5334 5335 CALL CCFOP_TRAN1_R(WORK(KBIOVVO),WORK(KBIOOVV), 5336 & WORK(KBIOOOO),WORK(KBIVVVV), 5337 & WORK(KLAMP0),WORK(KLAMH0), 5338 & WORK(KLAMP0),WORK(KLAMH0), 5339 & WORK(KLAMPB),WORK(KLAMHB), 5340 & WORK(KXINT),IDEL) 5341 5342C ---------------------------------------------- 5343C (OV|VO)-C (OO|VV)-C (OO|OO)-C (VV|VV)-C 5344C ---------------------------------------------- 5345 CALL CCFOP_TRAN1_R(WORK(KCIOVVO),WORK(KCIOOVV), 5346 & WORK(KCIOOOO),WORK(KCIVVVV), 5347 & WORK(KLAMP0),WORK(KLAMH0), 5348 & WORK(KLAMPC),WORK(KLAMHC), 5349 & WORK(KLAMP0),WORK(KLAMH0), 5350 & WORK(KXINT),IDEL) 5351 5352 CALL CCFOP_TRAN1_R(WORK(KCIOVVO),WORK(KCIOOVV), 5353 & WORK(KCIOOOO),WORK(KCIVVVV), 5354 & WORK(KLAMP0),WORK(KLAMH0), 5355 & WORK(KLAMP0),WORK(KLAMH0), 5356 & WORK(KLAMPC),WORK(KLAMHC), 5357 & WORK(KXINT),IDEL) 5358 5359C ---------------------------------------------- 5360C (OV|VO)-BC (OO|VV)-BC (OO|OO)-BC (VV|VV)-BC 5361C ---------------------------------------------- 5362 CALL CCFOP_TRAN1_R(WORK(KBCIOVVO),WORK(KBCIOOVV), 5363 & WORK(KBCIOOOO),WORK(KBCIVVVV), 5364 & WORK(KLAMP0),WORK(KLAMH0), 5365 & WORK(KLAMPB),WORK(KLAMHB), 5366 & WORK(KLAMPC),WORK(KLAMHC), 5367 & WORK(KXINT),IDEL) 5368 5369 CALL CCFOP_TRAN1_R(WORK(KBCIOVVO),WORK(KBCIOOVV), 5370 & WORK(KBCIOOOO),WORK(KBCIVVVV), 5371 & WORK(KLAMP0),WORK(KLAMH0), 5372 & WORK(KLAMPC),WORK(KLAMHC), 5373 & WORK(KLAMPB),WORK(KLAMHB), 5374 & WORK(KXINT),IDEL) 5375 5376 END DO 5377 END DO 5378 5379 if (locdbg) then 5380 write(lupri,*) 'norm^2(OVVO-B):', 5381 & ddot(NRHFT*NVIRT*NVIRT*NRHFT,WORK(KBIOVVO),1,WORK(KBIOVVO),1) 5382 write(lupri,*) 'norm^2(OOVV-B):', 5383 & ddot(NRHFT*NVIRT*NVIRT*NRHFT,WORK(KBIOOVV),1,WORK(KBIOOVV),1) 5384 write(lupri,*) 'norm^2(OOOO-B):', 5385 & ddot(NRHFT*NRHFT*NRHFT*NRHFT,WORK(KBIOOOO),1,WORK(KBIOOOO),1) 5386 write(lupri,*) 'norm^2(VVVV-B):', 5387 & ddot(NVIRT*NVIRT*NVIRT*NVIRT,WORK(KBIVVVV),1,WORK(KBIVVVV),1) 5388 5389 write(lupri,*) 'norm^2(OVVO-C):', 5390 & ddot(NRHFT*NVIRT*NVIRT*NRHFT,WORK(KCIOVVO),1,WORK(KCIOVVO),1) 5391 write(lupri,*) 'norm^2(OOVV-C):', 5392 & ddot(NRHFT*NVIRT*NVIRT*NRHFT,WORK(KCIOOVV),1,WORK(KCIOOVV),1) 5393 write(lupri,*) 'norm^2(OOOO-C):', 5394 & ddot(NRHFT*NRHFT*NRHFT*NRHFT,WORK(KCIOOOO),1,WORK(KCIOOOO),1) 5395 write(lupri,*) 'norm^2(VVVV-C):', 5396 & ddot(NVIRT*NVIRT*NVIRT*NVIRT,WORK(KCIVVVV),1,WORK(KCIVVVV),1) 5397 5398 write(lupri,*) 'norm^2(OVVO-BC):', 5399 & ddot(NRHFT*NVIRT*NVIRT*NRHFT,WORK(KBCIOVVO),1,WORK(KBCIOVVO),1) 5400 write(lupri,*) 'norm^2(OOVV-BC):', 5401 & ddot(NRHFT*NVIRT*NVIRT*NRHFT,WORK(KBCIOOVV),1,WORK(KBCIOOVV),1) 5402 write(lupri,*) 'norm^2(OOOO-BC):', 5403 & ddot(NRHFT*NRHFT*NRHFT*NRHFT,WORK(KBCIOOOO),1,WORK(KBCIOOOO),1) 5404 write(lupri,*) 'norm^2(VVVV-BC):', 5405 & ddot(NVIRT*NVIRT*NVIRT*NVIRT,WORK(KBCIVVVV),1,WORK(KBCIVVVV),1) 5406 5407 write(lupri,*) 'norm^2(XIAJB):', 5408 & ddot(NT1AMX*NT1AMX,WORK(KXIAJB),1,WORK(KXIAJB),1) 5409 write(lupri,*) 'norm^2(YIAJB):', 5410 & ddot(NT1AMX*NT1AMX,WORK(KYIAJB),1,WORK(KYIAJB),1) 5411 5412 write(lupri,*) 'norm^2(INT1T0):', 5413 & ddot(NT1AMX*NVIRT*NVIRT,WORK(KINT1T0),1,WORK(KINT1T0),1) 5414 write(lupri,*) 'norm^2(INT2T0):', 5415 & ddot(NRHFT*NRHFT*NT1AMX,WORK(KINT2T0),1,WORK(KINT2T0),1) 5416 5417 write(lupri,*) 'norm^2(INT1SBC):', 5418 & ddot(NT1AMX*NVIRT*NVIRT,WORK(KINT1SBC),1,WORK(KINT1SBC),1) 5419 write(lupri,*) 'norm^2(INT2SBC):', 5420 & ddot(NRHFT*NRHFT*NT1AMX,WORK(KINT2SBC),1,WORK(KINT2SBC),1) 5421 end if 5422 5423*---------------------------------------------------------------------* 5424* Compute L^0_3 multipliers: 5425*---------------------------------------------------------------------* 5426 ISYML = ILSTSYM(LISTL,IDLSTL) 5427 5428 CALL DZERO(WORK(KL3AM),NT1AMX*NT1AMX*NT1AMX) 5429 5430 IF (LISTL(1:3).EQ.'L0 ') THEN 5431 ! read left amplitudes from file and square up the doubles part 5432 IF (LWRK1 .LT. NT2AMX) 5433 * CALL QUIT('Insufficient space in CCSDT_GMAT_NODDY') 5434 IOPT = 3 5435 Call CC_RDRSP(LISTL,IDLSTL,ISYML,IOPT,MODEL, 5436 * WORK(KL1AM),WORK(KEND1)) 5437 CALL CC_T2SQ(WORK(KEND1),WORK(KL2AM),ISYML) 5438 5439 IF (NONHF .AND. LWRK1.LT.NT1AMX*NT1AMX*NT1AMX) 5440 * CALL QUIT('Insufficient space in CCSDT_GMAT_NODDY') 5441 5442 CALL CCSDT_L03AM(WORK(KL3AM),WORK(KINT1T0),WORK(KINT2T0), 5443 * WORK(KXIAJB),WORK(KFOCK0),WORK(KL1AM), 5444 * WORK(KL2AM),WORK(KSCR1),WORK(KFOCKD), 5445 * WORK(KFIELD),WORK(KEND1)) 5446 5447 ELSE IF (LISTL(1:3).EQ.'L1 ' .OR. LISTL(1:3).EQ.'LE ' .OR. 5448 & LISTL(1:3).EQ.'M1 ' .OR. LISTL(1:3).EQ.'N2 ' .OR. 5449 & LISTL(1:3).EQ.'E0 ' ) THEN 5450 5451 CALL CCSDT_TBAR31_NODDY(WORK(KL3AM),FREQL,LISTL,IDLSTL, 5452 & WORK(KLAMP0),WORK(KLAMH0), 5453 & WORK(KFOCK0),WORK(KFOCKD),WORK(KSCR1), 5454 & WORK(KXIAJB),WORK(KINT1T0),WORK(KINT2T0), 5455 & WORK(KEND1),LWRK1) 5456 5457 ELSE 5458 5459 CALL QUIT('CCSDT_GMAT_NODDY> LISTL NOT AVAILABLE:'//LISTL) 5460 5461 END IF 5462 5463 CALL DSCAL(NT1AMX*NT1AMX*NT1AMX,-1.0D0,WORK(KL3AM),1) 5464 5465 IF (LOCDBG) THEN 5466 WRITE(LUPRI,*) 'NORM^2(L3AM):', 5467 * DDOT(NT1AMX*NT1AMX*NT1AMX,WORK(KL3AM),1,WORK(KL3AM),1) 5468 END IF 5469 5470*---------------------------------------------------------------------* 5471* Compute contribution from <L_3|[[H^BC,T^0_2],\tau_nu_1]|HF> 5472* and <L_3|[[H^B, T^C_2],\tau_nu_1]|HF> 5473* and <L_3|[[H^C, T^B_2],\tau_nu_1]|HF> 5474*---------------------------------------------------------------------* 5475 KT2AM = KEND1 5476 KEND1 = KT2AM + NT1AMX*NT1AMX 5477 LWRK1 = LWORK - KEND1 5478 IF (LWRK1 .LT. NT2AMX) THEN 5479 CALL QUIT('Insufficient space in CCSDT_GMAT_NODDY') 5480 ENDIF 5481 5482 CALL DZERO(WORK(KOME1),NT1AMX) 5483 5484 ! read T^0 doubles amplitudes from file and square up 5485 IOPT = 2 5486 Call CC_RDRSP('R0',0,ISYM0,IOPT,MODEL,WORK(KDUM),WORK(KEND1)) 5487 CALL CC_T2SQ(WORK(KEND1),WORK(KT2AM),ISYM0) 5488 5489 CALL CC3_L3_OMEGA1_NODDY(WORK(KOME1),WORK(KL3AM), 5490 & WORK(KBCIOOOO),WORK(KBCIOVVO), 5491 & WORK(KBCIOOVV),WORK(KBCIVVVV), 5492 & WORK(KT2AM)) 5493 5494 if (locdbg) then 5495 write(lupri,*) 'norm^2(t2am)-0:', 5496 & ddot(nt1amx*nt1amx,work(kt2am),1,work(kt2am),1) 5497 write(lupri,*) 'norm^2(l3am):', 5498 * ddot(nt1amx*nt1amx*nt1amx,work(kl3am),1,work(kl3am),1) 5499 write(lupri,*) 'norm^2(ome1) after <l3|[[Hbc,T0],tau]|0>:', 5500 & ddot(nt1amx,work(kome1),1,work(kome1),1) 5501 end if 5502 5503 5504 ! read T^C doubles amplitudes from file and square up 5505 IOPT = 2 5506 Call CC_RDRSP(LISTC,IDLSTC,ISYMC,IOPT,MODEL, 5507 & WORK(KDUM),WORK(KEND1)) 5508 Call CCLR_DIASCL(WORK(KEND1),TWO,ISYMC) 5509 CALL CC_T2SQ(WORK(KEND1),WORK(KT2AM),ISYMC) 5510 5511 CALL CC3_L3_OMEGA1_NODDY(WORK(KOME1),WORK(KL3AM), 5512 & WORK(KBIOOOO),WORK(KBIOVVO), 5513 & WORK(KBIOOVV),WORK(KBIVVVV), 5514 & WORK(KT2AM)) 5515 5516 if (locdbg) then 5517 write(lupri,*) 'norm^2(t2am)-C:', 5518 & ddot(nt1amx*nt1amx,work(kt2am),1,work(kt2am),1) 5519 write(lupri,*) 'norm^2(l3am):', 5520 * ddot(nt1amx*nt1amx*nt1amx,work(kl3am),1,work(kl3am),1) 5521 write(lupri,*) 'norm^2(ome1) after <l3|[[Hb,T0],tau]|0>:', 5522 & ddot(nt1amx,work(kome1),1,work(kome1),1) 5523 end if 5524 5525 5526 ! read T^B doubles amplitudes from file and square up 5527 IOPT = 2 5528 Call CC_RDRSP(LISTB,IDLSTB,ISYMB,IOPT,MODEL, 5529 & WORK(KDUM),WORK(KEND1)) 5530 Call CCLR_DIASCL(WORK(KEND1),TWO,ISYMB) 5531 CALL CC_T2SQ(WORK(KEND1),WORK(KT2AM),ISYMB) 5532 5533 CALL CC3_L3_OMEGA1_NODDY(WORK(KOME1),WORK(KL3AM), 5534 & WORK(KCIOOOO),WORK(KCIOVVO), 5535 & WORK(KCIOOVV),WORK(KCIVVVV), 5536 & WORK(KT2AM)) 5537 5538 if (locdbg) then 5539 write(lupri,*) 'norm^2(t2am)-B:', 5540 & ddot(nt1amx*nt1amx,work(kt2am),1,work(kt2am),1) 5541 write(lupri,*) 'norm^2(l3am):', 5542 * ddot(nt1amx*nt1amx*nt1amx,work(kl3am),1,work(kl3am),1) 5543 write(lupri,*) 'norm^2(ome1) after <l3|[[Hc,T0],tau]|0>:', 5544 & ddot(nt1amx,work(kome1),1,work(kome1),1) 5545 end if 5546 5547 DO I = 1,NT1AMX 5548 OMEGA1(I) = OMEGA1(I) + WORK(KOME1+I-1) 5549 END DO 5550 5551*---------------------------------------------------------------------* 5552* Compute contribution from <L_3|[H^BC,\tau_nu_2]|HF> 5553*---------------------------------------------------------------------* 5554 CALL DZERO(WORK(KOME2),NT1AMX*NT1AMX) 5555 5556 CALL CC3_L3_OMEGA2_NODDY(WORK(KOME2),WORK(KL3AM), 5557 * WORK(KINT1SBC),WORK(KINT2SBC)) 5558 DO I = 1,NT1AMX 5559 DO J = 1,I 5560 IJ = NT1AMX*(I-1) + J 5561 NIJ = INDEX(I,J) 5562 OMEGA2(NIJ) = OMEGA2(NIJ) + WORK(KOME2+IJ-1) 5563 END DO 5564 END DO 5565 5566 IF (LOCDBG) THEN 5567 WRITE(LUPRI,*)'G matrix transformation after triples contrib:' 5568 Call CC_PRP(OMEGA1,OMEGA2,1,1,1) 5569 END IF 5570 5571 CALL QEXIT('CCSDT_GMAT_NODDY') 5572 5573 RETURN 5574 END 5575*---------------------------------------------------------------------* 5576* END OF SUBROUTINE CCSDT_GMAT_NODDY * 5577*---------------------------------------------------------------------* 5578