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*=====================================================================* 20 SUBROUTINE CC3_ETASD(LISTL,IDLSTL,FOCKY,ISYFKY, 21 * FREQ,FOCK0, ETA1,ETA2,ETA1EFF, 22 * ETA2EFF,ISYRES,WORK,LWORK,LUDKBC,FNDKBC, 23 * LUCKJD,FNCKJD, 24 * LUTOC,FNTOC,LU3VI,FN3VI,LUDKBC3,FNDKBC3, 25 * LU3FOPX,FN3FOPX,LU3FOP2X,FN3FOP2X) 26*---------------------------------------------------------------------* 27* 28* Purpose: compute triples component of eta vector 29* projected into the singles and doubles space 30* 31* ETA3 = ( <L2|[Y,tau3]|HF> + <L3|[Y^,tau3]|HF>)/ (w+epsiln(tau3)) 32* 33* ETA1EFF(ai) = ETA1(ai) + < ETA3|[[H^,T2],tau(ai)]|HF> 34* 35* ETA2EFF(aibj) = ETA2(aibj) + < ETA3|[H^,tau(aibj)]|HF> 36* 37* 38* To these we add 39* 40* ETA1EFF(ai) = ETA1EFF(ai) + <L3|[[X,tau_ai],T_3]|HF> 41* 42* ETA2EFF(aibj) = ETA2EFF(aibj) + <L3|[[X,E_aiE_bj],T_2]|HF> 43* 44* 45* Written by Poul Jorgensen and Filip Pawlowski, Fall 2002, Aarhus 46* 47*=====================================================================* 48C 49 IMPLICIT NONE 50C 51#include "priunit.h" 52#include "dummy.h" 53#include "iratdef.h" 54#include "ccsdsym.h" 55#include "ccorb.h" 56#include "ccsdinp.h" 57#include "ccinftap.h" 58#include "inftap.h" 59C 60 INTEGER ISYM0 61 PARAMETER(ISYM0 = 1) 62 CHARACTER LISTL*3, MODEL*10 63C 64 CHARACTER*6 FNGEI,FNFEI 65 CHARACTER*5 FNN1 66 PARAMETER( FNGEI = 'N1_GEI' , FNFEI = 'N1_FEI' , FNN1 = 'N1MAT' ) 67 INTEGER LUGEI,LUFEI,LUN1 68C 69 CHARACTER*(*) FNCKJD, FNTOC, FN3VI, FNDKBC3, FN3FOPX, FN3FOP2X 70 CHARACTER*(*) FNDKBC 71 CHARACTER*13 FN4V 72 CHARACTER*11 FNDKBC4 73 CHARACTER*1 CDUMMY 74 LOGICAL LOCDBG 75 INTEGER LU4V 76 INTEGER IDLSTL, LWORK 77 INTEGER ISYFKY, ISYRES 78 INTEGER LUCKJD, LUTOC, LU3VI, LUDKBC3, LU3FOPX, LU3FOP2X 79 INTEGER KLAMP0,KLAMH0,KFOCKD,KFCKBA,KT2TP,KL1AM,KL2TP 80 INTEGER KEND0,LWRK0,KT1AMP0,KEND1,LWRK1 81 INTEGER ISYCTR, ISYMC, ISYMK, KOFF1, KOFF2 82 INTEGER KXIAJB, KTROC01, KTROC21, KTROC03, KTROC23 83 INTEGER KINTOC, KEND2, LWRK2 84 INTEGER IOPT, LENGTH, ISYCKB, ISYMD, ISYMBD 85 INTEGER KTRVI4, KTRVI5, KTRVI7, KEND3, LWRK3 86 INTEGER KTRVI14, KTRVI15, KTRVI18, KTRVI19, KINTVI, KTRVI6 87 INTEGER KEND4, LWRK4, IOFF, ISYMB, ISYALJ, ISYALJ2, ISCKIJ 88 INTEGER ISYCKD, LENSQ 89 INTEGER KSMAT2, KUMAT2, KDIAG, KINDSQ, KINDEX, KINDEX2, KTMAT 90 INTEGER KTRVI16, KTRVI17, KTRVI20, KSMAT4, KUMAT4, KTRVI11 91 INTEGER KTRVI12, KTRVI13, KEND5, LWRK5 92 INTEGER KWMAT, ISWMAT, ISYMIM 93 INTEGER KINDSQW, LENSQW 94 INTEGER KOIOOOO, KOIOVVO, KOIOOVV, KVVVV, ISCKB2, KDIAGW 95 INTEGER LUDKBC4, LUDKBC 96 INTEGER ISINT1, ISINT2, KRBJIA, KTROC, KTROC1, ISCKB1 97 INTEGER KTRVI, KTRVI1 98 INTEGER KRAIJB,KFCKYCK,ISYMT2 99C 100 INTEGER ISYMN1,ISYMN2,KN2MAT,KINDSQN,LENSQN,ISGEI,ISFEI,KGEI,KFEI 101 INTEGER IADR 102C 103 INTEGER kx3am 104C 105 106#if defined (SYS_CRAY) 107 REAL FREQ, HALF, DDOT 108 REAL FOCKY(*), FOCK0(*), WORK(LWORK) 109 REAL ETA1(*), ETA2(*), ETA1EFF(*), ETA2EFF(*) 110 REAL XNORM 111#else 112 DOUBLE PRECISION FREQ, HALF, DDOT 113 DOUBLE PRECISION FOCKY(*), FOCK0(*), WORK(LWORK) 114 DOUBLE PRECISION ETA1(*), ETA2(*), ETA1EFF(*), ETA2EFF(*) 115 DOUBLE PRECISION XNORM 116#endif 117C 118 PARAMETER(HALF = 0.5D0) 119C 120 CALL QENTER('CC3_ETASD') 121C 122C 123C CALL DCOPY(NT1AM(ISYRES),ETA1,1,ETA1EFF,1) 124C CALL DCOPY(NT2AM(ISYRES),ETA2,1,ETA2EFF,1) 125C 126C---------------------------------------------------- 127C Initialise character strings and open files 128C---------------------------------------------------- 129C 130 CDUMMY = ' ' 131 132 LU4V = -1 133C 134 FN4V = 'CC3_ETASD_TMP' 135C 136 IF (LVVVV) THEN 137 CALL WOPEN2(LU4V,FN4V,64,0) 138 END IF 139C 140 LUDKBC4 = -1 141 FNDKBC4 = 'CC3_L3_TMP1' 142C 143 CALL WOPEN2(LUDKBC4,FNDKBC4,64,0) 144C 145 IF (.NOT.LVVVV) THEN 146 !Open files for N1MAT intermediates 147 LUGEI = -1 148 LUFEI = -1 149 LUN1 = -1 150 CALL WOPEN2(LUGEI,FNGEI,64,0) 151 CALL WOPEN2(LUFEI,FNFEI,64,0) 152 CALL WOPEN2(LUN1,FNN1,64,0) 153 END IF 154 155 156 157C 158C------------------------------------------------------------ 159C some initializations: 160C------------------------------------------------------------ 161C 162 163 IF (LISTL(1:3).EQ.'L0 ') THEN 164 ISYCTR = ISYM0 165 ELSE 166 ! ups, probably higher-order response, not yet implemented here 167 CALL QUIT('Unknown type of left vector in CC3_ETASD.') 168 END IF 169 170C 171C------------------------------------------------------- 172C initial allocations, orbital energy, fock matrix and T2 and L2 : 173C------------------------------------------------------- 174C 175C Symmetry of integrals in contraction: 176C 177 ISINT1 = 1 178 ISINT2 = 1 179 ISYMT2 = ISYM0 180 ISYMIM = MULD2H(ISYCTR,ISYFKY) 181 IF (.NOT.LVVVV) THEN 182 !Symmetries for N1 and N2 intermediates 183 ISYMN1 = MULD2H(ISYMIM,ISYMT2) 184 ISYMN2 = MULD2H(ISYMIM,ISYMT2) 185 END IF 186C 187 IF (LVVVV) THEN 188 KRBJIA = 1 189 ELSE 190 KN2MAT = 1 191 KRBJIA = KN2MAT + NCKIJ(ISYMN2) 192 END IF 193C 194 KRAIJB = KRBJIA + NT2SQ(ISYRES) 195 KLAMP0 = KRAIJB + NT2SQ(ISYRES) 196 KLAMH0 = KLAMP0 + NLAMDT 197 KFOCKD = KLAMH0 + NLAMDT 198 KFCKBA = KFOCKD + NORBTS 199 KFCKYCK = KFCKBA + NT1AMX 200 KT2TP = KFCKYCK + NT1AM(ISYFKY) 201 KL1AM = KT2TP + NT2SQ(ISYM0) 202 KL2TP = KL1AM + NT1AM(ISYCTR) 203 KEND0 = KL2TP + NT2SQ(ISYCTR) 204 LWRK0 = LWORK - KEND0 205C 206 IF (.NOT.LVVVV) THEN 207 KINDSQN = KEND0 208 KEND0 = KINDSQN + (6*NCKIJ(ISYMN2) - 1)/IRAT + 1 209 LWRK0 = LWORK - KEND0 210 END IF 211C 212 KT1AMP0 = KEND0 213 KEND1 = KT1AMP0 + NT1AMX 214 LWRK1 = LWORK - KEND1 215 216 CALL DZERO(WORK(KRBJIA),NT2SQ(ISYRES)) 217 CALL DZERO(WORK(KRAIJB),NT2SQ(ISYRES)) 218C 219 IF (.NOT.LVVVV) THEN 220 CALL DZERO(WORK(KN2MAT),NCKIJ(ISYMN2)) 221C 222 !index array for N2 223 LENSQN = NCKIJ(ISYMN2) 224 CALL CC3_INDSQ(WORK(KINDSQN),LENSQN,ISYMN2) 225 END IF 226C 227C------------------------------------- 228C Read in lamdap and lamdh 229C------------------------------------- 230C 231 IOPT = 1 232 CALL CC_RDRSP('R0',0,ISYM0,IOPT,MODEL,WORK(KT1AMP0),DUMMY) 233 234 CALL LAMMAT(WORK(KLAMP0),WORK(KLAMH0),WORK(KT1AMP0), 235 & WORK(KEND1),LWRK1) 236 237C 238C----------------------------------------------------------- 239C Calculate 2*C-E and store 240C FNDKBC3, FN3FOP and FN3FOP2 for f.o.p. later. 241C That's option IOPT = 1 242C 243C With option IOPT = 2 it just calculates LUDKBC4. 244C----------------------------------------------------------- 245C 246 CALL CC3_TCME(WORK(KLAMP0),ISINT1,WORK(KEND1),LWRK1,IDUMMY,CDUMMY, 247 * LUDKBC,FNDKBC,IDUMMY,CDUMMY,IDUMMY,CDUMMY, 248 * IDUMMY,CDUMMY,LUDKBC4,FNDKBC4,2) 249 250 251C 252C------------------------------------- 253C Read T2 amplitudes 254C------------------------------------- 255C 256 IOPT = 2 257 CALL CC_RDRSP('R0',0,ISYM0,IOPT,MODEL,DUMMY,WORK(KT2TP)) 258 259 CALL CC_T2SQ(WORK(KT2TP),WORK(KEND0),ISYM0) 260 261 CALL CC3_T2TP(WORK(KT2TP),WORK(KEND0),ISYM0) 262C 263C IF (LOCDBG) WRITE(LUPRI,*) 'Norm of T2TP ', 264C * DDOT(NT2SQ(ISYM0),WORK(KT2TP),1,WORK(KT2TP),1) 265C 266C------------------------------------- 267C Read L2 amplitudes 268C------------------------------------- 269C 270 IOPT = 3 271 CALL CC_RDRSP(LISTL,IDLSTL,ISYCTR,IOPT,MODEL, 272 * WORK(KL1AM),WORK(KL2TP)) 273 274 CALL CC_T2SQ(WORK(KL2TP),WORK(KEND0),ISYCTR) 275 276 CALL CC3_T2TP(WORK(KL2TP),WORK(KEND0),ISYCTR) 277 278C IF (LOCDBG) WRITE(LUPRI,*) 'Norm of L2TP ', 279C * DDOT(NT2SQ(ISYCTR),WORK(KL2TP),1,WORK(KL2TP),1) 280 281C 282C------------------------------------- 283C Read canonical orbital energies: 284C------------------------------------- 285C 286 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY, 287 & .FALSE.) 288 REWIND LUSIFC 289 290 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 291 READ (LUSIFC) 292 READ (LUSIFC) (WORK(KFOCKD+I-1), I=1,NORBTS) 293 294 CALL GPCLOSE(LUSIFC,'KEEP') 295C 296C--------------------------------------------- 297C Delete frozen orbitals in Fock diagonal. 298C--------------------------------------------- 299C 300 IF (FROIMP .OR. FROEXP) 301 * CALL CCSD_DELFRO(WORK(KFOCKD),WORK(KEND0),LWRK0) 302C 303COMMENT COMMENT COMMENT 304COMMENT COMMENT COMMENT If we want to sum the T3 amplitudes 305COMMENT COMMENT COMMENT 306 if (.false.) then 307 kx3am = kend0 308 kend0 = kx3am + nrhft*nrhft*nrhft*nvirt*nvirt*nvirt 309 call dzero(work(kx3am),nrhft*nrhft*nrhft*nvirt*nvirt*nvirt) 310 lwrk0 = lwork - kend0 311 if (lwrk0 .lt. 0) then 312 write(lupri,*) 'Memory available : ',lwork 313 write(lupri,*) 'Memory needed : ',kend0 314 call quit('Insufficient space (T3) in CC3_etasd') 315 END IF 316 endif 317COMMENT COMMENT COMMENT 318COMMENT COMMENT COMMENT 319COMMENT COMMENT COMMENT 320 321C----------------------------------------------------- 322C Sort the fock matrix 323C----------------------------------------------------- 324C 325 DO ISYMC = 1,NSYM 326 327 ISYMK = MULD2H(ISYMC,ISYM0) 328 329 DO K = 1,NRHF(ISYMK) 330 331 DO C = 1,NVIR(ISYMC) 332 333 KOFF1 = IFCVIR(ISYMK,ISYMC) + 334 * NORB(ISYMK)*(C - 1) + K 335 KOFF2 = IT1AM(ISYMC,ISYMK) 336 * + NVIR(ISYMC)*(K - 1) + C 337 338 WORK(KFCKBA-1+KOFF2) = FOCK0(KOFF1) 339 340 ENDDO 341 ENDDO 342 ENDDO 343 344C IF (LOCDBG) THEN 345C CALL AROUND('In CC3_ETASD: Fock MO matrix (sort)') 346C CALL CC_PRFCKMO(WORK(KFCKBA),ISYM0) 347C ENDIF 348 349C----------------------------------------------------- 350C Sort the FOCKY matrix 351C----------------------------------------------------- 352C 353 DO ISYMC = 1,NSYM 354 355 ISYMK = MULD2H(ISYMC,ISYFKY) 356 357 DO K = 1,NRHF(ISYMK) 358 359 DO C = 1,NVIR(ISYMC) 360 361 KOFF1 = IFCVIR(ISYMK,ISYMC) + 362 * NORB(ISYMK)*(C - 1) + K 363 KOFF2 = IT1AM(ISYMC,ISYMK) 364 * + NVIR(ISYMC)*(K - 1) + C 365 366 WORK(KFCKYCK-1+KOFF2) = FOCKY(KOFF1) 367 368 ENDDO 369 ENDDO 370 ENDDO 371 372* IF (LOCDBG) THEN 373* CALL AROUND('In CC3_ETASD: FOCKY MO matrix (sort)') 374* CALL CC_PRFCKMO(WORK(KFCKYCK),ISYFKY) 375* ENDIF 376 377 378C 379C-------------------------------------------------------------- 380C Calculate the normal g^0 integrals for 381C OOOO, OVVO and OOVV integrals. VVVV is stored on file 382C-------------------------------------------------------------- 383C 384 KOIOOOO = KEND0 385 KOIOVVO = KOIOOOO + N3ORHF(ISYM0) 386 KOIOOVV = KOIOVVO + NT2SQ(ISYM0) 387 KEND0 = KOIOOVV + NT2SQ(ISYM0) 388 LWRK0 = LWORK - KEND0 389C 390 IF (LWRK0 .LT. 0) THEN 391 CALL QUIT('Out of memory in CC3_FMAT (g^0[2o2v] kind)') 392 ENDIF 393C 394 CALL DZERO(WORK(KOIOVVO),NT2SQ(ISYMOP)) 395 CALL DZERO(WORK(KOIOOVV),NT2SQ(ISYMOP)) 396C 397 CALL CC3_INT(WORK(KOIOOOO),WORK(KOIOOVV),WORK(KOIOVVO), 398 * ISYMOP,LU4V,FN4V, 399 * .FALSE.,'DUMMY',IDUMMY,IDUMMY, 400 * .FALSE.,'DUMMY',IDUMMY,IDUMMY, 401 * WORK(KEND0),LWRK0) 402C 403C 404C----------------------------- 405C Memory allocation. 406C----------------------------- 407C 408 KTROC = KEND0 409 KTROC1 = KTROC + NTRAOC(ISINT2) 410 KXIAJB = KTROC1 + NTRAOC(ISINT2) 411 KEND0 = KXIAJB + NT2AM(ISYM0) 412 413 KTROC01 = KEND0 414 KTROC21 = KTROC01 + NTRAOC(ISYM0) 415 KTROC03 = KTROC21 + NTRAOC(ISYM0) 416 KTROC23 = KTROC03 + NTRAOC(ISYM0) 417 KEND0 = KTROC23 + NTRAOC(ISYM0) 418 LWRK0 = LWORK - KEND0 419 420 KINTOC = KEND0 421 KEND2 = KINTOC + MAX(NTOTOC(ISYM0),NTOTOC(ISINT2)) 422 LWRK2 = LWORK - KEND2 423 424 IF (LWRK2 .LT. 0) THEN 425 WRITE(LUPRI,*) 'Memory available : ',LWORK 426 WRITE(LUPRI,*) 'Memory needed : ',KEND2 427 CALL QUIT('Insufficient space in CC3_ETASD') 428 END IF 429C 430C------------------------ 431C Construct L(ia,jb). 432C------------------------ 433C 434 LENGTH = IRAT*NT2AM(ISYM0) 435 436 REWIND(LUIAJB) 437 CALL READI(LUIAJB,LENGTH,WORK(KXIAJB)) 438 439 CALL CCSD_TCMEPK(WORK(KXIAJB),1.0D0,ISYM0,1) 440C 441C------------------------------------------------------------ 442C Read in integrals used in contractions and transform. 443C------------------------------------------------------------ 444C 445 IOFF = 1 446 IF (NTOTOC(ISINT2) .GT. 0) THEN 447 CALL GETWA2(LUCKJD,FNCKJD,WORK(KINTOC),IOFF,NTOTOC(ISINT2)) 448 ENDIF 449C 450 CALL CC3_TROCC(WORK(KINTOC),WORK(KTROC),WORK(KLAMP0), 451 * WORK(KEND2),LWRK2,ISINT2) 452C 453 CALL CCFOP_SORT(WORK(KTROC),WORK(KTROC1),ISINT2,1) 454C 455 CALL CC3_LSORT1(WORK(KTROC),ISINT2,WORK(KEND2),LWRK2,5) 456 457C 458C------------------------ 459C Occupied integrals. 460C------------------------ 461C 462 IOFF = 1 463 IF (NTOTOC(ISYM0) .GT. 0) THEN 464 CALL GETWA2(LUCKJD,FNCKJD,WORK(KINTOC),IOFF,NTOTOC(ISYM0)) 465 ENDIF 466 467* IF (LOCDBG) WRITE(LUPRI,*) 'Norm of OCC-INT ', 468* * DDOT(NTOTOC(ISYM0),WORK(KINTOC),1,WORK(KINTOC),1) 469C 470C----------------------------------------------------------- 471C Construct 2*C-E of the integrals. 472C Have integral for both (ij,k,a) and (a,k,j,i) 473C----------------------------------------------------------- 474C 475 IOFF = 1 476 IF (NTOTOC(ISYM0) .GT. 0) THEN 477 CALL GETWA2(LUTOC,FNTOC,WORK(KINTOC),IOFF,NTOTOC(ISYM0)) 478 ENDIF 479 480 CALL CC3_TROCC(WORK(KINTOC),WORK(KTROC01),WORK(KLAMH0), 481 * WORK(KEND2),LWRK2,ISYM0) 482* IF (LOCDBG) WRITE(LUPRI,*) 'Norm of CKJDEL-INT ', 483* * DDOT(NTOTOC(ISYM0),WORK(KINTOC),1,WORK(KINTOC),1) 484 485 CALL CCSDT_TCMEOCC(WORK(KTROC01),WORK(KTROC21),ISYM0) 486 487 CALL CCFOP_SORT(WORK(KTROC01),WORK(KTROC03),ISYM0,1) 488 489 CALL CCFOP_SORT(WORK(KTROC21),WORK(KTROC23),ISYM0,1) 490 491C 492C----------------------------------------------------- 493C Calculate <L3|[[X,tau_ai],T_3]|HF> 494C----------------------------------------------------- 495C 496 CALL CC3_ETA_1(WORK(KFCKYCK),ISYFKY,ETA1,ISYRES, 497 * WORK(KEND2),LWRK2) 498 499C 500C---------------------------- 501C Loop over D 502C---------------------------- 503C 504 DO ISYMD = 1,NSYM 505 506 ISYCKB = MULD2H(ISYMD,ISYM0) 507 ISCKB1 = MULD2H(ISINT1,ISYMD) 508 ISCKB2 = MULD2H(ISYMD,ISYM0) 509C 510 IF (.NOT.LVVVV) THEN 511 !Symmetry of arrays needed to construct N1MAT 512 ISGEI = MULD2H(ISYMN1,ISYMD) 513 ISFEI = MULD2H(ISYMN1,ISYMD) 514 END IF 515C 516C IF (LOCDBG) WRITE(LUPRI,*) 'In CC3_ETASDPJ: ISYCKB :',ISYCKB 517C 518 KTRVI = KEND0 519 KTRVI1 = KTRVI + NCKATR(ISCKB1) 520 KEND1 = KTRVI1 + NCKATR(ISCKB1) 521 LWRK1 = LWORK - KEND1 522C 523 IF (LVVVV) THEN 524 KVVVV = KEND1 525 KEND1 = KVVVV + NMAABC(ISCKB2) 526 LWRK1 = LWORK - KEND1 527 END IF 528 529C 530 IF (.NOT.LVVVV) THEN 531 !Arrays needed to construct N1MAT 532 KGEI = KEND1 533 KFEI = KGEI + NCKATR(ISGEI) 534 KEND1 = KFEI + NCKATR(ISFEI) 535 LWRK1 = LWORK - KEND1 536 END IF 537C 538 IF (LWRK1 .LT. 0) THEN 539 WRITE(LUPRI,*) 'Memory available : ',LWORK 540 WRITE(LUPRI,*) 'Memory needed : ',KEND1 541 CALL QUIT('Insufficient space in CC3_ETASD') 542 END IF 543C 544 545 546 547 DO D = 1,NVIR(ISYMD) 548C 549 IF (.NOT.LVVVV) THEN 550 CALL DZERO(WORK(KGEI),NCKATR(ISGEI)) 551 CALL DZERO(WORK(KFEI),NCKATR(ISFEI)) 552 END IF 553 554C 555 IF (LVVVV) THEN 556C-------------------------------------------- 557C Read in g^0_{vvvv} for a given D 558C-------------------------------------------- 559C 560 IF (NMAABC(ISCKB2) .GT. 0) THEN 561 IOFF = I3VVIR(ISCKB2,ISYMD) 562 * + NMAABC(ISCKB2)*(D-1) 563 * + 1 564 CALL GETWA2(LU4V,FN4V,WORK(KVVVV),IOFF,NMAABC(ISCKB2)) 565 ENDIF 566C 567 END IF 568 569C 570C------------------------------------------------------------ 571C Read and transform integrals used in contraction. 572C------------------------------------------------------------ 573C 574 IF (NCKATR(ISCKB1) .GT. 0) THEN 575 IOFF = ICKBD(ISCKB1,ISYMD) + NCKATR(ISCKB1)*(D - 1) + 1 576 CALL GETWA2(LUDKBC,FNDKBC,WORK(KTRVI),IOFF, 577 & NCKATR(ISCKB1)) 578 ENDIF 579C 580 IF (LWRK1 .LT. NCKATR(ISCKB1)) THEN 581 CALL QUIT('Insufficient space for allocation in '// 582 & 'CC3_ETASD (TRVI)') 583 END IF 584C 585 CALL CCSDT_SRVIR3(WORK(KTRVI),WORK(KEND1),ISYMD,D,ISINT1) 586C 587C 588 IF (NCKATR(ISCKB1) .GT. 0) THEN 589 IOFF = ICKBD(ISCKB1,ISYMD) + NCKATR(ISCKB1)*(D - 1) + 1 590 CALL GETWA2(LUDKBC4,FNDKBC4,WORK(KTRVI1),IOFF, 591 & NCKATR(ISCKB1)) 592 ENDIF 593C 594 IF (LWRK1 .LT. NCKATR(ISCKB1)) THEN 595 CALL QUIT('Insufficient space for allocation in '// 596 & 'CC3_ETASD (TRVI1)') 597 END IF 598C 599 CALL CCSDT_SRVIR3(WORK(KTRVI1),WORK(KEND1),ISYMD,D,ISINT1) 600C 601 602C 603C 604C ------------------ 605C Memory allocation. 606C ------------------ 607 KTRVI4 = KEND1 608 KTRVI5 = KTRVI4 + NCKATR(ISYCKB) 609 KTRVI7 = KTRVI5 + NCKATR(ISYCKB) 610 KEND3 = KTRVI7 + NCKATR(ISYCKB) 611 LWRK3 = LWORK - KEND3 612 613 KTRVI14 = KEND3 614 KTRVI15 = KTRVI14 + NCKATR(ISYCKB) 615 KTRVI18 = KTRVI15 + NCKATR(ISYCKB) 616 KTRVI19 = KTRVI18 + NCKATR(ISYCKB) 617 KEND3 = KTRVI19 + NCKATR(ISYCKB) 618 LWRK3 = LWORK - KEND3 619 620 621 KINTVI = KEND3 622 KTRVI6 = KINTVI + MAX(NCKA(ISYMD),NCKA(ISYCKB)) 623 KEND4 = KTRVI6 + NCKATR(ISYCKB) 624 LWRK4 = LWORK - KEND4 625 626 IF (LWRK4 .LT. 0) THEN 627 WRITE(LUPRI,*) 'Memory available : ',LWORK 628 WRITE(LUPRI,*) 'Memory needed : ',KEND4 629 CALL QUIT('Insufficient space in CC3_ETASD') 630 END IF 631C 632C ------------------------------------------- 633C Read 2*C-E of integral used for t3-bar 634C ------------------------------------------- 635C 636 IOFF = ICKBD(ISYCKB,ISYMD) + NCKATR(ISYCKB)*(D - 1) + 1 637 IF (NCKATR(ISYCKB) .GT. 0) THEN 638 CALL GETWA2(LU3FOP2X,FN3FOP2X,WORK(KTRVI4),IOFF, 639 & NCKATR(ISYCKB)) 640 ENDIF 641C 642C ------------------------------------------------- 643C Integrals used for t3-bar for cc3 644C ------------------------------------------------- 645C 646 IOFF = ICKBD(ISYCKB,ISYMD) + NCKATR(ISYCKB)*(D - 1) + 1 647 IF (NCKATR(ISYCKB) .GT. 0) THEN 648 CALL GETWA2(LUDKBC3,FNDKBC3,WORK(KTRVI14),IOFF, 649 & NCKATR(ISYCKB)) 650 ENDIF 651 CALL CCSDT_SRVIR3(WORK(KTRVI14),WORK(KEND4), 652 * ISYMD,D,ISYM0) 653 CALL CCSDT_SRTVIR(WORK(KTRVI14),WORK(KTRVI15),WORK(KEND4) 654 * ,LWRK4,ISYMD,ISYM0) 655C 656C ------------------------------------------------ 657C Sort the integrals for t3-bar 658C ------------------------------------------------ 659C 660 661 CALL CCSDT_SRTVIR(WORK(KTRVI4),WORK(KTRVI5),WORK(KEND4), 662 * LWRK4,ISYMD,ISYM0) 663C 664C ---------------------------------------------------- 665C Read virtual integrals used in q3am/u3am for t3-bar. 666C ---------------------------------------------------- 667C 668 IOFF = ICKAD(ISYCKB,ISYMD) + NCKA(ISYCKB)*(D - 1) + 1 669 IF (NCKA(ISYCKB) .GT. 0) THEN 670 CALL GETWA2(LU3VI,FN3VI,WORK(KINTVI),IOFF, 671 * NCKA(ISYCKB)) 672 ENDIF 673 674 CALL CCSDT_TRVIR(WORK(KINTVI),WORK(KTRVI19),WORK(KLAMP0), 675 * ISYMD,D,ISYM0,WORK(KEND4),LWRK4) 676 677 IF (LWRK4 .LT. NCKATR(ISYCKB)) THEN 678 CALL QUIT('Insufficient space for allocation in '// 679 * 'CC3_ETASD (CC3 TRVI)') 680 END IF 681 682 CALL CCSDT_SRTVIR(WORK(KTRVI19),WORK(KTRVI18),WORK(KEND4) 683 * ,LWRK4,ISYMD,ISYM0) 684 685 IOFF = ICKBD(ISYCKB,ISYMD) + NCKATR(ISYCKB)*(D - 1) + 1 686 IF (NCKATR(ISYCKB) .GT. 0) THEN 687 CALL GETWA2(LU3FOPX,FN3FOPX,WORK(KTRVI6),IOFF, 688 * NCKATR(ISYCKB)) 689 ENDIF 690 691 IF (LWRK3 .LT. NCKATR(ISYCKB)) THEN 692 CALL QUIT('Insufficient space for allocation in '// 693 & 'CC3_ETASD (2)') 694 END IF 695 696 CALL DCOPY(NCKATR(ISYCKB),WORK(KTRVI6),1,WORK(KTRVI7),1) 697C 698C 699 DO ISYMB = 1,NSYM 700 701 ISYALJ = MULD2H(ISYMB,ISYM0) 702 ISYALJ2 = MULD2H(ISYMD,ISYM0) 703 ISYMBD = MULD2H(ISYMD,ISYMB) 704 ISCKIJ = MULD2H(ISYMBD,ISYCTR) 705 ISWMAT = MULD2H(ISCKIJ,ISYFKY) 706 ISYCKD = MULD2H(ISYM0,ISYMB) 707 708C Can use kend3 since we do not need the integrals anymore. 709 KSMAT2 = KEND3 710 KUMAT2 = KSMAT2 + NCKIJ(ISCKIJ) 711 KDIAG = KUMAT2 + NCKIJ(ISCKIJ) 712 KDIAGW = KDIAG + NCKIJ(ISCKIJ) 713 KINDSQ = KDIAGW + NCKIJ(ISWMAT) 714 KINDSQW = KINDSQ + (6*NCKIJ(ISCKIJ) - 1)/IRAT + 1 715 KINDEX = KINDSQW + (6*NCKIJ(ISWMAT) - 1)/IRAT + 1 716 KINDEX2 = KINDEX + (NCKI(ISYALJ) - 1)/IRAT + 1 717 KTMAT = KINDEX2 + (NCKI(ISYALJ2) - 1)/IRAT + 1 718 KWMAT = KTMAT + MAX(NCKIJ(ISCKIJ),NCKIJ(ISWMAT)) 719 KEND4 = KWMAT + NCKIJ(ISWMAT) 720 LWRK4 = LWORK - KEND4 721 722 KTRVI16 = KEND4 723 KTRVI17 = KTRVI16 + NCKATR(ISYCKD) 724 KTRVI20 = KTRVI17 + NCKATR(ISYCKD) 725 KEND4 = KTRVI20 + NCKATR(ISYCKD) 726 LWRK4 = LWORK - KEND4 727 728 KSMAT4 = KEND4 729 KUMAT4 = KSMAT4 + NCKIJ(ISCKIJ) 730 KTRVI11 = KUMAT4 + NCKIJ(ISCKIJ) 731 KTRVI12 = KTRVI11 + NCKATR(ISYCKD) 732 KTRVI13 = KTRVI12 + NCKATR(ISYCKD) 733 KEND4 = KTRVI13 + NCKATR(ISYCKD) 734 LWRK4 = LWORK - KEND4 735 736 737 KINTVI = KEND4 738 KEND5 = KINTVI + MAX(NCKA(ISYMB),NCKA(ISYCKD)) 739 LWRK5 = LWORK - KEND5 740 741 IF (LWRK5 .LT. 0) THEN 742 WRITE(LUPRI,*) 'Memory available : ',LWORK 743 WRITE(LUPRI,*) 'Memory needed : ',KEND5 744 CALL QUIT('Insufficient space in CC3_ETASD') 745 END IF 746C 747C 748C ------------------------------- 749C Construct part of the diagonal. 750C ------------------------------- 751C 752 CALL CC3_DIAG(WORK(KDIAG), WORK(KFOCKD),ISCKIJ) 753 CALL CC3_DIAG(WORK(KDIAGW),WORK(KFOCKD),ISWMAT) 754 755C IF (LOCDBG) THEN 756C WRITE(LUPRI,*) 'Norm of DIA ', 757C * DDOT(NCKIJ(ISCKIJ),WORK(KDIAG),1,WORK(KDIAG),1) 758C END IF 759C 760C ----------------------- 761C Construct index arrays. 762C ----------------------- 763C 764 LENSQ = NCKIJ(ISCKIJ) 765 CALL CC3_INDSQ(WORK(KINDSQ),LENSQ,ISCKIJ) 766 LENSQW = NCKIJ(ISWMAT) 767 CALL CC3_INDSQ(WORK(KINDSQW),LENSQW,ISWMAT) 768 CALL CC3_INDEX(WORK(KINDEX),ISYALJ) 769 CALL CC3_INDEX(WORK(KINDEX2),ISYALJ2) 770 771 DO B = 1,NVIR(ISYMB) 772C 773 CALL DZERO(WORK(KWMAT),NCKIJ(ISWMAT)) 774C 775C ---------------------------- 776C Read and transform integrals 777C ---------------------------- 778 IOFF = ICKBD(ISYCKD,ISYMB) 779 * + NCKATR(ISYCKD)*(B - 1) + 1 780 IF (NCKATR(ISYCKD) .GT. 0) THEN 781 CALL GETWA2(LUDKBC3,FNDKBC3,WORK(KTRVI16),IOFF, 782 * NCKATR(ISYCKD)) 783 ENDIF 784 CALL CCSDT_SRVIR3(WORK(KTRVI16),WORK(KEND5), 785 * ISYMB,B,ISYM0) 786 CALL CCSDT_SRTVIR(WORK(KTRVI16),WORK(KTRVI17), 787 * WORK(KEND4),LWRK4,ISYMB,ISYM0) 788C 789C ------------------------------------ 790C Read and transform integrals used in 791C second S-bar and U-bar 792C ------------------------------------ 793 IOFF = ICKBD(ISYCKD,ISYMB) 794 * + NCKATR(ISYCKD)*(B-1) + 1 795 IF (NCKATR(ISYCKD) .GT. 0) THEN 796 CALL GETWA2(LU3FOP2X,FN3FOP2X,WORK(KTRVI11), 797 * IOFF,NCKATR(ISYCKD)) 798 ENDIF 799 800 CALL CCSDT_SRTVIR(WORK(KTRVI11),WORK(KTRVI12), 801 * WORK(KEND5),LWRK5,ISYMB, 802 * ISYM0) 803 804 IOFF = ICKBD(ISYCKD,ISYMB) 805 * + NCKATR(ISYCKD)*(B - 1) + 1 806 IF (NCKATR(ISYCKD) .GT. 0) THEN 807 CALL GETWA2(LU3FOPX,FN3FOPX,WORK(KTRVI13),IOFF, 808 * NCKATR(ISYCKD)) 809 ENDIF 810 811 IOFF = ICKAD(ISYCKD,ISYMB) 812 * + NCKA(ISYCKD)*(B - 1) + 1 813 IF (NCKA(ISYCKD) .GT. 0) THEN 814 CALL GETWA2(LU3VI,FN3VI,WORK(KINTVI),IOFF, 815 * NCKA(ISYCKD)) 816 ENDIF 817 818 CALL CCSDT_TRVIR(WORK(KINTVI),WORK(KTRVI20), 819 * WORK(KLAMP0),ISYMB,B,ISYM0, 820 * WORK(KEND4),LWRK4) 821C 822C ---------------------------------------------------- 823C Calculate the S(ci,bk,dj) matrix for B,D for T3-BAR: 824C ---------------------------------------------------- 825 CALL DZERO(WORK(KSMAT2),NCKIJ(ISCKIJ)) 826 827 CALL CCFOP_SMAT(0.0D0,WORK(KL1AM),ISYCTR,WORK(KL2TP), 828 * ISYCTR,WORK(KTMAT), 829 * WORK(KFCKBA),WORK(KXIAJB),ISYM0, 830 * WORK(KTRVI14),WORK(KTRVI15), 831 * WORK(KTRVI4),WORK(KTRVI5), 832 * WORK(KTROC01),WORK(KTROC21), 833 * ISYM0,WORK(KFOCKD),WORK(KDIAG), 834 * WORK(KSMAT2),WORK(KEND4),LWRK4, 835 * WORK(KINDEX),WORK(KINDSQ),LENSQ, 836 * ISYMB,B,ISYMD,D) 837 838 CALL DSCAL(NCKIJ(ISCKIJ),HALF,WORK(KSMAT2),1) 839 840 CALL T3_FORBIDDEN(WORK(KSMAT2),ISYCTR,ISYMB,B,ISYMD,D) 841C 842C ---------------------------------------------------- 843C Calculate the S(ci,bk,dj) matrix for D,B for T3-BAR: 844C ---------------------------------------------------- 845 CALL DZERO(WORK(KSMAT4),NCKIJ(ISCKIJ)) 846 847 CALL CCFOP_SMAT(0.0D0,WORK(KL1AM),ISYCTR,WORK(KL2TP), 848 * ISYCTR,WORK(KTMAT),WORK(KFCKBA), 849 * WORK(KXIAJB),ISYM0, 850 * WORK(KTRVI16),WORK(KTRVI17), 851 * WORK(KTRVI11),WORK(KTRVI12), 852 * WORK(KTROC01),WORK(KTROC21), 853 * ISYM0,WORK(KFOCKD),WORK(KDIAG), 854 * WORK(KSMAT4),WORK(KEND4),LWRK4, 855 * WORK(KINDEX2),WORK(KINDSQ), 856 * LENSQ,ISYMD,D,ISYMB,B) 857 858 CALL DSCAL(NCKIJ(ISCKIJ),HALF,WORK(KSMAT4),1) 859C 860C 861 CALL T3_FORBIDDEN(WORK(KSMAT4),ISYCTR,ISYMD,D,ISYMB,B) 862C 863C ------------------------------------------------ 864C Calculate U(ci,jk) for fixed b,d for t3-bar. 865C ------------------------------------------------ 866 CALL DZERO(WORK(KUMAT2),NCKIJ(ISCKIJ)) 867 868 CALL CCFOP_UMAT(0.0D0,WORK(KL1AM),ISYCTR,WORK(KL2TP), 869 * ISYCTR, 870 * WORK(KXIAJB),ISYM0,WORK(KFCKBA), 871 * WORK(KTRVI19),WORK(KTRVI7), 872 * WORK(KTROC03),WORK(KTROC23),ISYM0, 873 * WORK(KFOCKD),WORK(KDIAG), 874 * WORK(KUMAT2), 875 * WORK(KTMAT),WORK(KEND4),LWRK4, 876 * WORK(KINDSQ),LENSQ,ISYMB,B,ISYMD,D) 877 878 CALL DSCAL(NCKIJ(ISCKIJ),HALF,WORK(KUMAT2),1) 879C 880C 881 CALL T3_FORBIDDEN(WORK(KUMAT2),ISYCTR,ISYMB,B,ISYMD,D) 882C 883C ------------------------------------------------ 884C Calculate U(ci,jk) for fixed d,b for t3-bar. 885C ------------------------------------------------ 886 CALL DZERO(WORK(KUMAT4),NCKIJ(ISCKIJ)) 887 888 CALL CCFOP_UMAT(0.0D0,WORK(KL1AM),ISYCTR,WORK(KL2TP), 889 * ISYCTR,WORK(KXIAJB),ISYM0, 890 * WORK(KFCKBA),WORK(KTRVI20), 891 * WORK(KTRVI13),WORK(KTROC03), 892 * WORK(KTROC23),ISYM0, 893 * WORK(KFOCKD),WORK(KDIAG), 894 * WORK(KUMAT4),WORK(KTMAT), 895 * WORK(KEND4),LWRK4,WORK(KINDSQ), 896 * LENSQ,ISYMD,D,ISYMB,B) 897 898 CALL DSCAL(NCKIJ(ISCKIJ),HALF,WORK(KUMAT4),1) 899 900 CALL T3_FORBIDDEN(WORK(KUMAT4),ISYCTR,ISYMD,D,ISYMB,B) 901C 902C ------------------------------------------- 903C Sum up the S-bar and U-bar to get a real T3-bar 904C ------------------------------------------- 905 CALL CC3_T3BD(ISCKIJ,WORK(KSMAT2),WORK(KSMAT4), 906 * WORK(KUMAT2),WORK(KUMAT4), 907 * WORK(KTMAT),WORK(KINDSQ),LENSQ) 908C 909C ------------------------------------- 910C 911C---------------------------------------------------------- 912C Calculate <L3|[[X,E_aiE_bj],T_2]|HF> 913C---------------------------------------------------------- 914C 915* ISYMT2 = ISYM0 916C 917 CALL CC3_ETA_2(WORK(KTMAT),ISCKIJ,WORK(KFCKYCK), 918 * ISYFKY,WORK(KT2TP),ISYMT2,ETA2, 919 * ISYRES,WORK(KRAIJB),WORK(KINDSQ), 920 * LENSQ,WORK(KINDEX2),ISYALJ2, 921 * ISYMB,B,ISYMD,D,WORK(KEND4),LWRK4) 922C 923C calculate <L3|[Y^,tau3]|HF> 924C 925 CALL WBARBD_V(WORK(KTMAT),ISCKIJ,FOCKY,ISYFKY, 926 * WORK(KWMAT),ISWMAT,WORK(KEND4),LWRK4) 927C 928 CALL WBARBD_O(WORK(KTMAT),ISCKIJ,FOCKY,ISYFKY, 929 * WORK(KWMAT),ISWMAT,WORK(KEND4),LWRK4) 930 931C calculate <L2|[Y,tau3]|HF> 932C 933 CALL WBARBD_T2(B,ISYMB,D,ISYMD,WORK(KL2TP),ISYCTR, 934 * FOCKY,ISYFKY,WORK(KWMAT),ISWMAT) 935C 936COMMENT COMMENT COMMENT 937C call sum_pt3(work(kwmat),isymb,b,isymd,d, 938C * iswmat,work(kx3am),4) 939COMMENT COMMENT COMMENT 940C 941C 942C------------------------------------------------ 943C Divide by the energy difference and 944C remove the forbidden elements 945C------------------------------------------------ 946C 947 CALL WBD_DIA(B,ISYMB,D,ISYMD,-FREQ, 948 * ISWMAT,WORK(KWMAT), 949 * WORK(KDIAGW),WORK(KFOCKD)) 950C 951 CALL T3_FORBIDDEN(WORK(KWMAT),ISYFKY,ISYMB,B,ISYMD,D) 952C 953 954COMMENT COMMENT COMMENT 955COMMENT COMMENT COMMENT 956C call sum_pt3(work(kwmat),isymb,b,isymd,d, 957C * iswmat,work(kx3am),4) 958COMMENT COMMENT COMMENT 959COMMENT COMMENT COMMENT 960C 961* ISYMIM = MULD2H(ISWMAT,ISYMBD) 962C 963C----------------------------------------- 964C Calculate <E_3|[H,\tau_nu_2]|HF> 965C----------------------------------------- 966C 967 CALL CC3_W3_CY2V(ETA2EFF,ISYRES,WORK(KRBJIA), 968 * WORK(KWMAT),ISWMAT, 969 * WORK(KTMAT),WORK(KTRVI),WORK(KTRVI1), 970 * ISINT1,WORK(KEND4),LWRK4, 971 * WORK(KINDSQW),LENSQW, 972 * ISYMB,B,ISYMD,D,.TRUE.) 973 974 CALL CC3_W3_CY2O(ETA2EFF,ISYRES,WORK(KWMAT),ISWMAT, 975 * WORK(KTMAT),WORK(KTROC),WORK(KTROC1), 976 * ISINT1,WORK(KEND4),LWRK4, 977 * WORK(KINDSQW),LENSQW,ISYMB,B,ISYMD,D, 978 * .TRUE.) 979 980C 981C----------------------------------------- 982C Calculate <E_3|[[H,T^0_2],\tau_nu_1]|HF> 983C----------------------------------------- 984C 985 IF (LVVVV) THEN 986 CALL CC3_W3_OMEGA1(ETA1EFF,ISYRES,WORK(KWMAT), 987 * WORK(KTMAT),ISYMIM, 988 * WORK(KOIOOOO),WORK(KOIOVVO), 989 * WORK(KOIOOVV),WORK(KVVVV),ISYM0, 990 * WORK(KT2TP),ISYM0, 991 * WORK(KEND4),LWRK4, 992 * LENSQW,WORK(KINDSQW), 993 * ISYMB,B,ISYMD,D,.TRUE.) 994 ELSE 995 CALL DSCAL(NCKIJ(ISWMAT),-1.0D0,WORK(KWMAT),1) 996 997 !Construct N1 and N2 intermediates 998 CALL WT2_N1N2(WORK(KWMAT),ISYMIM, 999 * WORK(KT2TP),ISYM0, 1000 * WORK(KGEI),WORK(KFEI), 1001 * ISYMN1, 1002 * WORK(KN2MAT),ISYMN2, 1003 * B,ISYMB,D,ISYMD, 1004 * WORK(KINDSQW),LENSQW, 1005 * WORK(KINDSQN),LENSQN, 1006 * WORK(KEND4),LWRK4, 1007 * .TRUE.) 1008 1009 END IF 1010 1011 ENDDO ! B 1012 ENDDO ! ISYMB 1013 1014 IF (.NOT.LVVVV) THEN 1015C 1016C ---------------------------------------------------------- 1017C Put KGEI(ge,i)^F and KFEI(fe,i)^G (which are intermediates 1018C for N1MAT(fge,i) ) to files (for fixed F=D and G=D). 1019C ---------------------------------------------------------- 1020 1021 !Put KGEI to file as (gei,F) (fixed F corresponds to D) 1022 IADR = ICKBD(ISGEI,ISYMD) + NCKATR(ISGEI)*(D-1) + 1 1023 CALL PUTWA2(LUGEI,FNGEI,WORK(KGEI),IADR,NCKATR(ISGEI)) 1024C 1025 !Put KFEI to file as (fei,G) (fixed G corresponds to D) 1026 IADR = ICKBD(ISFEI,ISYMD) + NCKATR(ISFEI)*(D-1) + 1 1027 CALL PUTWA2(LUFEI,FNFEI,WORK(KFEI),IADR,NCKATR(ISFEI)) 1028C 1029 END IF 1030 1031 1032 ENDDO ! D 1033 ENDDO ! ISYMD 1034C------------------------------------------------------ 1035C Accumulate RBJIA from <mu2|[H,W^BD(3)]|HF> ( Vccupied cont ) 1036C in ETA2EFF 1037C 1038C Accumulate RAIJB from <L3|[[X,E_aiE_bj],T_2]|HF> 1039C in ETA2EFF 1040C------------------------------------------------------ 1041C 1042 CALL CC3_RBJIA(ETA2EFF,ISYRES,WORK(KRBJIA)) 1043 CALL CC3_RBJIA(ETA2EFF,ISYRES,WORK(KRAIJB)) 1044 1045C 1046 IF (.NOT.LVVVV) THEN 1047C 1048 !Read (gei,F) and (fei,G) intermediates from files 1049 !add them and put the result to a file as (fge,I) 1050 CALL N1_RESORT(ISYMN1,LUN1,FNN1,LUGEI,FNGEI,LUFEI,FNFEI, 1051 * WORK(KEND0),LWRK0,.FALSE.) 1052C 1053 !Calculate <T3|[[H,T2],tau_ai]|HF> except VVVV contribution 1054 CALL N1N2_G(LUN1,FNN1, 1055 * ISYMN1, 1056 * WORK(KN2MAT),ISYMN2, 1057 * WORK(KOIOVVO),WORK(KOIOOVV), 1058 * WORK(KOIOOOO),ISYM0, 1059 * ETA1EFF,ISYRES, 1060 * WORK(KINDSQN),LENSQN, 1061 * WORK(KEND0),LWRK0) 1062C 1063 !Calculate VVVV contribution to <T3|[[H,T2],tau_ai]|HF> 1064 IOPT = 0 !normal Lambda matrices used in backtransformation 1065 CALL N1_GV4(IOPT, 1066 * LUN1,FNN1, 1067 * ISYMN1, 1068 * WORK(KLAMP0),1, 1069 * WORK(KLAMP0),1, 1070 * WORK(KLAMH0),1, 1071 * WORK(KLAMH0),1, 1072 * ETA1EFF,ISYRES, 1073 * WORK(KEND0),LWRK0) 1074C 1075 END IF 1076C 1077COMMENT COMMENT 1078C write(lupri,*) 'WMAT in CC3_ETASD : ' 1079C call print_pt3(work(kx3am),ISYFKY,4) 1080COMMENT COMMENT 1081C 1082C------------------------------- 1083C Close and delete files 1084C------------------------------- 1085C 1086 IF (LVVVV) THEN 1087 CALL WCLOSE2(LU4V,FN4V,'DELETE') 1088 END IF 1089C 1090 CALL WCLOSE2(LUDKBC4,FNDKBC4,'DELETE') 1091C 1092 IF (.NOT.LVVVV) THEN 1093 !Close files for N1MAT intermediates 1094 CALL WCLOSE2(LUGEI,FNGEI,'DELETE') 1095 CALL WCLOSE2(LUFEI,FNFEI,'DELETE') 1096 CALL WCLOSE2(LUN1,FNN1,'DELETE') 1097 END IF 1098C 1099C------------------------------------------ 1100C Accumulate CCSD and CC3 contributions 1101C------------------------------------------ 1102C 1103 DO I = 1,NT2AM(ISYRES) 1104 ETA2EFF(I) = ETA2EFF(I) + ETA2(I) 1105 END DO 1106C 1107 DO I = 1,NT1AM(ISYRES) 1108 ETA1EFF(I) = ETA1EFF(I) + ETA1(I) 1109 END DO 1110C 1111 1112C 1113C------------- 1114C End 1115C------------- 1116C 1117C 1118 CALL QEXIT('CC3_ETASD') 1119C 1120 RETURN 1121 END 1122 1123C /* Deck wbarbd_v */ 1124 SUBROUTINE WBARBD_V(TMAT,ISTMAT,FOCKY,ISYFKY, 1125 * WMAT,ISWMAT,WRK,LWRK) 1126C 1127C WBD(a,i,k,j) = WBD(a,i,k,j) + sum (f) focky(f,a)*tmatBD(f,i,k,j) 1128C 1129C 1130C Written by P. Jorgensen and F. Pawlowski, Spring 2002. 1131C 1132 1133C 1134 IMPLICIT NONE 1135C 1136 INTEGER LWRK, KFCAF, KEND0, LWRK0, KOFF1, KOFF2 1137 INTEGER NF, KOFFY, KOFFT, KOFFW 1138 INTEGER ISTMAT, ISYFKY, ISWMAT, ISFIKJ, ISYFIK 1139 INTEGER ISYMA, ISYAI, ISYAIK, NA 1140 INTEGER ISYMF, ISYMJ, ISYMK, ISYMI, ISYFI 1141#if defined (SYS_CRAY) 1142 REAL TMAT(*), FOCKY(*), WMAT(*), WRK(*) 1143 REAL HALF, ONE 1144#else 1145 DOUBLE PRECISION TMAT(*), FOCKY(*), WMAT(*), WRK(*) 1146 DOUBLE PRECISION HALF, ONE 1147#endif 1148C 1149#include "priunit.h" 1150#include "dummy.h" 1151#include "iratdef.h" 1152#include "ccsdsym.h" 1153#include "inftap.h" 1154#include "ccinftap.h" 1155#include "ccorb.h" 1156#include "ccsdinp.h" 1157C 1158 PARAMETER (ONE = 1.0D0) 1159C 1160 CALL QENTER('WBARBD_V') 1161C 1162C 1163C RESORT VIR-VIR FOCKY ELEMENTS (A,F) 1164C 1165C 1166 KFCAF = 1 1167 KEND0 = KFCAF + NMATAB(ISYFKY) 1168 LWRK0 = LWRK - KEND0 1169C 1170 IF (LWRK0 .LT. 0) THEN 1171 WRITE(LUPRI,*) 'Memory available : ',LWRK0 1172 WRITE(LUPRI,*) 'Memory needed : ',KEND0 1173 CALL QUIT('Insufficient space in WBARBD_V') 1174 END IF 1175C 1176 DO ISYMF = 1,NSYM 1177 ISYMA = MULD2H(ISYMF,ISYFKY) 1178 DO F = 1,NVIR(ISYMF) 1179 KOFF1 = IFCVIR(ISYMA,ISYMF) + NORB(ISYMA)*(F - 1) 1180 * + NRHF(ISYMA) + 1 1181 KOFF2 = KFCAF + IMATAB(ISYMA,ISYMF) + NVIR(ISYMA)*(F - 1) 1182 CALL DCOPY(NVIR(ISYMA),FOCKY(KOFF1),1,WRK(KOFF2),1) 1183 END DO 1184 END DO 1185C 1186C CARRY OUT MATRIX MULTIPLICATION 1187C WBD(a,i,k,j) = WBD(a,i,k,j) + sum (f) focky(f,a)*tmatBD(f,i,k,j) 1188C 1189 ISFIKJ = ISTMAT 1190 DO ISYMJ = 1,NSYM 1191 ISYFIK =MULD2H(ISYMJ,ISFIKJ) 1192 DO J = 1,NRHF(ISYMJ) 1193 DO ISYMK = 1,NSYM 1194 ISYFI = MULD2H(ISYMK,ISYFIK) 1195 DO K = 1,NRHF(ISYMK) 1196 DO ISYMI = 1,NSYM 1197 ISYMF = MULD2H(ISYFI,ISYMI) 1198 ISYMA = MULD2H(ISYFKY,ISYMF) 1199 ISYAIK = MULD2H(ISWMAT,ISYMJ) 1200 ISYAI = MULD2H(ISYAIK,ISYMK) 1201 NA = MAX(1,NVIR(ISYMA)) 1202 NF = MAX(1,NVIR(ISYMF)) 1203 KOFFY = KFCAF + IMATAB(ISYMF,ISYMA) 1204 KOFFT = ISAIKJ(ISYFIK,ISYMJ) 1205 * + NCKI(ISYFIK)*(J-1) 1206 * + ISAIK(ISYFI,ISYMK) 1207 * + NT1AM(ISYFI)*(K-1) 1208 * + IT1AM(ISYMF,ISYMI) + 1 1209 KOFFW = ISAIKJ(ISYAIK,ISYMJ) 1210 * + NCKI(ISYAIK)*(J-1) 1211 * + ISAIK(ISYAI,ISYMK) 1212 * + NT1AM(ISYAI)*(K-1) 1213 * + IT1AM(ISYMA,ISYMI) + 1 1214C 1215C SYMMETRY BETWEEN BJ AND CK INTRODUCE A FACTOR TWO 1216C DENOTE t3 IS CALCULATED WITH NEGATIVE SIGN 1217C 1218C WBD(a,i,k,j) = WBD(a,i,k,j) + sum (f) focky(f,a)*tmatBD(f,i,k,j) 1219C 1220 CALL DGEMM('T','N',NVIR(ISYMA),NRHF(ISYMI), 1221 * NVIR(ISYMF),-ONE,WRK(KOFFY),NF, 1222 * TMAT(KOFFT),NF,ONE,WMAT(KOFFW),NA) 1223C 1224 END DO 1225 END DO 1226 END DO 1227 END DO 1228 END DO 1229C 1230 CALL QEXIT('WBARBD_V') 1231C 1232 RETURN 1233 END 1234C 1235C /* Deck wbarbd_o */ 1236 SUBROUTINE WBARBD_O(TMAT,ISTMAT,FOCKY,ISYFKY, 1237 * WMAT,ISWMAT,WRK,LWRK) 1238C 1239C WBD(a,i,k,j) = WBD(a,i,k,j) - sum (l) tmatBD(a,l,k,j)*focky(i,l) 1240C 1241C 1242C Written by P. Jorgensen and F. Pawlowski, Spring 2002. 1243C 1244 1245 IMPLICIT NONE 1246C 1247 INTEGER LWRK, KFCLI, KEND0, LWRK0, KOFF1, KOFF2 1248 INTEGER NI, KOFFY, KOFFT, KOFFW 1249 INTEGER ISTMAT, ISYFKY, ISWMAT, ISALKJ 1250 INTEGER ISYMA, ISYAI, ISYAIK, ISYALK, ISYAL, NA 1251 INTEGER ISYMJ, ISYMK, ISYMI, ISYML, ISYFI 1252C 1253#if defined (SYS_CRAY) 1254 REAL TMAT(*), FOCKY(*), WMAT(*), WRK(*) 1255 REAL MHALF, ONE 1256#else 1257 DOUBLE PRECISION TMAT(*), FOCKY(*), WMAT(*), WRK(*) 1258 DOUBLE PRECISION MHALF, ONE 1259#endif 1260C 1261#include "priunit.h" 1262#include "dummy.h" 1263#include "iratdef.h" 1264#include "ccsdsym.h" 1265#include "inftap.h" 1266#include "ccinftap.h" 1267#include "ccorb.h" 1268#include "ccsdinp.h" 1269C 1270 PARAMETER ( ONE = 1.0D0) 1271C 1272 CALL QENTER('WBARBD_O') 1273C 1274 1275C 1276C RESORT OCC-OCC FOCKY ELEMENTS (L,I) 1277C 1278C 1279 KFCLI = 1 1280 KEND0 = KFCLI + NMATIJ(ISYFKY) 1281 LWRK0 = LWRK - KEND0 1282C 1283 IF (LWRK0 .LT. 0) THEN 1284 WRITE(LUPRI,*) 'Memory available : ',LWRK0 1285 WRITE(LUPRI,*) 'Memory needed : ',KEND0 1286 CALL QUIT('Insufficient space in WBARBD_O') 1287 END IF 1288C 1289 DO ISYMI = 1,NSYM 1290 ISYML = MULD2H(ISYMI,ISYFKY) 1291 DO I = 1,NRHF(ISYMI) 1292 KOFF1 = IFCRHF(ISYML,ISYMI) + NORB(ISYML)*(I - 1) + 1 1293 KOFF2 = KFCLI + IMATIJ(ISYML,ISYMI) + NRHF(ISYML)*(I - 1) 1294 CALL DCOPY(NRHF(ISYML),FOCKY(KOFF1),1,WRK(KOFF2),1) 1295 END DO 1296 END DO 1297C 1298C CARRY OUT MATRIX MULTIPLICATION 1299C WBD(a,i,k,j) = WBD(a,i,k,j) - sum (l) tmatBD(a,l,k,j)*focky(i,l) 1300C 1301 ISALKJ = ISTMAT 1302 DO ISYMJ = 1,NSYM 1303 ISYALK =MULD2H(ISYMJ,ISALKJ) 1304 DO J = 1,NRHF(ISYMJ) 1305 DO ISYMK = 1,NSYM 1306 ISYAL = MULD2H(ISYMK,ISYALK) 1307 DO K = 1,NRHF(ISYMK) 1308 DO ISYML = 1,NSYM 1309 ISYMA = MULD2H(ISYAL,ISYML) 1310 ISYMI = MULD2H(ISYFKY,ISYML) 1311 ISYAIK = MULD2H(ISWMAT,ISYMJ) 1312 ISYAI = MULD2H(ISYAIK,ISYMK) 1313 NA = MAX(1,NVIR(ISYMA)) 1314 NI = MAX(1,NRHF(ISYMI)) 1315 KOFFY = KFCLI + IMATIJ(ISYMI,ISYML) 1316 KOFFT = ISAIKJ(ISYALK,ISYMJ) 1317 * + NCKI(ISYALK)*(J-1) 1318 * + ISAIK(ISYAL,ISYMK) 1319 * + NT1AM(ISYAL)*(K-1) 1320 * + IT1AM(ISYMA,ISYML) + 1 1321 KOFFW = ISAIKJ(ISYAIK,ISYMJ) 1322 * + NCKI(ISYAIK)*(J-1) 1323 * + ISAIK(ISYAI,ISYMK) 1324 * + NT1AM(ISYAI)*(K-1) 1325 * + IT1AM(ISYMA,ISYMI) + 1 1326C 1327C SYMMETRY BETWEEN BJ AND CK INTRODUCE A FACTOR TWO 1328C DENOTE t3 IS CALCULATED WITH NEGATIVE SIGN 1329C 1330 CALL DGEMM('N','T',NVIR(ISYMA),NRHF(ISYMI), 1331 * NRHF(ISYML),ONE,TMAT(KOFFT),NA, 1332 * WRK(KOFFY),NI,ONE,WMAT(KOFFW),NA) 1333 END DO 1334 END DO 1335 END DO 1336 END DO 1337 END DO 1338C 1339 CALL QEXIT('WBARBD_O') 1340C 1341 RETURN 1342 END 1343C 1344C /* Deck wbarbd_t2 */ 1345 SUBROUTINE WBARBD_T2(B,ISYMB,D,ISYMD,T2TP,ISYMT2,FOCKY,ISYFKY, 1346 * WMAT,ISWMAT) 1347C 1348C WBD(a,i,k,j) = WBD(a,i,k,j) + 1349C 1350C focky(j,B)*t2(ai,Dk) - focky(k,B)*t2(ai,Dj) 1351C focky(k,D)*t2(ai,Bj) - focky(j,D)*t2(ai,Bk) 1352C 1353C Written by P. Jorgensen and F. Pawlowski, Spring 2002. 1354C 1355 1356 IMPLICIT NONE 1357C 1358 INTEGER ISYMB, ISYMD, ISYMT2, ISYFKY, ISWMAT 1359 INTEGER ISYMJ, KJB, KJD, ISYMK, KKB, KKD, ISYMI, ISYIJ, ISYIK 1360 INTEGER ISYMA, ISYAI, ISYAIK, ISYAIJ, KAIKD, KAIJD, KAIJB 1361 INTEGER KAIKB, KAIKJ 1362C 1363#if defined (SYS_CRAY) 1364 REAL T2TP(*), FOCKY(*), WMAT(*) 1365#else 1366 DOUBLE PRECISION T2TP(*), FOCKY(*), WMAT(*) 1367#endif 1368C 1369#include "priunit.h" 1370#include "ccsdsym.h" 1371#include "ccorb.h" 1372#include "ccsdinp.h" 1373C 1374 CALL QENTER('WBARBD_T2') 1375C 1376C 1377C focky(j,B)*t2(ai,Dk) - focky(k,B)*t2(ai,Dj) 1378C focky(k,D)*t2(ai,Bj) - focky(j,D)*t2(ai,Bk) 1379C 1380 1381C 1382C (1) wmat(aikj) = wmat(aikj) + focky(j,B)*t2(ai,Dk) 1383C 1384 ISYMJ = MULD2H(ISYFKY,ISYMB) 1385 ISYAIK = MULD2H(ISYMT2,ISYMD) 1386 DO ISYMK = 1,NSYM 1387 ISYAI = MULD2H(ISYAIK,ISYMK) 1388 DO ISYMI = 1,NSYM 1389 ISYMA = MULD2H(ISYAI,ISYMI) 1390 DO J = 1,NRHF(ISYMJ) 1391 KJB = IFCVIR(ISYMJ,ISYMB) + NORB(ISYMJ)*(B - 1) + J 1392 DO K = 1,NRHF(ISYMK) 1393 DO I = 1,NRHF(ISYMI) 1394 DO A = 1,NVIR(ISYMA) 1395 KAIKD = IT2SP(ISYAIK,ISYMD) 1396 * + NCKI(ISYAIK)*(D-1) 1397 * + ISAIK(ISYAI,ISYMK) 1398 * + NT1AM(ISYAI)*(K-1) 1399 * + IT1AM(ISYMA,ISYMI) 1400 * + NVIR(ISYMA)*(I-1) 1401 * + A 1402 1403 KAIKJ = ISAIKJ(ISYAIK,ISYMJ) 1404 * + NCKI(ISYAIK)*(J-1) 1405 * + ISAIK(ISYAI,ISYMK) 1406 * + NT1AM(ISYAI)*(K-1) 1407 * + IT1AM(ISYMA,ISYMI) 1408 * + NVIR(ISYMA)*(I-1) 1409 * + A 1410 1411 WMAT(KAIKJ) = WMAT(KAIKJ) 1412 * + FOCKY(KJB)*T2TP(KAIKD) 1413 END DO 1414 END DO 1415 END DO 1416 END DO 1417 END DO 1418 END DO 1419 1420C 1421C (2) wmat(aikj) = wmat(aikj) - focky(k,B)*t2(ai,Dj) 1422C 1423 ISYMK = MULD2H(ISYFKY,ISYMB) 1424 ISYAIJ = MULD2H(ISYMT2,ISYMD) 1425 DO ISYMJ = 1,NSYM 1426 ISYAI = MULD2H(ISYAIJ,ISYMJ) 1427 ISYAIK = MULD2H(ISYAI,ISYMK) 1428 DO ISYMI = 1,NSYM 1429 ISYMA = MULD2H(ISYAI,ISYMI) 1430 DO J = 1,NRHF(ISYMJ) 1431 DO K = 1,NRHF(ISYMK) 1432 KKB = IFCVIR(ISYMK,ISYMB) + NORB(ISYMK)*(B - 1) + K 1433 DO I = 1,NRHF(ISYMI) 1434 DO A = 1,NVIR(ISYMA) 1435 1436 KAIJD = IT2SP(ISYAIJ,ISYMD) 1437 * + NCKI(ISYAIJ)*(D-1) 1438 * + ISAIK(ISYAI,ISYMJ) 1439 * + NT1AM(ISYAI)*(J-1) 1440 * + IT1AM(ISYMA,ISYMI) 1441 * + NVIR(ISYMA)*(I-1) 1442 * + A 1443 1444 KAIKJ = ISAIKJ(ISYAIK,ISYMJ) 1445 * + NCKI(ISYAIK)*(J-1) 1446 * + ISAIK(ISYAI,ISYMK) 1447 * + NT1AM(ISYAI)*(K-1) 1448 * + IT1AM(ISYMA,ISYMI) 1449 * + NVIR(ISYMA)*(I-1) 1450 * + A 1451 1452 WMAT(KAIKJ) = WMAT(KAIKJ) 1453 * - FOCKY(KKB)*T2TP(KAIJD) 1454 END DO 1455 END DO 1456 END DO 1457 END DO 1458 END DO 1459 END DO 1460 1461 1462C 1463C (3) wmat(aikj) = wmat(aikj) + focky(k,D)*t2(ai,Bj) 1464C 1465 ISYMK = MULD2H(ISYFKY,ISYMD) 1466 ISYAIJ = MULD2H(ISYMT2,ISYMB) 1467 DO ISYMJ = 1,NSYM 1468 ISYAI = MULD2H(ISYAIJ,ISYMJ) 1469 ISYAIK = MULD2H(ISYAI,ISYMK) 1470 DO ISYMI = 1,NSYM 1471 ISYMA = MULD2H(ISYAI,ISYMI) 1472 DO J = 1,NRHF(ISYMJ) 1473 DO K = 1,NRHF(ISYMK) 1474 KKD = IFCVIR(ISYMK,ISYMD) + NORB(ISYMK)*(D - 1) + K 1475 DO I = 1,NRHF(ISYMI) 1476 DO A = 1,NVIR(ISYMA) 1477 1478 KAIJB = IT2SP(ISYAIJ,ISYMB) 1479 * + NCKI(ISYAIJ)*(B-1) 1480 * + ISAIK(ISYAI,ISYMJ) 1481 * + NT1AM(ISYAI)*(J-1) 1482 * + IT1AM(ISYMA,ISYMI) 1483 * + NVIR(ISYMA)*(I-1) 1484 * + A 1485 1486 KAIKJ = ISAIKJ(ISYAIK,ISYMJ) 1487 * + NCKI(ISYAIK)*(J-1) 1488 * + ISAIK(ISYAI,ISYMK) 1489 * + NT1AM(ISYAI)*(K-1) 1490 * + IT1AM(ISYMA,ISYMI) 1491 * + NVIR(ISYMA)*(I-1) 1492 * + A 1493 1494 WMAT(KAIKJ) = WMAT(KAIKJ) 1495 * + FOCKY(KKD)*T2TP(KAIJB) 1496 END DO 1497 END DO 1498 END DO 1499 END DO 1500 END DO 1501 END DO 1502 1503C 1504C (4) wmat(aikj) = wmat(aikj) - focky(j,D)*t2(ai,Bk) 1505C 1506 ISYMJ = MULD2H(ISYFKY,ISYMD) 1507 ISYAIK = MULD2H(ISYMT2,ISYMB) 1508 DO ISYMK = 1,NSYM 1509 ISYAI = MULD2H(ISYAIK,ISYMK) 1510 DO ISYMI = 1,NSYM 1511 ISYMA = MULD2H(ISYAI,ISYMI) 1512 DO J = 1,NRHF(ISYMJ) 1513 KJD = IFCVIR(ISYMJ,ISYMD) + NORB(ISYMJ)*(D - 1) + J 1514 DO K = 1,NRHF(ISYMK) 1515 DO I = 1,NRHF(ISYMI) 1516 DO A = 1,NVIR(ISYMA) 1517 1518 KAIKB = IT2SP(ISYAIK,ISYMB) 1519 * + NCKI(ISYAIK)*(B-1) 1520 * + ISAIK(ISYAI,ISYMK) 1521 * + NT1AM(ISYAI)*(K-1) 1522 * + IT1AM(ISYMA,ISYMI) 1523 * + NVIR(ISYMA)*(I-1) 1524 * + A 1525 1526 KAIKJ = ISAIKJ(ISYAIK,ISYMJ) 1527 * + NCKI(ISYAIK)*(J-1) 1528 * + ISAIK(ISYAI,ISYMK) 1529 * + NT1AM(ISYAI)*(K-1) 1530 * + IT1AM(ISYMA,ISYMI) 1531 * + NVIR(ISYMA)*(I-1) 1532 * + A 1533 1534 WMAT(KAIKJ) = WMAT(KAIKJ) 1535 * - FOCKY(KJD)*T2TP(KAIKB) 1536 END DO 1537 END DO 1538 END DO 1539 END DO 1540 END DO 1541 END DO 1542 1543 CALL QEXIT('WBARBD_T2') 1544C 1545 RETURN 1546 END 1547C /* Deck cc3_w3_omega1 */ 1548 SUBROUTINE CC3_W3_OMEGA1(OMEGA1,ISYRES,WMAT,TMAT,ISYMIM, 1549 * XOOOO,XOVVO,XOOVV,XVVVV,ISYINT, 1550 * T2TP,ISYMT2,WORK,LWORK,LENSQ,INDSQ, 1551 * ISYMIB,IB,ISYMID,ID,W3X) 1552C 1553C Written by K. Hald, Spring 2002. 1554C 1555C Calculate the L3 contributions to omega1. 1556C 1557 IMPLICIT NONE 1558C 1559 LOGICAL W3X 1560C 1561 INTEGER ISYMIM, ISYINT, ISYMT2, LWORK, ISYMIB, IB, ISYMID, ID 1562 INTEGER LENSQ, INDSQ(LENSQ,6) 1563 INTEGER ISYMBD, ISCKIJ, ISYCKM, LENGTH, ISYTMP, KSCR1 1564 INTEGER KEND1, LWRK1, KEND2, LWRK2, ISYMCK, ISYMIJ, ISYMDM 1565 INTEGER ISYMM, KOFF1, KOFF2, KOFF3, NTOTIJ, NTOTCK 1566 INTEGER KT2TMP, ISYMCM, ISYMK, ISYMC 1567 INTEGER ISYRES, ISYMI, ISYOOO, NTOIJK, NTOTB, NTOTI, NBI 1568 INTEGER ISYVVV, ISYEIJ, ISYMKM, ISYME, KSCR2, ISYMEK 1569 INTEGER ISYMCE, ISYMAC, ISYMA, NTOTCE, NTOTA, ISYENI, ISYMEN 1570 INTEGER ISYDLM, ISYMN, NTODLM, NTOTE, ISYMDN, ISYDNI 1571 INTEGER NTOTEN, ISYENF, ISYELM, ISYMLM, ISYML, ISYMFN, ISYFNI 1572 INTEGER ISYMEL, ISYLMI, NTOTFN, NTOTLM, ISYMEI, KSCR3, ISYMF 1573 INTEGER ISYMFI, ISYMDL, ISYMIN, NTOTDL, NTOTIN, ISYMD, ISYTMP2 1574 INTEGER ISYMMN, ISYAMN, ISYDMN, NTOTMN, ISYBMN, ISYMBN, ISYMDI 1575 INTEGER KMIMAT, KRMAT, KSORT, ISWMAT, ISYAON, ISAONM, KOFFTM 1576 INTEGER KOFFTP, KOFFRE, NL, NAON, ISONM, KOFFI, KOFFAI, NONM, NA 1577 INTEGER ISRMAT, ISYMBA, ISYMDLM, NTOTDLM, ISYMT2W, ISYMBAE 1578 INTEGER ISYMAE, KT2W 1579 INTEGER ISYMMLE, ISYMENI, ISYMLNI, KGD, KT2B, KMLNI, ISYMDNI 1580 INTEGER ISYMNI, ISYMMLN, ISYMML, NTOTML, NTOTN, NTOTMLN 1581 INTEGER ISYMEML, ISYMEM, ISYMBIN, ISYMBI 1582 INTEGER ISYMBLM, ISYFINM, ISYMMI, ISYBMI, ISYMBM, ISYMFIN 1583 INTEGER NTOTFIN, NTOTFNK, NTOTL, IOPT, ISYMFMN, NTOTFMN, ISYMB 1584 INTEGER ISYFNIM, ISYMFNI, NTOTFNI, ISYMFNM, NTOTFNM, KSCR4 1585C 1586#if defined (SYS_CRAY) 1587 REAL ZERO, ONE 1588 REAL OMEGA1(*), WMAT(*), TMAT(*) 1589 REAL XOOOO(*), XOVVO(*), XOOVV(*), XVVVV(*) 1590 REAL T2TP(*), WORK(LWORK) 1591 REAL ZERO, ONE, DDOT, XNORM 1592#else 1593 DOUBLE PRECISION OMEGA1(*), WMAT(*), TMAT(*) 1594 DOUBLE PRECISION XOOOO(*), XOVVO(*), XOOVV(*), XVVVV(*) 1595 DOUBLE PRECISION T2TP(*), WORK(LWORK) 1596 DOUBLE PRECISION ZERO, ONE, DDOT, XNORM 1597#endif 1598C 1599#include "priunit.h" 1600#include "ccorb.h" 1601#include "ccsdinp.h" 1602#include "ccsdsym.h" 1603C 1604 PARAMETER (ZERO= 0.0D0, ONE= 1.0D0) 1605C 1606C If W3X = .TRUE., then WMAT contains wbar_X and all terms in 1607C 1-6 are included 1608C else 1609C we assume that we get tbar_0 in as WMAT and calculate only 1610C those terms in 1-6 which are marked with star (*). 1611C 1612C Note Wbar_X and Tbar_0 is stored in the same way for fixed IB and ID 1613C 1614C 1615C CALCULATE THE TERMS 1-6 USING W INTERMEDIATE 1616C 1617C 1. 1618C + L^{dae}_{lno} t^{de}_{lm} g_{inmo} 1619C 2. 1620C + L^{dfg}_{lim} t^{de}_{lm} g_{fage} 1621C 3. 1622C - L^{daf}_{lmn} t^{de}_{lm} g_{iefn} 1623C 4. 1624C - L^{daf}_{lnm} t^{de}_{lm} g_{infe} 1625C 5. 1626C - L^{def}_{lin} t^{de}_{lm} g_{mafn} 1627C 6. 1628C - L^{def}_{lni} t^{de}_{lm} g_{mnfa} 1629C 1630C 1. Calculate contribution from g_{oooo} 1631C ( Wae(dlon)(*) + Wad(eoln) + Wed(anlo) )* t(dl,em) * g(inmo) 1632C 1633C 2. Calculate contribution from g_{vvvv} 1634C ( Wdg(fiml) + Wdf(gmil) + Wfg(dlmi)(*) )* t(dl,em) * g(fage) 1635C 1636C 3. Calculate contributions from g_{ovvo} 1637C -( Waf(dlnm)(*) + Wad(fnlm) + Wdf(amnl) )* t(dl,em) * g(iefn) 1638C 1639C 4. Calculate contributions from g_{oovv} 1640C -( Waf(dlmn)(*) + Wad(fmln) + Wfd(anlm) )* t(dl,em) * g(infe) 1641C 1642C 5. Calculate contributions from g_{ovvo} 1643C -( Wfe(dlin)(*) + Wfd(eiln) + Wde(fnil) )* t(dl,em) * g(mafn) 1644C 1645C 6. Calculate contributions from g_{oovv} 1646C -( Wfe(dlni)(*) + Wfd(enli) + Wde(finl) )* t(dl,em) * g(mnfa) 1647C 1648 CALL QENTER('CC3_W3_OMEGA1') 1649C 1650 ISYTMP = MULD2H(ISYMIM,ISYMT2) 1651 ISYTMP2 = MULD2H(ISYTMP,ISYINT) 1652 IF (ISYRES .NE. ISYTMP2) THEN 1653 CALL QUIT('Symmetry mimatch in CC3_W3_OMEGA1') 1654 ENDIF 1655C 1656 ISYMBD = MULD2H(ISYMIB,ISYMID) 1657 ISCKIJ = MULD2H(ISYMIM,ISYMBD) 1658 ISWMAT = ISCKIJ 1659C 1660 LENGTH = NCKIJ(ISCKIJ) 1661 IF (LENSQ.NE.LENGTH) THEN 1662 WRITE(LUPRI,*)'CC3_W3_OMEGA1, SYMMETRY MISMATCH : LENGTH 1663 * NE LENSQ' 1664 CALL QUIT(' CC3_W3_OMEGA1, LENGTH NE LENSQ') 1665 END IF 1666C 1667C================================================ 1668C 1. Calculate contribution from g_{oooo} 1669C + L^{dae}_{lno} t^{de}_{lm} g_{inmo} = 1670C ( Wae(dlon) + Wad(eoln) + Wed(anlo) )* t(dl,em) * g(inmo) 1671C================================================ 1672C 1673 ISYCKM = MULD2H(ISYMT2,ISYMID) 1674 ISYTMP = MULD2H(ISCKIJ,ISYCKM) ! Symmetry of intermediate 1675C 1676 KSCR1 = 1 1677 KEND1 = KSCR1 + NMAIJK(ISYTMP) 1678 LWRK1 = LWORK - KEND1 1679C 1680 IF (LWRK1 .LT. 0) THEN 1681 CALL QUIT('Out of memory in CC3_W3_OMEGA1') 1682 ENDIF 1683C 1684 CALL DZERO(WORK(KSCR1),NMAIJK(ISYTMP)) 1685C 1686C-------------------------------------------- 1687C First contribution to intermediate 1688C-------------------------------------------- 1689C 1690C 1691C Wae(dlon) * t(dl,em) * g(inmo) 1692C 1693C WAD(ckij) * t2tp(ckmD) * g(Ijmi) 1694C GAD(ij,m) 1695C 1696 DO I = 1, LENGTH 1697 TMAT(I) = WMAT(I) 1698 ENDDO 1699C 1700 IF (NSYM .GT. 1) THEN 1701 IF (LWRK1 .LT. LENGTH) THEN 1702 CALL QUIT('Out of memory in CC3_W3_OMEGA1 (CC_GATHER-1)') 1703 ENDIF 1704 CALL CC_GATHER(LENGTH,WORK(KEND1),TMAT,INDSQ(1,6)) 1705 CALL DCOPY(LENGTH,WORK(KEND1),1,TMAT,1) 1706 ENDIF 1707C 1708 DO ISYMCK = 1, NSYM 1709C 1710 ISYMIJ = MULD2H(ISCKIJ,ISYMCK) 1711 ISYMM = MULD2H(ISYCKM,ISYMCK) 1712C 1713 KOFF1 = ISAIKL(ISYMCK,ISYMIJ) 1714 * + 1 1715 KOFF2 = IT2SP(ISYCKM,ISYMID) 1716 * + NCKI(ISYCKM)*(ID-1) 1717 * + ICKI(ISYMCK,ISYMM) 1718 * + 1 1719 KOFF3 = KSCR1 1720 * + IMAIJK(ISYMIJ,ISYMM) 1721C 1722 NTOTIJ = MAX(1,NMATIJ(ISYMIJ)) 1723 NTOTCK = MAX(1,NT1AM(ISYMCK)) 1724C 1725 CALL DGEMM('T','N',NMATIJ(ISYMIJ),NRHF(ISYMM),NT1AM(ISYMCK), 1726 * ONE,TMAT(KOFF1),NTOTCK,T2TP(KOFF2),NTOTCK, 1727 * ONE,WORK(KOFF3),NTOTIJ) 1728C 1729 ENDDO 1730c 1731* write(lupri,*)'kscr1(on,m) in cc3_w3_lhtr ib,id ',ib,id 1732* call output(work(kscr1),1,nmatij(1),1,nrhf(1), 1733* * nmatij(1),nrhf(1), 1734* * 1,lupri) 1735C 1736C-------------------------------------------- 1737C Second contribution to intermediate 1738C-------------------------------------------- 1739C 1740C Wad(eoln) t(dl,em) * g(inmo) 1741C 1742C WAD(ckij * tD(ckm) * g(Ijmi) 1743C GAD(ij,m) 1744C 1745 ! Calculate only when contracting with first-order multipliers 1746 IF (W3X) THEN 1747C 1748 DO I = 1, LENGTH 1749 TMAT(I) = WMAT(INDSQ(I,1)) 1750 ENDDO 1751C 1752 IF (NSYM .GT. 1) THEN 1753 IF (LWRK1 .LT. LENGTH) THEN 1754 CALL QUIT('Out of memory in CC3_W3_OMEGA1 (CC_GATHER-2)') 1755 ENDIF 1756 CALL CC_GATHER(LENGTH,WORK(KEND1),TMAT,INDSQ(1,6)) 1757 CALL DCOPY(LENGTH,WORK(KEND1),1,TMAT,1) 1758 ENDIF 1759C 1760 KT2TMP = KEND1 1761 KEND2 = KT2TMP + NCKI(ISYCKM) 1762 LWRK2 = LWORK - KEND2 1763C 1764 IF (LWRK2 .LT. 0) THEN 1765 CALL QUIT('Memory exceeded in CC3_W3_OMEGA1') 1766 ENDIF 1767C 1768C sort t2tp(ckmD) as TD(cmk) 1769C 1770 DO ISYMK = 1, NSYM 1771 ISYMCM = MULD2H(ISYCKM,ISYMK) 1772 DO ISYMC = 1, NSYM 1773 ISYMM = MULD2H(ISYMCM,ISYMC) 1774 ISYMCK = MULD2H(ISYMC,ISYMK) 1775C 1776 DO K = 1, NRHF(ISYMK) 1777 DO M = 1, NRHF(ISYMM) 1778C 1779 KOFF1 = IT2SP(ISYCKM,ISYMID) 1780 * + NCKI(ISYCKM)*(ID-1) 1781 * + ICKI(ISYMCK,ISYMM) 1782 * + NT1AM(ISYMCK)*(M-1) 1783 * + IT1AM(ISYMC,ISYMK) 1784 * + NVIR(ISYMC)*(K-1) 1785 * + 1 1786 KOFF2 = KT2TMP 1787 * + ICKI(ISYMCM,ISYMK) 1788 * + NT1AM(ISYMCM)*(K-1) 1789 * + IT1AM(ISYMC,ISYMM) 1790 * + NVIR(ISYMC)*(M-1) 1791C 1792 CALL DCOPY(NVIR(ISYMC),T2TP(KOFF1),1,WORK(KOFF2),1) 1793C 1794 ENDDO 1795 ENDDO 1796 ENDDO 1797 ENDDO 1798C 1799 DO ISYMCK = 1, NSYM 1800C 1801 ISYMIJ = MULD2H(ISCKIJ,ISYMCK) 1802 ISYMM = MULD2H(ISYCKM,ISYMCK) 1803C 1804 KOFF1 = ISAIKL(ISYMCK,ISYMIJ) 1805 * + 1 1806 KOFF2 = KT2TMP 1807 * + ICKI(ISYMCK,ISYMM) 1808 KOFF3 = KSCR1 1809 * + IMAIJK(ISYMIJ,ISYMM) 1810C 1811 NTOTIJ = MAX(1,NMATIJ(ISYMIJ)) 1812 NTOTCK = MAX(1,NT1AM(ISYMCK)) 1813C 1814 CALL DGEMM('T','N',NMATIJ(ISYMIJ),NRHF(ISYMM),NT1AM(ISYMCK), 1815 * ONE,TMAT(KOFF1),NTOTCK,WORK(KOFF2),NTOTCK, 1816 * ONE,WORK(KOFF3),NTOTIJ) 1817C 1818 ENDDO 1819C 1820 END IF 1821c 1822* write(lupri,*)'kscr1(on,m) in cc3_w3_lhtr ib,id ',ib,id 1823* call output(work(kscr1),1,nmatij(1),1,nrhf(1), 1824* * nmatij(1),nrhf(1), 1825* * 1,lupri) 1826C 1827C------------------------------------------------ 1828C Contract the intermediate with g_{oooo} 1829C------------------------------------------------ 1830C 1831 ISYMI = MULD2H(ISYRES,ISYMIB) 1832 ISYOOO = MULD2H(ISYINT,ISYMI) 1833 DO I = 1, NRHF(ISYMI) 1834 NBI = IT1AM(ISYMIB,ISYMI) + NVIR(ISYMIB)*(I-1) + IB 1835 KOFF1 = I3ORHF(ISYOOO,ISYMI) 1836 * + NMAIJK(ISYOOO)*(I-1) 1837 * + 1 1838C 1839C 1840C 1.1 (+ 1.2) omega contribution addomega 1841C 1842 OMEGA1(NBI) = OMEGA1(NBI) 1843 * - DDOT(NMAIJK(ISYOOO),XOOOO(KOFF1),1,WORK(KSCR1),1) 1844 ENDDO 1845C 1846C----------------------------------------- 1847C Wed(anlo) * t(dl,em) * g(inmo) 1848C 1849C WBD(anlo) * t2t2(DlmB) * g(inmo) 1850C 1851C TMAT(aonl) * T(lm) * G(onmi) 1852C 1853C R(aonm) 1854C----------------------------------------- 1855C 1856C TMAT(aonl) = WBD(anlo) 1857C 1858 ! Calculate only when contracting with first-order multipliers 1859 IF (W3X) THEN 1860C 1861 DO I = 1, LENGTH 1862 TMAT(I) = WMAT(INDSQ(I,2)) 1863 END DO 1864C 1865C ------------------------------- 1866C Sort T^{DB}_{lm} as T_{lm} 1867C ------------------------------- 1868C 1869 ISYMBD = MULD2H(ISYMIB,ISYMID) 1870 ISYMLM = MULD2H(ISYMBD,ISYMT2) 1871 ISYDLM = MULD2H(ISYMLM,ISYMID) 1872 ISRMAT = MULD2H(ISWMAT,ISYMLM) 1873C 1874 KMIMAT = 1 1875 KRMAT = KMIMAT + NMATIJ(ISYMLM) 1876 KSORT = KRMAT + N3VOOO(ISRMAT) 1877 KEND1 = KSORT + N3VOOO(ISRMAT) 1878 LWRK1 = LWORK - KEND1 1879C 1880 IF (LWRK1 .LT.0 ) THEN 1881 CALL QUIT('Out of memory in CC3_W3_OMEGA1 (sort)') 1882 ENDIF 1883C 1884 DO ISYMM = 1,NSYM 1885 ISYML = MULD2H(ISYMLM,ISYMM) 1886 ISYMDL = MULD2H(ISYDLM,ISYMM) 1887 DO M = 1,NRHF(ISYMM) 1888 DO L = 1,NRHF(ISYML) 1889 KOFF1 = IT2SP(ISYDLM,ISYMIB) 1890 * + NCKI(ISYDLM)*(IB-1) 1891 * + ICKI(ISYMDL,ISYMM) 1892 * + NT1AM(ISYMDL)*(M-1) 1893 * + IT1AM(ISYMID,ISYML) 1894 * + NVIR(ISYMID)*(L-1) 1895 * + ID 1896 KOFF2 = IMATIJ(ISYML,ISYMM) 1897 * + NRHF(ISYML)*(M-1) 1898 * + L 1899 WORK(KMIMAT-1+KOFF2) = T2TP(KOFF1) 1900 ENDDO 1901 ENDDO 1902 ENDDO 1903C 1904C R(aonm) = sum_l ( TMAT(aonl) * T(lm) ) 1905C 1906 DO ISYMM = 1,NSYM 1907 ISYML = MULD2H(ISYMLM,ISYMM) 1908 ISYAON = MULD2H(ISWMAT,ISYML) 1909 ISAONM = MULD2H(ISYAON,ISYMM) 1910 KOFFTM = ISAIKJ(ISYAON,ISYML) + 1 1911 KOFFTP = IMATIJ(ISYML,ISYMM) + KMIMAT 1912 KOFFRE = ISAIKJ(ISYAON,ISYMM) + KRMAT 1913 NL = MAX(1,NRHF(ISYML)) 1914 NAON = MAX(1,NCKI(ISYAON)) 1915 1916C 1917 CALL DGEMM('N','N',NCKI(ISYAON),NRHF(ISYMM), 1918 * NRHF(ISYML),ONE,TMAT(KOFFTM),NAON, 1919 * WORK(KOFFTP),NL,ZERO,WORK(KOFFRE),NAON) 1920C 1921 END DO 1922* write(lupri,*)'krmat(aon,m) in cc3_w3_lhtr ib,id ',ib,id 1923* call output(work(krmat),1,ncki(1),1,nrhf(1), 1924* * ncki(1),nrhf(1), 1925* * 1,lupri) 1926 1927C 1928C omega(ai) = sum_onm ( R(aonm) * G(onmi) ) 1929C 1930C 1931 IF (NSYM .GT. 1) THEN 1932 CALL CC3_SRTVOOO(WORK(KSORT),WORK(KRMAT),ISRMAT) 1933 CALL DCOPY(N3VOOO(ISRMAT),WORK(KSORT),1,WORK(KRMAT),1) 1934 ENDIF 1935C 1936 1937 DO ISYMI = 1,NSYM 1938 ISYMA = MULD2H(ISYRES,ISYMI) 1939 ISONM = MULD2H(ISYINT,ISYMI) 1940 KOFFRE = I3VOOO(ISYMA,ISONM) + KRMAT 1941 KOFFI = I3ORHF(ISONM,ISYMI) + 1 1942 KOFFAI = IT1AM(ISYMA,ISYMI) + 1 1943 NONM = MAX(1,NMAIJK(ISONM)) 1944 NA = MAX(1,NVIR(ISYMA)) 1945 1946C 1947C 1.3 omega contribution addomega 1948C 1949 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI), 1950 * NMAIJK(ISONM),-ONE,WORK(KOFFRE),NA, 1951 * XOOOO(KOFFI),NONM,ONE,OMEGA1(KOFFAI),NA) 1952 END DO 1953C 1954 END IF 1955C 1956C============================================= 1957C 2. Calculate contribution from g_{vvvv} 1958C + L^{dfg}_{lim} t^{de}_{lm} g_{fage} = 1959C ( Wdg(fiml) + Wdf(gmil) + Wfg(dlmi) )* t(dl,em) * g(fage) 1960C============================================= 1961C 1962 ! Calculate only when contracting with first-order multipliers 1963 IF (W3X) THEN 1964C 1965 ISYCKM = MULD2H(ISYMT2,ISYMIB) 1966 ISYEIJ = ISYCKM 1967 ISYTMP = MULD2H(ISCKIJ,ISYEIJ) 1968C 1969 KSCR1 = 1 1970 KT2TMP = KSCR1 + NCKATR(ISYTMP) 1971 KEND1 = KT2TMP + NCKI(ISYEIJ) 1972 LWRK1 = LWORK - KEND1 1973C 1974 IF (LWRK1 .LT. 0) THEN 1975 CALL QUIT('Out of memory in CC3_W3_OMEGA1 (VVVV T2-sort)') 1976 ENDIF 1977C 1978 CALL DZERO(WORK(KSCR1),NCKATR(ISYTMP)) 1979C 1980C---------------- 1981C Sort T2 tB(km,c) = t2tp(ckmB) 1982C---------------- 1983C 1984 DO ISYMK = 1, NSYM 1985 ISYMCM = MULD2H(ISYCKM,ISYMK) 1986 DO ISYMC = 1, NSYM 1987 ISYMM = MULD2H(ISYMCM,ISYMC) 1988 ISYMKM = MULD2H(ISYMK,ISYMM) 1989 ISYMCK = MULD2H(ISYMK,ISYMC) 1990C 1991 DO K = 1, NRHF(ISYMK) 1992 DO M = 1, NRHF(ISYMM) 1993 KOFF1 = IT2SP(ISYCKM,ISYMIB) 1994 * + NCKI(ISYCKM)*(IB-1) 1995 * + ICKI(ISYMCK,ISYMM) 1996 * + NT1AM(ISYMCK)*(M-1) 1997 * + IT1AM(ISYMC,ISYMK) 1998 * + NVIR(ISYMC)*(K-1) 1999 * + 1 2000 KOFF2 = KT2TMP - 1 2001 * + IMAIJA(ISYMKM,ISYMC) 2002 * + IMATIJ(ISYMK,ISYMM) 2003 * + NRHF(ISYMK)*(M-1) 2004 * + K 2005C 2006 CALL DCOPY(NVIR(ISYMC),T2TP(KOFF1),1, 2007 * WORK(KOFF2),NMATIJ(ISYMKM)) 2008 ENDDO 2009 ENDDO 2010 ENDDO 2011 ENDDO 2012C 2013C-------------------------------- 2014C First intermediate 2015C-------------------------------- 2016C 2017C Wdg(fiml) * t(dl,em) * g(fage) 2018C 2019C WBD(ck,ij) * tB(ij,e) * g(caDe) 2020C 2021C TMAT(ck,ij) * tB(ij,e) * gD(cae) 2022C G(ck,e) * GD(ce,a) 2023C 2024 DO I = 1, LENGTH 2025 TMAT(I) = WMAT(I) 2026 ENDDO 2027C 2028 IF (NSYM .GT. 1) THEN 2029 IF (LWRK1 .LT. LENGTH) THEN 2030 CALL QUIT('Out of memory in CC3_W3_OMEGA1 (CC_GATHER-3)') 2031 ENDIF 2032 CALL CC_GATHER(LENGTH,WORK(KEND1),TMAT,INDSQ(1,6)) 2033 CALL DCOPY(LENGTH,WORK(KEND1),1,TMAT,1) 2034 ENDIF 2035C 2036 DO ISYMCK = 1, NSYM 2037 ISYMIJ = MULD2H(ISCKIJ,ISYMCK) 2038 ISYME = MULD2H(ISYEIJ,ISYMIJ) 2039C 2040 KOFF1 = ISAIKL(ISYMCK,ISYMIJ) 2041 * + 1 2042 KOFF2 = KT2TMP 2043 * + IMAIJA(ISYMIJ,ISYME) 2044 KOFF3 = KSCR1 2045 * + ICKATR(ISYMCK,ISYME) 2046C 2047 NTOTCK = MAX(1,NT1AM(ISYMCK)) 2048 NTOTIJ = MAX(1,NMATIJ(ISYMIJ)) 2049C 2050 CALL DGEMM('N','N',NT1AM(ISYMCK),NVIR(ISYME), 2051 * NMATIJ(ISYMIJ),ONE,TMAT(KOFF1),NTOTCK, 2052 * WORK(KOFF2),NTOTIJ,ONE,WORK(KOFF3), 2053 * NTOTCK) 2054 ENDDO 2055C 2056C--------------------- 2057C Sort result. 2058C--------------------- 2059C 2060 KSCR2 = KEND1 2061 KEND2 = KSCR2 + NCKATR(ISYTMP) 2062 LWRK2 = LWORK - KEND2 2063C 2064 IF (LWRK2 .LT. 0) THEN 2065 CALL QUIT('Out of memory in CC3_W3_OMEGA1 (Sorting-1)') 2066 ENDIF 2067C 2068 DO ISYMC = 1, NSYM 2069 ISYMEK = MULD2H(ISYMC,ISYTMP) 2070 DO ISYMK = 1, NSYM 2071 ISYME = MULD2H(ISYMK,ISYMEK) 2072 ISYMCE = MULD2H(ISYMC,ISYME) 2073 ISYMCK = MULD2H(ISYMC,ISYMK) 2074C 2075 DO K = 1, NRHF(ISYMK) 2076 DO E = 1, NVIR(ISYME) 2077C 2078 KOFF1 = KSCR1 2079 * + ICKATR(ISYMCK,ISYME) 2080 * + NT1AM(ISYMCK)*(E-1) 2081 * + IT1AM(ISYMC,ISYMK) 2082 * + NVIR(ISYMC)*(K-1) 2083 KOFF2 = KSCR2 2084 * + ICKASR(ISYMCE,ISYMK) 2085 * + NMATAB(ISYMCE)*(K-1) 2086 * + IMATAB(ISYMC,ISYME) 2087 * + NVIR(ISYMC)*(E-1) 2088C 2089 CALL DCOPY(NVIR(ISYMC),WORK(KOFF1),1, 2090 * WORK(KOFF2),1) 2091C 2092 ENDDO 2093 ENDDO 2094 ENDDO 2095 ENDDO 2096C 2097 CALL DCOPY(NCKATR(ISYTMP),WORK(KSCR2),1,WORK(KSCR1),1) 2098C 2099C---------------------------------------- 2100C Sort and contract with integral. 2101C---------------------------------------- 2102C 2103 ISYVVV = MULD2H(ISYINT,ISYMID) 2104C 2105 KSCR2 = KEND1 2106 KEND2 = KSCR2 + NMAABC(ISYVVV) 2107 LWRK2 = LWORK - KEND2 2108C 2109 IF (LWRK2 .LT. 0) THEN 2110 CALL QUIT('Out of memory in CC3_W3_OMEGA1 (Sort/Contract)') 2111 ENDIF 2112C 2113 DO ISYME = 1, NSYM 2114 ISYMAC = MULD2H(ISYVVV,ISYME) 2115 DO ISYMC = 1, NSYM 2116 ISYMA = MULD2H(ISYMAC,ISYMC) 2117 ISYMCE = MULD2H(ISYMC,ISYME) 2118C 2119 DO A = 1, NVIR(ISYMA) 2120 DO E = 1, NVIR(ISYME) 2121C 2122 KOFF1 = IMAABC(ISYMAC,ISYME) 2123 * + NMATAB(ISYMAC)*(E-1) 2124 * + IMATAB(ISYMC,ISYMA) 2125 * + NVIR(ISYMC)*(A-1) 2126 * + 1 2127 KOFF2 = KSCR2 2128 * + IMAABC(ISYMCE,ISYMA) 2129 * + NMATAB(ISYMCE)*(A-1) 2130 * + IMATAB(ISYMC,ISYME) 2131 * + NVIR(ISYMC)*(E-1) 2132C 2133 CALL DCOPY(NVIR(ISYMC),XVVVV(KOFF1),1, 2134 * WORK(KOFF2),1) 2135C 2136 ENDDO 2137 ENDDO 2138 ENDDO 2139 ENDDO 2140C 2141 DO ISYMA = 1, NSYM 2142 ISYMI = MULD2H(ISYMA,ISYRES) 2143 ISYMCE = MULD2H(ISYMA,ISYVVV) 2144C 2145 KOFF1 = KSCR2 2146 * + IMAABC(ISYMCE,ISYMA) 2147 KOFF2 = KSCR1 2148 * + ICKASR(ISYMCE,ISYMI) 2149 KOFF3 = IT1AM(ISYMA,ISYMI) 2150 * + 1 2151C 2152 NTOTCE = MAX(1,NMATAB(ISYMCE)) 2153 NTOTA = MAX(1,NVIR(ISYMA)) 2154C 2155C 2156C 2157C 2.1 omega contribution addomega 2158C 2159 CALL DGEMM('T','N',NVIR(ISYMA),NRHF(ISYMI),NMATAB(ISYMCE), 2160 * -ONE,WORK(KOFF1),NTOTCE,WORK(KOFF2),NTOTCE, 2161 * ONE,OMEGA1(KOFF3),NTOTA) 2162 ENDDO 2163C 2164C--------------------------------------- 2165C Second contribution. 2166C--------------------------------------- 2167C 2168C Wdf(gmil) * t(dl,em) * g(fage) 2169C 2170C WBD(gmil) * t(emlB) * g(geDa) 2171C 2172C TMAT(giml) * tB(mle) * gD(gea) 2173C 2174C TMAT(ckij) * tB(ij,e) * gD(ce,a) 2175C G(ck,e) 2176C sort G(ci,e) as M(ce,i) 2177C 2178C eta(ai) = eta(ai) + gD(ce,a) * M(ce,i) 2179C 2180 ISYCKM = MULD2H(ISYMT2,ISYMIB) 2181 ISYEIJ = ISYCKM 2182 ISYTMP = MULD2H(ISCKIJ,ISYEIJ) 2183C 2184 KSCR1 = 1 2185 KT2TMP = KSCR1 + NCKATR(ISYTMP) 2186 KEND1 = KT2TMP + NCKI(ISYCKM) 2187 LWRK1 = LWORK - KEND1 2188C 2189 IF (LWRK1 .LT. 0) THEN 2190 CALL QUIT('Out of memory in CC3_W3_OMEGA1 (T2-sort-2)') 2191 ENDIF 2192C 2193 CALL DZERO(WORK(KSCR1),NCKATR(ISYTMP)) 2194C 2195C---------------- 2196C Sort T2 tB(km,c) = t2tp(ckmB) 2197C---------------- 2198C 2199 DO ISYMK = 1, NSYM 2200 ISYMCM = MULD2H(ISYCKM,ISYMK) 2201 DO ISYMC = 1, NSYM 2202 ISYMM = MULD2H(ISYMCM,ISYMC) 2203 ISYMKM = MULD2H(ISYMK,ISYMM) 2204 ISYMCK = MULD2H(ISYMK,ISYMC) 2205C 2206 DO K = 1, NRHF(ISYMK) 2207 DO M = 1, NRHF(ISYMM) 2208 KOFF1 = IT2SP(ISYCKM,ISYMIB) 2209 * + NCKI(ISYCKM)*(IB-1) 2210 * + ICKI(ISYMCK,ISYMM) 2211 * + NT1AM(ISYMCK)*(M-1) 2212 * + IT1AM(ISYMC,ISYMK) 2213 * + NVIR(ISYMC)*(K-1) 2214 * + 1 2215 KOFF2 = KT2TMP - 1 2216 * + IMAIJA(ISYMKM,ISYMC) 2217 * + IMATIJ(ISYMK,ISYMM) 2218 * + NRHF(ISYMK)*(M-1) 2219 * + K 2220C 2221 CALL DCOPY(NVIR(ISYMC),T2TP(KOFF1),1, 2222 * WORK(KOFF2),NMATIJ(ISYMKM)) 2223 ENDDO 2224 ENDDO 2225 ENDDO 2226 ENDDO 2227C 2228C-------------------------------- 2229C Second intermediate 2230C-------------------------------- 2231C 2232 DO I = 1, LENGTH 2233 TMAT(I) = WMAT(INDSQ(I,1)) 2234 ENDDO 2235C 2236 IF (NSYM .GT. 1) THEN 2237 IF (LWRK1 .LT. LENGTH) THEN 2238 CALL QUIT('Out of memory in CC3_W3_OMEGA1 (CC_GATHER-3)') 2239 ENDIF 2240 CALL CC_GATHER(LENGTH,WORK(KEND1),TMAT,INDSQ(1,6)) 2241 CALL DCOPY(LENGTH,WORK(KEND1),1,TMAT,1) 2242 ENDIF 2243C 2244 DO ISYMCK = 1, NSYM 2245 ISYMIJ = MULD2H(ISCKIJ,ISYMCK) 2246 ISYME = MULD2H(ISYEIJ,ISYMIJ) 2247C 2248 KOFF1 = ISAIKL(ISYMCK,ISYMIJ) 2249 * + 1 2250 KOFF2 = KT2TMP 2251 * + IMAIJA(ISYMIJ,ISYME) 2252 KOFF3 = KSCR1 2253 * + ICKATR(ISYMCK,ISYME) 2254C 2255 NTOTCK = MAX(1,NT1AM(ISYMCK)) 2256 NTOTIJ = MAX(1,NMATIJ(ISYMIJ)) 2257C 2258 CALL DGEMM('N','N',NT1AM(ISYMCK),NVIR(ISYME), 2259 * NMATIJ(ISYMIJ),ONE,TMAT(KOFF1),NTOTCK, 2260 * WORK(KOFF2),NTOTIJ,ONE,WORK(KOFF3), 2261 * NTOTCK) 2262 ENDDO 2263C 2264C--------------------- 2265C Sort result. 2266C--------------------- 2267C 2268 KSCR2 = KEND1 2269 KEND2 = KSCR2 + NCKATR(ISYTMP) 2270 LWRK2 = LWORK - KEND2 2271C 2272 IF (LWRK2 .LT. 0) THEN 2273 CALL QUIT('Out of memory in CC3_W3_OMEGA1 (Sorting-1)') 2274 ENDIF 2275C 2276 DO ISYMC = 1, NSYM 2277 ISYMEK = MULD2H(ISYMC,ISYTMP) 2278 DO ISYMK = 1, NSYM 2279 ISYME = MULD2H(ISYMK,ISYMEK) 2280 ISYMCE = MULD2H(ISYMC,ISYME) 2281 ISYMCK = MULD2H(ISYMC,ISYMK) 2282C 2283 DO K = 1, NRHF(ISYMK) 2284 DO E = 1, NVIR(ISYME) 2285C 2286 KOFF1 = KSCR1 - 1 2287 * + ICKATR(ISYMCK,ISYME) 2288 * + NT1AM(ISYMCK)*(E-1) 2289 * + IT1AM(ISYMC,ISYMK) 2290 * + NVIR(ISYMC)*(K-1) 2291 * + 1 2292 KOFF2 = KSCR2 - 1 2293 * + ICKASR(ISYMCE,ISYMK) 2294 * + NMATAB(ISYMCE)*(K-1) 2295 * + IMATAB(ISYMC,ISYME) 2296 * + NVIR(ISYMC)*(E-1) 2297 * + 1 2298C 2299 CALL DCOPY(NVIR(ISYMC),WORK(KOFF1),1, 2300 * WORK(KOFF2),1) 2301C 2302 ENDDO 2303 ENDDO 2304 ENDDO 2305 ENDDO 2306C 2307 CALL DCOPY(NCKATR(ISYTMP),WORK(KSCR2),1,WORK(KSCR1),1) 2308C 2309C---------------------------------------- 2310C Contract with integral. 2311C---------------------------------------- 2312C 2313 ISYVVV = MULD2H(ISYINT,ISYMID) 2314C 2315 DO ISYMA = 1, NSYM 2316 ISYMI = MULD2H(ISYMA,ISYRES) 2317 ISYMCE = MULD2H(ISYMA,ISYVVV) 2318C 2319 KOFF1 = IMAABC(ISYMCE,ISYMA) 2320 * + 1 2321 KOFF2 = KSCR1 2322 * + ICKASR(ISYMCE,ISYMI) 2323 KOFF3 = IT1AM(ISYMA,ISYMI) 2324 * + 1 2325C 2326 NTOTCE = MAX(1,NMATAB(ISYMCE)) 2327 NTOTA = MAX(1,NVIR(ISYMA)) 2328C 2329C 2.2 omega contribution addomega 2330C 2331 CALL DGEMM('T','N',NVIR(ISYMA),NRHF(ISYMI),NMATAB(ISYMCE), 2332 * -ONE,XVVVV(KOFF1),NTOTCE,WORK(KOFF2),NTOTCE, 2333 * ONE,OMEGA1(KOFF3),NTOTA) 2334 ENDDO 2335C 2336 END IF 2337C 2338C third contribution 2339C 2340C----------------------------------------- 2341C Wfg(dlmi) * t(dl,em) * g(fage) 2342C 2343C WBD(dlmi) * t2tp(dlme) * g(BaDe) 2344C 2345C t2tp(dlme) * WBD(dlmi) * gD(Bae) 2346C 2347C eta (ai) = eta(ai) + gDB(ae) * t2w(ei) 2348C---------------------------------------- 2349C 2350C sort gD(Bae) as gDB(ae) 2351C 2352 ISWMAT = ISCKIJ 2353 ISYMT2W = MULD2H(ISWMAT,ISYMT2) 2354 ISYMBAE = MULD2H(ISYINT,ISYMID) 2355 ISYMAE = MULD2H(ISYMBAE,ISYMIB) 2356C 2357 KSCR2 = 1 2358 KT2W = KSCR2 + NMATAB(ISYMAE) 2359 KEND2 = KT2W + NT1AM(ISYMT2W) 2360 LWRK2 = LWORK - KEND2 2361C 2362 CALL DZERO(WORK(KT2W),NT1AM(ISYMT2W)) 2363C 2364 IF (LWRK2 .LT. 0) THEN 2365 CALL QUIT('Out of memory in CC3_W3_OMEGA1 (Sort/Contract)') 2366 ENDIF 2367C 2368 DO ISYME = 1, NSYM 2369 ISYMBA = MULD2H(ISYMBAE,ISYME) 2370 ISYMA = MULD2H(ISYMBA,ISYMIB) 2371C 2372 DO E = 1, NVIR(ISYME) 2373 DO A = 1, NVIR(ISYMA) 2374C 2375 KOFF1 = IMAABC(ISYMBA,ISYME) 2376 * + NMATAB(ISYMBA)*(E-1) 2377 * + IMATAB(ISYMIB,ISYMA) 2378 * + NVIR(ISYMIB)*(A-1) 2379 * + IB 2380 KOFF2 = KSCR2 2381 * + IMATAB(ISYMA,ISYME) 2382 * + NVIR(ISYMA)*(E-1) 2383 * + A -1 2384C 2385 WORK(KOFF2) = XVVVV(KOFF1) 2386C 2387 ENDDO 2388 ENDDO 2389 ENDDO 2390C 2391C t2w(ei) = t2tp(dlme) * WBD(dlmi) 2392C 2393 DO I = 1, LENGTH 2394 TMAT(I) = WMAT(I) 2395 ENDDO 2396C 2397 DO ISYMI = 1,NSYM 2398 ISYMDLM = MULD2H(ISWMAT,ISYMI) 2399 ISYME = MULD2H(ISYMDLM,ISYMT2) 2400C 2401 KOFF1 = IT2SP(ISYMDLM,ISYME) + 1 2402 KOFF2 = ISAIKJ(ISYMDLM,ISYMI) + 1 2403 KOFF3 = KT2W + IT1AM(ISYME,ISYMI) 2404C 2405 NTOTDLM = MAX(1,NCKI(ISYMDLM)) 2406 NTOTE = MAX(1,NVIR(ISYME)) 2407C 2408 CALL DGEMM('T','N',NVIR(ISYME),NRHF(ISYMI), 2409 * NCKI(ISYMDLM),ONE,T2TP(KOFF1),NTOTDLM, 2410 * TMAT(KOFF2),NTOTDLM,ONE,WORK(KOFF3), 2411 * NTOTE) 2412 2413C 2414 END DO 2415C 2416C eta (ai) = eta(ai) + gDB(ae) * t2w(ei) 2417C 2418 DO ISYMI = 1,NSYM 2419 ISYME = MULD2H(ISYMT2W,ISYMI) 2420 ISYMA = MULD2H(ISYMAE,ISYME) 2421C 2422 KOFF1 = KSCR2 + IMATAB(ISYMA,ISYME) 2423 KOFF2 = KT2W + IT1AM(ISYME,ISYMI) 2424 KOFF3 = IT1AM(ISYMA,ISYMI) + 1 2425C 2426 NTOTA = MAX(1,NVIR(ISYMA)) 2427 NTOTE = MAX(1,NVIR(ISYME)) 2428C 2429C 2.3 omega contribution addomega 2430C 2431 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI), 2432 * NVIR(ISYME),-ONE,WORK(KOFF1),NTOTA, 2433 * WORK(KOFF2),NTOTE,ONE,OMEGA1(KOFF3), 2434 * NTOTA) 2435 END DO 2436C 2437C================================================ 2438C 3. Calculate contributions from g_{ovvo} 2439C - L^{daf}_{lmn} t^{de}_{lm} g_{iefn} = 2440C -( Waf(dlnm) + Wad(fnlm) + Wdf(amnl) )* t(dl,em) * g(iefn) 2441C 2442C================================================ 2443C 2444C-------------------------------- 2445C First contribution 2446C-------------------------------- 2447C - Waf(dlnm) * t(dl,em) * g(iefn) 2448C 2449C t2tp(dlme) * Tmat(dlmn) xovvo(Dnie) 2450C 2451C omega(ib,i) = omega(ib,i) + t2T(en) * gD(eni) 2452C 2453 ISYMEN = MULD2H(ISCKIJ,ISYMT2) 2454 ISYMI = MULD2H(ISYRES,ISYMIB) 2455 ISYENI = MULD2H(ISYMEN,ISYMI) 2456C 2457 KSCR1 = 1 2458 KSCR2 = KSCR1 + NT1AM(ISYMEN) 2459 KEND1 = KSCR2 + NCKI(ISYENI) 2460 LWRK1 = LWORK - KEND1 2461C 2462 IF (LWRK1 .LT. 0) THEN 2463 CALL QUIT('Out of memory in CC3_W3_OMEGA1 (OVVO-1)') 2464 ENDIF 2465C 2466C t2T(en) = t2tp(dlme) * Tmat(dlmn) 2467C 2468 CALL DZERO(WORK(KSCR1),NT1AM(ISYMEN)) 2469 CALL DZERO(WORK(KSCR2),NCKI(ISYENI)) 2470C 2471 DO I = 1, LENGTH 2472 TMAT(I) = WMAT(INDSQ(I,3)) 2473 ENDDO 2474C 2475 DO ISYME = 1, NSYM 2476 ISYDLM = MULD2H(ISYME,ISYMT2) 2477 ISYMN = MULD2H(ISYMEN,ISYME) 2478C 2479 KOFF1 = IT2SP(ISYDLM,ISYME) 2480 * + 1 2481 KOFF2 = ISAIKJ(ISYDLM,ISYMN) 2482 * + 1 2483 KOFF3 = KSCR1 2484 * + IT1AM(ISYME,ISYMN) 2485C 2486 NTODLM = MAX(1,NCKI(ISYDLM)) 2487 NTOTE = MAX(1,NVIR(ISYME)) 2488C 2489 CALL DGEMM('T','N',NVIR(ISYME),NRHF(ISYMN),NCKI(ISYDLM), 2490 * -ONE,T2TP(KOFF1),NTODLM,TMAT(KOFF2),NTODLM, 2491 * ONE,WORK(KOFF3),NTOTE) 2492C 2493 ENDDO 2494C 2495C sort integrals xovvo(Dnie) as gD(eni) 2496C 2497 DO ISYME = 1, NSYM 2498 ISYMN = MULD2H(ISYMEN,ISYME) 2499 ISYMDN = MULD2H(ISYMN,ISYMID) 2500 ISYDNI = MULD2H(ISYMDN,ISYMI) 2501C 2502 DO E = 1, NVIR(ISYME) 2503 DO N = 1, NRHF(ISYMN) 2504C 2505 KOFF1 = IT2SP(ISYDNI,ISYME) 2506 * + NCKI(ISYDNI)*(E-1) 2507 * + ICKI(ISYMDN,ISYMI) 2508 * + IT1AM(ISYMID,ISYMN) 2509 * + NVIR(ISYMID)*(N-1) 2510 * + ID 2511C 2512 KOFF2 = KSCR2 - 1 2513 * + ICKI(ISYMEN,ISYMI) 2514 * + IT1AM(ISYME,ISYMN) 2515 * + NVIR(ISYME)*(N-1) 2516 * + E 2517C 2518 CALL DCOPY(NRHF(ISYMI),XOVVO(KOFF1),NT1AM(ISYMDN), 2519 * WORK(KOFF2),NT1AM(ISYMEN)) 2520C 2521 ENDDO 2522 ENDDO 2523 ENDDO 2524C 2525C omega(ib,i) = omega(ib,i) + t2T(en) * gD(eni) 2526C 2527 KOFF1 = KSCR2 2528 * + ICKI(ISYMEN,ISYMI) 2529 KOFF3 = IT1AM(ISYMIB,ISYMI) 2530 * + IB 2531C 2532 NTOTEN = MAX(1,NT1AM(ISYMEN)) 2533 NTOTB = MAX(1,NVIR(ISYMIB)) 2534C 2535C 3.1 omega contribution addomega 2536C 2537 CALL DGEMV('T',NT1AM(ISYMEN),NRHF(ISYMI),-ONE,WORK(KOFF1), 2538 * NTOTEN,WORK(KSCR1),1,ONE,OMEGA1(KOFF3),NTOTB) 2539C 2540C 2541C-------------------------------- 2542C Second contribution 2543C 2544C - Wad(fnlm) * t(dl,em) * g(iefn) 2545C 2546C TMAT(fnlm) * xovvo(fnie) t2tp(emlD) 2547C 2548C omega1(Bi) = omega(Bi) - Tx(lmie) * t2D(lme) 2549C-------------------------------- 2550C 2551 ! Calculate only when contracting with first-order multipliers 2552 IF (W3X) THEN 2553C 2554 2555 ISYTMP = MULD2H(ISCKIJ,ISYINT) 2556 ISYMI = MULD2H(ISYRES,ISYMIB) 2557 ISYENF = MULD2H(ISYINT,ISYMI) 2558 ISYELM = MULD2H(ISYMT2,ISYMID) 2559C 2560 KSCR1 = 1 2561 KEND1 = KSCR1 + NMAIJA(ISYELM) 2562 LWRK1 = LWORK - KEND1 2563C 2564 IF (LWRK1 .LT. 0) THEN 2565 CALL QUIT('Out of memory in CC3_W3_OMEGA1 (OVVO-2)') 2566 ENDIF 2567C 2568 DO I = 1, LENGTH 2569 TMAT(I) = WMAT(I) 2570 ENDDO 2571C 2572 IF (NSYM .GT. 1) THEN 2573 IF (LWORK .LT. LENGTH) THEN 2574 CALL QUIT('Out of memory in CC3_W3_OMEGA1 (CC_GATHER-4)') 2575 ENDIF 2576 CALL CC_GATHER(LENGTH,WORK,TMAT,INDSQ(1,6)) 2577 CALL DCOPY(LENGTH,WORK,1,TMAT,1) 2578 ENDIF 2579C 2580C sort t2 amplitudes t2tp(elmD) as t2D(mle) 2581C 2582 DO ISYME = 1, NSYM 2583 ISYMLM = MULD2H(ISYELM,ISYME) 2584 DO ISYML = 1, NSYM 2585 ISYMM = MULD2H(ISYMLM,ISYML) 2586 ISYMEL = MULD2H(ISYME,ISYML) 2587C 2588 DO L = 1, NRHF(ISYML) 2589 DO M = 1, NRHF(ISYMM) 2590C 2591 KOFF1 = IT2SP(ISYELM,ISYMID) 2592 * + NCKI(ISYELM)*(ID-1) 2593 * + ICKI(ISYMEL,ISYMM) 2594 * + NT1AM(ISYMEL)*(M-1) 2595 * + IT1AM(ISYME,ISYML) 2596 * + NVIR(ISYME)*(L-1) 2597 * + 1 2598C 2599 KOFF2 = KSCR1 - 1 2600 * + IMAIJA(ISYMLM,ISYME) 2601 * + IMATIJ(ISYMM,ISYML) 2602 * + NRHF(ISYMM)*(L-1) 2603 * + M 2604C 2605 CALL DCOPY(NVIR(ISYME),T2TP(KOFF1),1, 2606 * WORK(KOFF2),NMATIJ(ISYMLM)) 2607C 2608 ENDDO 2609 ENDDO 2610 ENDDO 2611 ENDDO 2612C 2613 DO ISYME = 1, NSYM 2614C 2615 ISYMFN = MULD2H(ISYME,ISYENF) 2616 ISYMLM = MULD2H(ISCKIJ,ISYMFN) 2617 ISYFNI = MULD2H(ISYMI,ISYMFN) 2618 ISYLMI = MULD2H(ISYMI,ISYMLM) 2619C 2620 KSCR2 = KEND1 2621 KEND2 = KSCR2 + NMAIJK(ISYLMI) 2622 LWRK2 = LWORK - KEND2 2623C 2624 IF (LWRK2 .LT. 0) THEN 2625 CALL QUIT('Out of memory in CC3_W3_OMEGA1 (OVVO-3)') 2626 ENDIF 2627C 2628 DO E = 1, NVIR(ISYME) 2629C 2630 KOFF1 = ISAIKL(ISYMFN,ISYMLM) 2631 * + 1 2632 KOFF2 = IT2SP(ISYFNI,ISYME) 2633 * + NCKI(ISYFNI)*(E-1) 2634 * + ICKI(ISYMFN,ISYMI) 2635 * + 1 2636 KOFF3 = KSCR2 2637 * + IMAIJK(ISYMLM,ISYMI) 2638C 2639 NTOTFN = MAX(1,NT1AM(ISYMFN)) 2640 NTOTLM = MAX(1,NMATIJ(ISYMLM)) 2641C 2642C TMAT(fnlm) * xovvo(fnie) 2643C 2644C 2645 CALL DGEMM('T','N',NMATIJ(ISYMLM),NRHF(ISYMI), 2646 * NT1AM(ISYMFN), 2647 * -ONE,TMAT(KOFF1),NTOTFN,XOVVO(KOFF2),NTOTFN, 2648 * ZERO,WORK(KOFF3),NTOTLM) 2649C 2650 KOFF1 = KSCR2 2651 * + IMAIJK(ISYMLM,ISYMI) 2652 KOFF2 = KSCR1 2653 * + IMAIJA(ISYMLM,ISYME) 2654 * + NMATIJ(ISYMLM)*(E-1) 2655 KOFF3 = IT1AM(ISYMIB,ISYMI) 2656 * + IB 2657C 2658 NTOTB = MAX(1,NVIR(ISYMIB)) 2659 NTOTIJ = MAX(1,NMATIJ(ISYMLM)) 2660C 2661C omega1(Bi) = omega(Bi) - Tx(lmie) * t2D(lme) 2662C 2663C 2664C 3.2 omega contribution addomega 2665C 2666 CALL DGEMV('T',NMATIJ(ISYMLM),NRHF(ISYMI),-ONE, 2667 * WORK(KOFF1), 2668 * NTOTIJ,WORK(KOFF2),1,ONE,OMEGA1(KOFF3),NTOTB) 2669C 2670C 2671 ENDDO 2672 ENDDO 2673C 2674C---------------------------------------------- 2675C third contribution 2676C 2677C - Wdf(amnl) * t(dl,em) * g(iefn) 2678C 2679C - WBD(amnl) * t2tp(emlB) * xovvo(Dnie) 2680C 2681C - TMAT(amln) t2B(mle) * gD(eni) 2682C 2683C omega1(ai) = omega1(ai) - TMAT(amln) * t2g(mlni) 2684C---------------------------------------------- 2685C 2686 ISYMMLE = MULD2H(ISYMT2,ISYMIB) 2687 ISYMENI = MULD2H(ISYINT,ISYMID) 2688 ISYMLNI = MULD2H(ISYMMLE,ISYMENI) 2689 KGD = 1 2690 KT2B = KGD + NCKI(ISYMENI) 2691 KMLNI = KT2B + NCKI(ISYMMLE) 2692 KSCR1 = KMLNI + N3ORHF(ISYMLNI) 2693 KEND1 = KSCR1 + NCKIJ(ISWMAT) 2694 LWRK1 = LWORK - KEND1 2695C 2696 CALL DZERO(WORK(KMLNI),N3ORHF(ISYMLNI)) 2697C 2698 DO I = 1, LENGTH 2699 TMAT(I) = WMAT(INDSQ(I,3)) 2700 ENDDO 2701C 2702 IF (LWRK1 .LT. 0) THEN 2703 CALL QUIT('Out of memory in CC3_W3_OMEGA1 3.3') 2704 ENDIF 2705C 2706 IF (NSYM .GT. 1) THEN 2707 CALL CC3_SRTVOOO(WORK(KSCR1),TMAT,ISWMAT) 2708 CALL DCOPY(LENGTH,WORK(KSCR1),1,TMAT,1) 2709 END IF 2710C 2711 2712C 2713C---------------- 2714C Sort T2 tB(kmc) = t2tp(ckmB) 2715C---------------- 2716C 2717 ISYCKM = MULD2H(ISYMT2,ISYMIB) 2718 DO ISYMK = 1, NSYM 2719 ISYMCM = MULD2H(ISYCKM,ISYMK) 2720 DO ISYMC = 1, NSYM 2721 ISYMM = MULD2H(ISYMCM,ISYMC) 2722 ISYMKM = MULD2H(ISYMK,ISYMM) 2723 ISYMCK = MULD2H(ISYMK,ISYMC) 2724C 2725 DO K = 1, NRHF(ISYMK) 2726 DO M = 1, NRHF(ISYMM) 2727 KOFF1 = IT2SP(ISYCKM,ISYMIB) 2728 * + NCKI(ISYCKM)*(IB-1) 2729 * + ICKI(ISYMCK,ISYMM) 2730 * + NT1AM(ISYMCK)*(M-1) 2731 * + IT1AM(ISYMC,ISYMK) 2732 * + NVIR(ISYMC)*(K-1) 2733 * + 1 2734 KOFF2 = KT2B - 1 2735 * + IMAIJA(ISYMKM,ISYMC) 2736 * + IMATIJ(ISYMK,ISYMM) 2737 * + NRHF(ISYMK)*(M-1) 2738 * + K 2739C 2740 CALL DCOPY(NVIR(ISYMC),T2TP(KOFF1),1, 2741 * WORK(KOFF2),NMATIJ(ISYMKM)) 2742 ENDDO 2743 ENDDO 2744 ENDDO 2745 ENDDO 2746C 2747C--------------------------------------- 2748C sort xovvo(Dnie) as gD(eni) 2749C--------------------------------------- 2750C 2751 DO ISYME = 1,NSYM 2752 ISYMDNI = MULD2H(ISYINT,ISYME) 2753 ISYMNI = MULD2H(ISYMDNI,ISYMID) 2754 DO ISYMI = 1,NSYM 2755 ISYMN = MULD2H(ISYMNI,ISYMI) 2756 ISYMEN = MULD2H(ISYMENI,ISYMI) 2757 ISYMDN = MULD2H(ISYMDNI,ISYMI) 2758 DO E = 1,NVIR(ISYME) 2759 DO I = 1,NRHF(ISYMI) 2760 DO N = 1,NRHF(ISYMN) 2761 KOFF1 = IT2SP(ISYMDNI,ISYME) 2762 * + NCKI(ISYMDNI)*(E-1) 2763 * + ICKI(ISYMDN,ISYMI) 2764 * + NT1AM(ISYMDN)*(I-1) 2765 * + IT1AM(ISYMID,ISYMN) 2766 * + NVIR(ISYMID)*(N-1) 2767 * + ID 2768 KOFF2 = KGD - 1 2769 * + ICKI(ISYMEN,ISYMI) 2770 * + NT1AM(ISYMEN)*(I-1) 2771 * + IT1AM(ISYME,ISYMN) 2772 * + NVIR(ISYME)*(N-1) 2773 * + E 2774 WORK(KOFF2) = XOVVO(KOFF1) 2775 END DO 2776 END DO 2777 END DO 2778 END DO 2779 END DO 2780C 2781C - WBD(amnl) * t2tp(emlB) * xovvo(Dnie) 2782C 2783C - TMAT(amln) t2B(mle) * gD(eni) 2784C 2785 DO ISYMI = 1,NSYM 2786 ISYMMLN = MULD2H(ISYMLNI,ISYMI) 2787 ISYMEN = MULD2H(ISYMENI,ISYMI) 2788 DO ISYMN = 1,NSYM 2789 ISYME = MULD2H(ISYMEN,ISYMN) 2790 ISYMML = MULD2H(ISYMMLE,ISYME) 2791 DO I = 1,NRHF(ISYMI) 2792C 2793 KOFF1 = KT2B 2794 * + IMAIJA(ISYMML,ISYME) 2795 KOFF2 = KGD 2796 * + ICKI(ISYMEN,ISYMI) 2797 * + NT1AM(ISYMEN)*(I-1) 2798 * + IT1AM(ISYME,ISYMN) 2799 2800 KOFF3 = KMLNI 2801 * + I3ORHF(ISYMMLN,ISYMI) 2802 * + NMAIJK(ISYMMLN)*(I-1) 2803 * + IMAIJK(ISYMML,ISYMN) 2804C 2805 NTOTML = MAX(1,NMATIJ(ISYMML)) 2806 NTOTE = MAX(1,NVIR(ISYME)) 2807C 2808C 2809 CALL DGEMM('N','N',NMATIJ(ISYMML),NRHF(ISYMN), 2810 * NVIR(ISYME), 2811 * ONE,WORK(KOFF1),NTOTML,WORK(KOFF2),NTOTE, 2812 * ONE,WORK(KOFF3),NTOTML) 2813C 2814 2815 END DO 2816 END DO 2817 END DO 2818C 2819C omega1(ai) = omega1(ai) - TMAT(amln) * t2Bgd(mlni) 2820C 2821C t2B(mle) * gD(eni) 2822C 2823 DO ISYMI = 1,NSYM 2824 ISYMMLN = MULD2H(ISYMLNI,ISYMI) 2825 ISYMA = MULD2H(ISWMAT,ISYMMLN) 2826C 2827 KOFF1 = I3VOOO(ISYMA,ISYMMLN) + 1 2828 KOFF2 = KMLNI + I3ORHF(ISYMMLN,ISYMI) 2829 KOFF3 = IT1AM(ISYMA,ISYMI) + 1 2830C 2831 NTOTA = MAX(1,NVIR(ISYMA)) 2832 NTOTMLN = MAX(1,NMAIJK(ISYMMLN)) 2833C 2834C 3.3 omega contribution addomega 2835C 2836 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI),NMAIJK(ISYMMLN), 2837 * ONE,TMAT(KOFF1),NTOTA,WORK(KOFF2),NTOTMLN, 2838 * ONE,OMEGA1(KOFF3),NTOTA) 2839 END DO 2840C 2841 END IF 2842C 2843C================================================ 2844C 4. Calculate contributions from g_{oovv} 2845C - L^{daf}_{lnm} t^{de}_{lm} g_{infe} = 2846C -( Waf(dlmn) + Wad(fmln) + Wfd(anlm) )* t(dl,em) * g(infe) 2847C================================================ 2848C 2849C-------------------------------- 2850C First contribution 2851C-------------------------------- 2852C - Waf(dlmn) * t(dl,em) * g(infe) 2853C 2854C t2T(en) = t2tp(glme) * TMAT(glmn) 2855C 2856C g_{inDe} = xoovv(Dine) = gD(eni) 2857C 2858C wmega1(Ai) = wmega1(Ai) + gD(eni) * t2T(en) 2859C 2860 ISYMEN = MULD2H(ISCKIJ,ISYMT2) 2861 ISYMI = MULD2H(ISYRES,ISYMIB) 2862 ISYENI = MULD2H(ISYMEN,ISYMI) 2863C 2864 KSCR1 = 1 2865 KSCR2 = KSCR1 + NT1AM(ISYMEN) 2866 KEND1 = KSCR2 + NCKI(ISYENI) 2867 LWRK1 = LWORK - KEND1 2868C 2869 IF (LWRK1 .LT. 0) THEN 2870 CALL QUIT('Out of memory in CC3_W3_OMEGA1 (OOVV-1)') 2871 ENDIF 2872C 2873 CALL DZERO(WORK(KSCR1),NT1AM(ISYMEN)) 2874 CALL DZERO(WORK(KSCR2),NCKI(ISYENI)) 2875C 2876C 2877C t2T(en) = t2tp(glme) * TMAT(glmn) 2878C 2879 DO I = 1, LENGTH 2880 TMAT(I) = WMAT(I) 2881 ENDDO 2882C 2883 DO ISYME = 1, NSYM 2884 ISYDLM = MULD2H(ISYME,ISYMT2) 2885 ISYMN = MULD2H(ISYMEN,ISYME) 2886C 2887 KOFF1 = IT2SP(ISYDLM,ISYME) 2888 * + 1 2889 KOFF2 = ISAIKJ(ISYDLM,ISYMN) 2890 * + 1 2891 KOFF3 = KSCR1 2892 * + IT1AM(ISYME,ISYMN) 2893C 2894 NTODLM = MAX(1,NCKI(ISYDLM)) 2895 NTOTE = MAX(1,NVIR(ISYME)) 2896C 2897 CALL DGEMM('T','N',NVIR(ISYME),NRHF(ISYMN),NCKI(ISYDLM), 2898 * -ONE,T2TP(KOFF1),NTODLM,TMAT(KOFF2),NTODLM, 2899 * ONE,WORK(KOFF3),NTOTE) 2900C 2901 ENDDO 2902C 2903C g_{inDe} = xoovv(Dine) = gD(eni) 2904C 2905 DO ISYME = 1, NSYM 2906 ISYMN = MULD2H(ISYMEN,ISYME) 2907 ISYMEI = MULD2H(ISYME,ISYMI) 2908 ISYMDN = MULD2H(ISYMN,ISYMID) 2909 ISYMDI = MULD2H(ISYMI,ISYMID) 2910 ISYDNI = MULD2H(ISYMDN,ISYMI) 2911C 2912 DO E = 1, NVIR(ISYME) 2913 DO N = 1, NRHF(ISYMN) 2914C 2915 KOFF1 = IT2SP(ISYDNI,ISYME) 2916 * + NCKI(ISYDNI)*(E-1) 2917 * + ICKI(ISYMDI,ISYMN) 2918 * + NT1AM(ISYMDI)*(N-1) 2919 * + IT1AM(ISYMID,ISYMI) 2920 * + ID 2921C 2922 KOFF2 = KSCR2 - 1 2923 * + ICKI(ISYMEN,ISYMI) 2924 * + IT1AM(ISYME,ISYMN) 2925 * + NVIR(ISYME)*(N-1) 2926 * + E 2927CC 2928 CALL DCOPY(NRHF(ISYMI),XOOVV(KOFF1),NVIR(ISYMID), 2929 * WORK(KOFF2),NT1AM(ISYMEN)) 2930C 2931 ENDDO 2932 ENDDO 2933 ENDDO 2934C 2935C wmega1(Ai) = wmega1(Ai) + gD(eni) * t2T(en) 2936C 2937 KOFF1 = KSCR2 2938 * + ICKI(ISYMEN,ISYMI) 2939 KOFF3 = IT1AM(ISYMIB,ISYMI) 2940 * + IB 2941C 2942 NTOTEN = MAX(1,NT1AM(ISYMEN)) 2943 NTOTB = MAX(1,NVIR(ISYMIB)) 2944C 2945C 4.1 omega contribution addomega 2946C 2947 CALL DGEMV('T',NT1AM(ISYMEN),NRHF(ISYMI),-ONE,WORK(KOFF1), 2948 * NTOTEN,WORK(KSCR1),1,ONE,OMEGA1(KOFF3),NTOTB) 2949C 2950C-------------------------------- 2951C Second contribution 2952C-------------------------------- 2953C 2954C - Wad(fmln) * t(dl,em) * g(infe) 2955C 2956C g_{infe} = xoovv(fine) = ge(fni) 2957C 2958C t2tp(elmD) = tD(mle) 2959C 2960C TX(lmie) = TMAT(fnlm) * ge(fni) tD(mle) 2961C 2962C wmega1(Ai) = wmega1(Ai) - TX(lmie) * tD(mle) 2963C 2964 ! Calculate only when contracting with first-order multipliers 2965 IF (W3X) THEN 2966C 2967 2968 ISYTMP = MULD2H(ISCKIJ,ISYINT) 2969 ISYMI = MULD2H(ISYRES,ISYMIB) 2970 ISYENF = MULD2H(ISYINT,ISYMI) 2971 ISYELM = MULD2H(ISYMT2,ISYMID) 2972C 2973 KSCR1 = 1 2974 KEND1 = KSCR1 + NMAIJA(ISYELM) 2975 LWRK1 = LWORK - KEND1 2976C 2977 IF (LWRK1 .LT. 0) THEN 2978 CALL QUIT('Out of memory in CC3_W3_OMEGA1 (OOVV-2)') 2979 ENDIF 2980C 2981 DO I = 1, LENGTH 2982 TMAT(I) = WMAT(INDSQ(I,5)) 2983 ENDDO 2984C 2985 IF (NSYM .GT. 1) THEN 2986 IF (LWORK .LT. LENGTH) THEN 2987 CALL QUIT('Out of memory in CC3_W3_OMEGA1 (CC_GATHER-5)') 2988 ENDIF 2989 CALL CC_GATHER(LENGTH,WORK,TMAT,INDSQ(1,6)) 2990 CALL DCOPY(LENGTH,WORK,1,TMAT,1) 2991 ENDIF 2992C 2993C t2tp(elmD) = tD(mle) 2994C 2995 DO ISYME = 1, NSYM 2996 ISYMLM = MULD2H(ISYELM,ISYME) 2997 DO ISYML = 1, NSYM 2998 ISYMM = MULD2H(ISYMLM,ISYML) 2999 ISYMEL = MULD2H(ISYME,ISYML) 3000C 3001 DO L = 1, NRHF(ISYML) 3002 DO M = 1, NRHF(ISYMM) 3003C 3004 KOFF1 = IT2SP(ISYELM,ISYMID) 3005 * + NCKI(ISYELM)*(ID-1) 3006 * + ICKI(ISYMEL,ISYMM) 3007 * + NT1AM(ISYMEL)*(M-1) 3008 * + IT1AM(ISYME,ISYML) 3009 * + NVIR(ISYME)*(L-1) 3010 * + 1 3011C 3012 KOFF2 = KSCR1 - 1 3013 * + IMAIJA(ISYMLM,ISYME) 3014 * + IMATIJ(ISYMM,ISYML) 3015 * + NRHF(ISYMM)*(L-1) 3016 * + M 3017C 3018 CALL DCOPY(NVIR(ISYME),T2TP(KOFF1),1, 3019 * WORK(KOFF2),NMATIJ(ISYMLM)) 3020C 3021 ENDDO 3022 ENDDO 3023 ENDDO 3024 ENDDO 3025C 3026 DO ISYME = 1, NSYM 3027C 3028 ISYMFN = MULD2H(ISYME,ISYENF) 3029 ISYMLM = MULD2H(ISCKIJ,ISYMFN) 3030 ISYFNI = MULD2H(ISYMI,ISYMFN) 3031 ISYLMI = MULD2H(ISYMI,ISYMLM) 3032C 3033 KSCR2 = KEND1 3034 KSCR3 = KSCR2 + NMAIJK(ISYLMI) 3035 KEND2 = KSCR3 + NCKI(ISYFNI) 3036 LWRK2 = LWORK - KEND2 3037C 3038 IF (LWRK2 .LT. 0) THEN 3039 CALL QUIT('Out of memory in CC3_W3_OMEGA1 (OOVV-3)') 3040 ENDIF 3041C 3042C g_{infe} = xoovv(fine) = ge(fni) 3043C 3044 DO E = 1, NVIR(ISYME) 3045C 3046 DO ISYMF = 1, NSYM 3047 ISYMN = MULD2H(ISYMFN,ISYMF) 3048 ISYMFI = MULD2H(ISYMI,ISYMF) 3049 DO F = 1, NVIR(ISYMF) 3050 DO N = 1, NRHF(ISYMN) 3051C 3052 KOFF1 = IT2SP(ISYFNI,ISYME) 3053 * + NCKI(ISYFNI)*(E-1) 3054 * + ICKI(ISYMFI,ISYMN) 3055 * + NT1AM(ISYMFI)*(N-1) 3056 * + IT1AM(ISYMF,ISYMI) 3057 * + F 3058C 3059 KOFF2 = KSCR3 - 1 3060 * + ICKI(ISYMFN,ISYMI) 3061 * + IT1AM(ISYMF,ISYMN) 3062 * + NVIR(ISYMF)*(N-1) 3063 * + F 3064C 3065 CALL DCOPY(NRHF(ISYMI),XOOVV(KOFF1),NVIR(ISYMF), 3066 * WORK(KOFF2),NT1AM(ISYMFN)) 3067 ENDDO 3068 ENDDO 3069 ENDDO 3070C 3071 KOFF1 = ISAIKL(ISYMFN,ISYMLM) 3072 * + 1 3073 KOFF2 = KSCR3 3074 * + ICKI(ISYMFN,ISYMI) 3075 KOFF3 = KSCR2 3076 * + IMAIJK(ISYMLM,ISYMI) 3077C 3078 NTOTFN = MAX(1,NT1AM(ISYMFN)) 3079 NTOTLM = MAX(1,NMATIJ(ISYMLM)) 3080C 3081C TX(lmie) = TMAT(fnlm) * ge(fni) 3082C 3083 CALL DGEMM('T','N',NMATIJ(ISYMLM),NRHF(ISYMI), 3084 * NT1AM(ISYMFN), 3085 * -ONE,TMAT(KOFF1),NTOTFN,WORK(KOFF2),NTOTFN, 3086 * ZERO,WORK(KOFF3),NTOTLM) 3087C 3088 KOFF1 = KSCR2 3089 * + IMAIJK(ISYMLM,ISYMI) 3090 KOFF2 = KSCR1 3091 * + IMAIJA(ISYMLM,ISYME) 3092 * + NMATIJ(ISYMLM)*(E-1) 3093 KOFF3 = IT1AM(ISYMIB,ISYMI) 3094 * + IB 3095C 3096 NTOTB = MAX(1,NVIR(ISYMIB)) 3097 NTOTIJ = MAX(1,NMATIJ(ISYMLM)) 3098C 3099C wmega1(Ai) = wmega1(Ai) - TX(lmie) * tD(mle) 3100C 3101C 3102C 4.2 omega contribution addomega 3103C 3104 CALL DGEMV('T',NMATIJ(ISYMLM),NRHF(ISYMI),-ONE, 3105 * WORK(KOFF1), 3106 * NTOTIJ,WORK(KOFF2),1,ONE,OMEGA1(KOFF3),NTOTB) 3107C 3108 3109 ENDDO 3110 ENDDO 3111C 3112C-------------------------------- 3113C Third contribution 3114C-------------------------------- 3115C 3116C - Wfd(anlm) * t(dl,em) * g(infe) 3117C 3118C - WBD(anlm) * t2tp(emlD) * xoovv(Bine) 3119C 3120C tD(mle) * xB(eni) = tdxB(mlni) 3121C 3122C omega1(ai) = omega1(ai) + TMAT(amln) * tdxB(mlni) 3123C 3124C 3125 ISYMMLE = MULD2H(ISYMT2,ISYMID) 3126 ISYMENI = MULD2H(ISYINT,ISYMIB) 3127 ISYMLNI = MULD2H(ISYMMLE,ISYMENI) 3128C 3129 KSCR1 = 1 3130 KSCR2 = KSCR1 + NMAIJA(ISYMMLE) 3131 KSCR3 = KSCR2 + NMAIJA(ISYMENI) 3132 KSCR4 = KSCR3 + N3ORHF(ISYMLNI) 3133 KEND1 = KSCR4 + NCKIJ(ISWMAT) 3134 LWRK1 = LWORK - KEND1 3135C 3136 CALL DZERO(WORK(KSCR3),N3ORHF(ISYMLNI)) 3137 IF (LWRK1 .LT. 0) THEN 3138 CALL QUIT('Out of memory in CC3_W3_OMEGA1 4.3') 3139 ENDIF 3140C 3141 DO I = 1, LENGTH 3142 TMAT(I) = WMAT(INDSQ(I,5)) 3143 ENDDO 3144C 3145 IF (LWORK .LT. NCKIJ(ISCKIJ)) THEN 3146 CALL QUIT('Out of memory in CC3_W3_OMEGA1 4.3.1') 3147 ENDIF 3148C 3149 IF (NSYM .GT. 1) THEN 3150 CALL CC3_SRTVOOO(WORK(KSCR4),TMAT,ISWMAT) 3151 CALL DCOPY(LENGTH,WORK(KSCR4),1,TMAT,1) 3152 END IF 3153C 3154C 3155C t2tp(emlD) = tD(mle) 3156C 3157 ISYMEML = MULD2H(ISYMT2,ISYMID) 3158 DO ISYME = 1, NSYM 3159 ISYMML = MULD2H(ISYMEML,ISYME) 3160 DO ISYML = 1, NSYM 3161 ISYMM = MULD2H(ISYMML,ISYML) 3162 ISYMEM = MULD2H(ISYME,ISYMM) 3163C 3164 DO E = 1,NVIR(ISYME) 3165 DO L = 1, NRHF(ISYML) 3166 DO M = 1, NRHF(ISYMM) 3167C 3168 KOFF1 = IT2SP(ISYMEML,ISYMID) 3169 * + NCKI(ISYMEML)*(ID-1) 3170 * + ICKI(ISYMEM,ISYML) 3171 * + NT1AM(ISYMEM)*(L-1) 3172 * + IT1AM(ISYME,ISYMM) 3173 * + NVIR(ISYME)*(M-1) 3174 * + E 3175C 3176 KOFF2 = KSCR1 - 1 3177 * + IMAIJA(ISYMML,ISYME) 3178 * + NMATIJ(ISYMML)*(E-1) 3179 * + IMATIJ(ISYMM,ISYML) 3180 * + NRHF(ISYMM)*(L-1) 3181 * + M 3182C 3183 WORK(KOFF2) = T2TP(KOFF1) 3184C 3185 ENDDO 3186 ENDDO 3187 ENDDO 3188 ENDDO 3189 ENDDO 3190C 3191C 3192C xoovv(Bine) = xB(eni) 3193C 3194C 3195 DO ISYME = 1, NSYM 3196 ISYMBIN = MULD2H(ISYINT,ISYME) 3197 DO ISYMN = 1,NSYM 3198 ISYMEN = MULD2H(ISYME,ISYMN) 3199 ISYMBI = MULD2H(ISYMBIN,ISYMN) 3200 ISYMI = MULD2H(ISYMBI,ISYMIB) 3201 DO E = 1, NVIR(ISYME) 3202 DO N = 1, NRHF(ISYMN) 3203 DO I = 1, NRHF(ISYMI) 3204C 3205 KOFF1 = IT2SP(ISYMBIN,ISYME) 3206 * + NCKI(ISYMBIN)*(E-1) 3207 * + ICKI(ISYMBI,ISYMN) 3208 * + NT1AM(ISYMBI)*(N-1) 3209 * + IT1AM(ISYMIB,ISYMI) 3210 * + NVIR(ISYMIB)*(I-1) 3211 * + IB 3212C 3213 KOFF2 = KSCR2 - 1 3214 * + ICKI(ISYMEN,ISYMI) 3215 * + NT1AM(ISYMEN)*(I-1) 3216 * + IT1AM(ISYME,ISYMN) 3217 * + NVIR(ISYME)*(N-1) 3218 * + E 3219C 3220 WORK(KOFF2) = XOOVV(KOFF1) 3221 ENDDO 3222 ENDDO 3223 ENDDO 3224 ENDDO 3225 ENDDO 3226C 3227C tD(mle) * xB(eni) = tdxB(mlni) 3228C 3229 DO ISYMI = 1,NSYM 3230 ISYMEN = MULD2H(ISYMENI,ISYMI) 3231 DO ISYMN = 1,NSYM 3232 ISYME = MULD2H(ISYMEN,ISYMN) 3233 ISYMML = MULD2H(ISYMMLE,ISYME) 3234 ISYMMLN = MULD2H(ISYMML,ISYMN) 3235 DO I = 1,NRHF(ISYMI) 3236 3237 KOFF1 = KSCR1 3238 * + IMAIJA(ISYMML,ISYME) 3239 KOFF2 = KSCR2 3240 * + ICKI(ISYMEN,ISYMI) 3241 * + NT1AM(ISYMEN)*(I-1) 3242 * + IT1AM(ISYME,ISYMN) 3243 KOFF3 = KSCR3 3244 * + I3ORHF(ISYMMLN,ISYMI) 3245 * + NMAIJK(ISYMMLN) * (I-1) 3246 * + IMAIJK(ISYMML,ISYMN) 3247C 3248 NTOTML = MAX(1,NMATIJ(ISYMML)) 3249 NTOTE = MAX(1,NVIR(ISYME)) 3250C 3251 CALL DGEMM('N','N',NMATIJ(ISYMML),NRHF(ISYMN), 3252 * NVIR(ISYME),-ONE,WORK(KOFF1), 3253 * NTOTML,WORK(KOFF2),NTOTE, 3254 * ONE,WORK(KOFF3),NTOTML) 3255C 3256 END DO 3257 END DO 3258 END DO 3259C 3260C omega1(ai) = omega1(ai) + TMAT(amln) * tdxB(mlni) 3261C 3262 DO ISYMI = 1,NSYM 3263 ISYMMLN = MULD2H(ISYMLNI,ISYMI) 3264 ISYMA = MULD2H(ISYRES,ISYMI) 3265C 3266 KOFF1 = I3VOOO(ISYMA,ISYMMLN) + 1 3267 KOFF2 = KSCR3 + I3ORHF(ISYMMLN,ISYMI) 3268 KOFF3 = IT1AM(ISYMA,ISYMI) + 1 3269C 3270 NTOTMLN = MAX(1,NMAIJK(ISYMMLN)) 3271 NTOTA = MAX(1,NVIR(ISYMA)) 3272C 3273C 4.3 omega contribution addomega 3274C 3275 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI), 3276 * NMAIJK(ISYMMLN),-ONE,TMAT(KOFF1), 3277 * NTOTA,WORK(KOFF2),NTOTMLN, 3278 * ONE,OMEGA1(KOFF3),NTOTA) 3279C 3280 END DO 3281C 3282 END IF 3283C 3284C================================================ 3285C 5. Calculate contributions from g_{ovvo} 3286C - L^{def}_{lin} t^{de}_{lm} g_{mafn} 3287C -( Wfe(dlin) + Wfd(eiln) + Wde(fnil) )* t(dl,em) * g(mafn) 3288C================================================ 3289C 3290 ISYDLM = MULD2H(ISYMT2,ISYMID) 3291 ISYTMP = MULD2H(ISYDLM,ISCKIJ) 3292C 3293 KSCR1 = 1 3294 KSCR2 = KSCR1 + NMAIJK(ISYTMP) 3295 KEND1 = KSCR2 + NCKI(ISYDLM) 3296 LWRK1 = LWORK - KEND1 3297C 3298 IF (LWRK1 .LT. 0) THEN 3299 CALL QUIT('Out of memory CC3_W3_OMEGA1 (OVVO-4)') 3300 ENDIF 3301C 3302 CALL DZERO(WORK(KSCR1),NMAIJK(ISYTMP)) 3303C 3304C---------------------------------------------- 3305C -( Wfe(dlin) + Wfd(eiln) )* t(dl,em) * g(mafn) 3306C---------------------------------------------- 3307C 3308C - Wfe(dlin) * t(dl,em) * g(mafn) 3309C 3310C Tt2(inm) = TMAT(dlin) * t2tp(dlmD) 3311C 3312C - Wfd(eiln) * t(dl,em) * g(mafn) 3313C 3314C Tt2(inm) = Tt2(inm) + TMAT(elin) * t2tp(emlD) 3315C tD(elm) 3316C g(maBn) = xovvo(Bnma) = IB(mna) 3317C 3318C Tt2(inm) = T(mni) 3319C 3320C omega(ai) = omega(ai) + IB(mna) * T(mni) 3321C 3322C----------------------------------------------- 3323C First contribution to intermediate 3324C----------------------------------------------- 3325C 3326C - Wfe(dlin) * t(dl,em) * g(mafn) 3327 DO I = 1, LENGTH 3328 TMAT(I) = WMAT(I) 3329 ENDDO 3330C 3331 IF (NSYM .GT. 1) THEN 3332 IF (LWRK1 .LT. LENGTH) THEN 3333 CALL QUIT('Out of memory in CC3_W3_OMEGA1 (CC_GATHER-6)') 3334 ENDIF 3335 CALL CC_GATHER(LENGTH,WORK(KEND1),TMAT,INDSQ(1,6)) 3336 CALL DCOPY(LENGTH,WORK(KEND1),1,TMAT,1) 3337 ENDIF 3338C 3339C Tt2(inm) = TMAT(dlin) * t2tp(dlmD) 3340C 3341 DO ISYMDL = 1, NSYM 3342 ISYMM = MULD2H(ISYDLM,ISYMDL) 3343 ISYMIN = MULD2H(ISCKIJ,ISYMDL) 3344C 3345 KOFF1 = ISAIKL(ISYMDL,ISYMIN) 3346 * + 1 3347 KOFF2 = IT2SP(ISYDLM,ISYMID) 3348 * + NCKI(ISYDLM)*(ID-1) 3349 * + ICKI(ISYMDL,ISYMM) 3350 * + 1 3351 KOFF3 = KSCR1 3352 * + IMAIJK(ISYMIN,ISYMM) 3353C 3354 NTOTDL = MAX(1,NT1AM(ISYMDL)) 3355 NTOTIN = MAX(1,NMATIJ(ISYMIN)) 3356C 3357 CALL DGEMM('T','N',NMATIJ(ISYMIN),NRHF(ISYMM),NT1AM(ISYMDL), 3358 * -ONE,TMAT(KOFF1),NTOTDL,T2TP(KOFF2),NTOTDL, 3359 * ONE,WORK(KOFF3),NTOTIN) 3360 ENDDO 3361c 3362* write(lupri,*)'kscr1(on,m) in cc3_w3_lhtr ib,id ',ib,id 3363* call output(work(kscr1),1,nmatij(1),1,nrhf(1), 3364* * nmatij(1),nrhf(1), 3365* * 1,lupri) 3366 3367C 3368C----------------------------------------------- 3369C Second contribution to intermediate 3370C----------------------------------------------- 3371C 3372C Tt2(inm) = Tt2(inm) + TMAT(elin) * t2tp(emlD) 3373C tD(elm) 3374C 3375 ! Calculate only when contracting with first-order multipliers 3376 IF (W3X) THEN 3377C 3378 DO I = 1, LENGTH 3379 TMAT(I) = WMAT(INDSQ(I,1)) 3380 ENDDO 3381C 3382 IF (NSYM .GT. 1) THEN 3383 IF (LWRK1 .LT. LENGTH) THEN 3384 CALL QUIT('Out of memory in CC3_W3_OMEGA1 (CC_GATHER-6)') 3385 ENDIF 3386 CALL CC_GATHER(LENGTH,WORK(KEND1),TMAT,INDSQ(1,6)) 3387 CALL DCOPY(LENGTH,WORK(KEND1),1,TMAT,1) 3388 ENDIF 3389C 3390C t2tp(emlD) = tD(elm) 3391C 3392 DO ISYMM = 1, NSYM 3393 ISYMDL = MULD2H(ISYDLM,ISYMM) 3394 DO ISYMD = 1, NSYM 3395 ISYML = MULD2H(ISYMDL,ISYMD) 3396 ISYMDM = MULD2H(ISYMD,ISYMM) 3397 DO M = 1, NRHF(ISYMM) 3398 DO L = 1, NRHF(ISYML) 3399C 3400 KOFF1 = IT2SP(ISYDLM,ISYMID) 3401 * + NCKI(ISYDLM)*(ID-1) 3402 * + ICKI(ISYMDL,ISYMM) 3403 * + NT1AM(ISYMDL)*(M-1) 3404 * + IT1AM(ISYMD,ISYML) 3405 * + NVIR(ISYMD)*(L-1) 3406 * + 1 3407C 3408 KOFF2 = KSCR2 3409 * + ICKI(ISYMDM,ISYML) 3410 * + NT1AM(ISYMDM)*(L-1) 3411 * + IT1AM(ISYMD,ISYMM) 3412 * + NVIR(ISYMD)*(M-1) 3413C 3414 CALL DCOPY(NVIR(ISYMD),T2TP(KOFF1),1, 3415 * WORK(KOFF2),1) 3416C 3417 ENDDO 3418 ENDDO 3419 ENDDO 3420 ENDDO 3421C 3422C Tt2(inm) = Tt2(inm) + TMAT(elin) * t2tp(emlD) 3423C tD(elm) 3424 DO ISYMDL = 1, NSYM 3425 ISYMM = MULD2H(ISYDLM,ISYMDL) 3426 ISYMIN = MULD2H(ISCKIJ,ISYMDL) 3427C 3428 KOFF1 = ISAIKL(ISYMDL,ISYMIN) 3429 * + 1 3430 KOFF2 = KSCR2 3431 * + ICKI(ISYMDL,ISYMM) 3432 KOFF3 = KSCR1 3433 * + IMAIJK(ISYMIN,ISYMM) 3434C 3435 NTOTDL = MAX(1,NT1AM(ISYMDL)) 3436 NTOTIN = MAX(1,NMATIJ(ISYMIN)) 3437C 3438 CALL DGEMM('T','N',NMATIJ(ISYMIN),NRHF(ISYMM),NT1AM(ISYMDL), 3439 * -ONE,TMAT(KOFF1),NTOTDL,WORK(KOFF2),NTOTDL, 3440 * ONE,WORK(KOFF3),NTOTIN) 3441 ENDDO 3442c 3443* write(lupri,*)'kscr1(on,m) in cc3_w3_lhtr ib,id ',ib,id 3444* call output(work(kscr1),1,nmatij(1),1,nrhf(1), 3445* * nmatij(1),nrhf(1), 3446* * 1,lupri) 3447 3448C 3449 END IF 3450C 3451C------------------------------------------------------- 3452C Sort intermediate and integrals and contract 3453C------------------------------------------------------- 3454C 3455 KEND1 = KSCR2 + NMAIJK(ISYTMP) 3456 LWRK1 = LWORK - KEND1 3457C 3458 IF (LWRK1 .LT. 0) THEN 3459 CALL QUIT('Out of memory CC3_W3_OMEGA1 (OVVO-5)') 3460 ENDIF 3461C 3462C Tt2(inm) = T(mni) 3463C 3464 DO ISYMM = 1, NSYM 3465 ISYMIN = MULD2H(ISYTMP,ISYMM) 3466 DO ISYMN = 1, NSYM 3467 ISYMI = MULD2H(ISYMIN,ISYMN) 3468 ISYMMN = MULD2H(ISYMM,ISYMN) 3469C 3470 DO M = 1, NRHF(ISYMM) 3471 DO N = 1, NRHF(ISYMN) 3472C 3473 KOFF1 = KSCR1 3474 * + IMAIJK(ISYMIN,ISYMM) 3475 * + NMATIJ(ISYMIN)*(M-1) 3476 * + IMATIJ(ISYMI,ISYMN) 3477 * + NRHF(ISYMI)*(N-1) 3478C 3479 KOFF2 = KSCR2 - 1 3480 * + IMAIJK(ISYMMN,ISYMI) 3481 * + IMATIJ(ISYMM,ISYMN) 3482 * + NRHF(ISYMM)*(N-1) 3483 * + M 3484C 3485 CALL DCOPY(NRHF(ISYMI),WORK(KOFF1),1, 3486 * WORK(KOFF2),NMATIJ(ISYMMN)) 3487C 3488 ENDDO 3489 ENDDO 3490 ENDDO 3491 ENDDO 3492C 3493 CALL DCOPY(NMAIJK(ISYTMP),WORK(KSCR2),1,WORK(KSCR1),1) 3494C 3495 ISYAMN = MULD2H(ISYINT,ISYMIB) 3496C 3497C g(maBn) = xovvo(Bnma) = IB(mna) 3498C 3499 DO ISYMA = 1, NSYM 3500 ISYMMN = MULD2H(ISYMA,ISYAMN) 3501 ISYBMN = MULD2H(ISYINT,ISYMA) 3502 DO ISYMM = 1, NSYM 3503 ISYMN = MULD2H(ISYMMN,ISYMM) 3504 ISYMBN = MULD2H(ISYBMN,ISYMM) 3505C 3506 DO M = 1, NRHF(ISYMM) 3507 DO N = 1, NRHF(ISYMN) 3508C 3509 KOFF1 = IT2SP(ISYBMN,ISYMA) 3510 * + ICKI(ISYMBN,ISYMM) 3511 * + NT1AM(ISYMBN)*(M-1) 3512 * + IT1AM(ISYMIB,ISYMN) 3513 * + NVIR(ISYMIB)*(N-1) 3514 * + IB 3515C 3516 KOFF2 = KSCR2 - 1 3517 * + IMAIJA(ISYMMN,ISYMA) 3518 * + IMATIJ(ISYMM,ISYMN) 3519 * + NRHF(ISYMM)*(N-1) 3520 * + M 3521C 3522 CALL DCOPY(NVIR(ISYMA),XOVVO(KOFF1),NCKI(ISYBMN), 3523 * WORK(KOFF2),NMATIJ(ISYMMN)) 3524C 3525 ENDDO 3526 ENDDO 3527 ENDDO 3528 ENDDO 3529C 3530C omega(ai) = omega(ai) + IB(mna) * T(mni) 3531C 3532 DO ISYMA = 1, NSYM 3533 ISYMI = MULD2H(ISYRES,ISYMA) 3534 ISYMMN = MULD2H(ISYAMN,ISYMA) 3535C 3536 KOFF1 = KSCR2 3537 * + IMAIJA(ISYMMN,ISYMA) 3538 KOFF2 = KSCR1 3539 * + IMAIJK(ISYMMN,ISYMI) 3540 KOFF3 = IT1AM(ISYMA,ISYMI) 3541 * + 1 3542C 3543 NTOTMN = MAX(1,NMATIJ(ISYMMN)) 3544 NTOTA = MAX(1,NVIR(ISYMA)) 3545C 3546C 5.1 + 5.2 omega contribution addomega 3547C 3548 CALL DGEMM('T','N',NVIR(ISYMA),NRHF(ISYMI),NMATIJ(ISYMMN), 3549 * -ONE,WORK(KOFF1),NTOTMN,WORK(KOFF2),NTOTMN, 3550 * ONE,OMEGA1(KOFF3),NTOTA) 3551 ENDDO 3552C 3553C------------------------------------- 3554C Third contribution 3555C 3556C - Wde(fnil) * t(dl,em) * g(mafn) 3557C------------------------------------- 3558C 3559C TMAT(fnil) * t2tp(BlmD xovvo(fnma) 3560C 3561C Tt2(fnim) * xovvo(fnma) 3562C tt2(fnmi) 3563C 3564C omega(ai) = omega(ai) + xoovv(fnma) * tt2(fnmi) 3565C 3566C 3567 ! Calculate only when contracting with first-order multipliers 3568 IF (W3X) THEN 3569C 3570 3571 ISYMBLM = MULD2H(ISYMT2,ISYMID) 3572 ISYMLM = MULD2H(ISYMBLM,ISYMIB) 3573 ISYFNIM = MULD2H(ISYMLM,ISCKIJ) 3574C 3575 KSCR1 = 1 3576 KSCR2 = KSCR1 + NMATIJ(ISYMLM) 3577 KEND1 = KSCR2 + NCKIJ(ISYFNIM) 3578 LWRK1 = LWORK - KEND1 3579C 3580 IF (LWRK1 .LT. 0) THEN 3581 CALL QUIT('Out of memory CC3_W3_OMEG 6.)') 3582 ENDIF 3583C 3584C 3585 CALL DZERO(WORK(KSCR2),NCKIJ(ISYFNIM)) 3586C 3587 DO I = 1, LENGTH 3588 TMAT(I) = WMAT(I) 3589 ENDDO 3590 3591C 3592C ------------------------------- 3593C Sort T^{BD}_{mi} as T_{mi} 3594C ------------------------------- 3595C 3596 ISYMBD = MULD2H(ISYMIB,ISYMID) 3597 ISYMMI = MULD2H(ISYMBD,ISYMT2) 3598 ISYBMI = MULD2H(ISYMMI,ISYMIB) 3599C 3600 DO ISYMI = 1,NSYM 3601 ISYMM = MULD2H(ISYMMI,ISYMI) 3602 ISYMBM = MULD2H(ISYBMI,ISYMI) 3603 DO I = 1,NRHF(ISYMI) 3604 DO M = 1,NRHF(ISYMM) 3605 KOFF1 = IT2SP(ISYBMI,ISYMID) 3606 * + NCKI(ISYBMI)*(ID-1) 3607 * + ICKI(ISYMBM,ISYMI) 3608 * + NT1AM(ISYMBM)*(I-1) 3609 * + IT1AM(ISYMIB,ISYMM) 3610 * + NVIR(ISYMIB)*(M-1) 3611 * + IB 3612 KOFF2 = IMATIJ(ISYMM,ISYMI) 3613 * + NRHF(ISYMM)*(I-1) 3614 * + M 3615 WORK(KSCR1-1+KOFF2) = T2TP(KOFF1) 3616 ENDDO 3617 ENDDO 3618C 3619 KOFF3 = KSCR1 + IMATIJ(ISYMI,ISYMM) 3620C 3621 ENDDO 3622C 3623C - TMAT(fnil) * t2tp(BlmD xovvo(fnma) 3624C 3625 DO ISYML = 1,NSYM 3626 ISYMFNI = MULD2H(ISCKIJ,ISYML) 3627 ISYMM = MULD2H(ISYMLM,ISYML) 3628C 3629 KOFF1 = ISAIKJ(ISYMFNI,ISYML) + 1 3630 KOFF2 = KSCR1 + IMATIJ(ISYML,ISYMM) 3631 KOFF3 = KSCR2 + ISAIKJ(ISYMFNI,ISYMM) 3632C 3633 NTOTL = MAX(1,NRHF(ISYML)) 3634 NTOTFNI = MAX(1,NCKI(ISYMFNI)) 3635C 3636 CALL DGEMM('N','N',NCKI(ISYMFNI),NRHF(ISYMM),NRHF(ISYML), 3637 * -ONE,TMAT(KOFF1),NTOTFNI,WORK(KOFF2),NTOTL, 3638 * ONE,WORK(KOFF3),NTOTFNI) 3639C 3640 END DO 3641C 3642C sort indices of result vedtor FNIM as FNMI 3643C 3644 IOPT = 2 3645 CALL CC3_LSORT2(WORK(KSCR2),ISYFNIM,WORK(KEND1),LWRK1,IOPT) 3646C 3647C omega(ai) = omega(ai) + xoovv(fnma) * tt2(fnmi) 3648C 3649 3650 DO ISYMI = 1,NSYM 3651 ISYMFNM = MULD2H(ISYFNIM,ISYMI) 3652 ISYMA = MULD2H(ISYRES,ISYMI) 3653C 3654 KOFF1 = IT2SP(ISYMFNM,ISYMA) + 1 3655 KOFF2 = KSCR2 + ISAIKJ(ISYMFNM,ISYMI) 3656 KOFF3 = IT1AM(ISYMA,ISYMI) + 1 3657C 3658 NTOTFNM = MAX(1,NCKI(ISYMFNM)) 3659 NTOTA = MAX(1,NVIR(ISYMA)) 3660C 3661C 5.3 omega contribution addomega 3662C 3663 CALL DGEMM('T','N',NVIR(ISYMA),NRHF(ISYMI),NCKI(ISYMFNM), 3664 * -ONE,XOVVO(KOFF1),NTOTFNM,WORK(KOFF2),NTOTFNM, 3665 * ONE,OMEGA1(KOFF3),NTOTA) 3666C 3667 END DO 3668C 3669 END IF 3670C 3671C================================================ 3672C 6. Calculate contributions from g_{oovv} 3673C - L^{def}_{lni} t^{de}_{lm} g_{mnfa} 3674C -( Wfe(dlni) + Wfd(enli) + Wde(finl) )* t(dl,em) * g(mnfa) 3675C================================================ 3676C 3677 ISYDLM = MULD2H(ISYMT2,ISYMID) 3678 ISYTMP = MULD2H(ISYDLM,ISCKIJ) 3679C 3680 KSCR1 = 1 3681 KSCR2 = KSCR1 + NMAIJK(ISYTMP) 3682 KEND1 = KSCR2 + NCKI(ISYDLM) 3683 LWRK1 = LWORK - KEND1 3684C 3685 IF (LWRK1 .LT. 0) THEN 3686 CALL QUIT('Out of memory CC3_W3_OMEGA1 (OVVO-4)') 3687 ENDIF 3688C 3689 CALL DZERO(WORK(KSCR1),NMAIJK(ISYTMP)) 3690C 3691C----------------------------------------------- 3692C Calculate first and second contribution 3693C -( Wfe(dlni) + Wfd(enli) )* t(dl,em) * g(mnfa) 3694C----------------------------------------------- 3695C 3696C - Wfe(dlni) * t(dl,em) * g(mnfa) 3697C 3698C Tt2(inm) = TMAT(dlin) * t2tp(dlmD) 3699C 3700C - Wfd(enli) * t(dl,em) * g(mnfa) 3701C 3702C Tt2(inm) = Tt2(inm) + TMAT(elin) * t2tp(emlD) 3703C tD(eml) 3704C Tt2(inm) = T(mni) 3705C 3706C g(mnBa) = xoovv(Bmna) = IB(mna) 3707C 3708C omega(ai) = omega(ai) + IB(mna) * T(mni) 3709C 3710C----------------------------------------------- 3711C First contribution to intermediate 3712C----------------------------------------------- 3713C 3714 DO I = 1, LENGTH 3715 TMAT(I) = WMAT(INDSQ(I,3)) 3716 ENDDO 3717C 3718 IF (NSYM .GT. 1) THEN 3719 IF (LWRK1 .LT. LENGTH) THEN 3720 CALL QUIT('Out of memory in CC3_W3_OMEGA1 (CC_GATHER-6)') 3721 ENDIF 3722 CALL CC_GATHER(LENGTH,WORK(KEND1),TMAT,INDSQ(1,6)) 3723 CALL DCOPY(LENGTH,WORK(KEND1),1,TMAT,1) 3724 ENDIF 3725C 3726C Tt2(inm) = TMAT(dlin) * t2tp(dlmD) 3727C 3728 DO ISYMDL = 1, NSYM 3729 ISYMM = MULD2H(ISYDLM,ISYMDL) 3730 ISYMIN = MULD2H(ISCKIJ,ISYMDL) 3731C 3732 KOFF1 = ISAIKL(ISYMDL,ISYMIN) 3733 * + 1 3734 KOFF2 = IT2SP(ISYDLM,ISYMID) 3735 * + NCKI(ISYDLM)*(ID-1) 3736 * + ICKI(ISYMDL,ISYMM) 3737 * + 1 3738 KOFF3 = KSCR1 3739 * + IMAIJK(ISYMIN,ISYMM) 3740C 3741 NTOTDL = MAX(1,NT1AM(ISYMDL)) 3742 NTOTIN = MAX(1,NMATIJ(ISYMIN)) 3743C 3744 CALL DGEMM('T','N',NMATIJ(ISYMIN),NRHF(ISYMM),NT1AM(ISYMDL), 3745 * -ONE,TMAT(KOFF1),NTOTDL,T2TP(KOFF2),NTOTDL, 3746 * ONE,WORK(KOFF3),NTOTIN) 3747 ENDDO 3748C 3749C----------------------------------------------- 3750C Second contribution to intermediate 3751C----------------------------------------------- 3752C 3753 ! Calculate only when contracting with first-order multipliers 3754 IF (W3X) THEN 3755C 3756 3757 DO I = 1, LENGTH 3758 TMAT(I) = WMAT(INDSQ(I,4)) 3759 ENDDO 3760C 3761C Tt2(inm) = Tt2(inm) + TMAT(elin) * t2tp(emlD) 3762C tD(eml) 3763 IF (NSYM .GT. 1) THEN 3764 IF (LWRK1 .LT. LENGTH) THEN 3765 CALL QUIT('Out of memory in CC3_W3_OMEGA1 (CC_GATHER-6)') 3766 ENDIF 3767 CALL CC_GATHER(LENGTH,WORK(KEND1),TMAT,INDSQ(1,6)) 3768 CALL DCOPY(LENGTH,WORK(KEND1),1,TMAT,1) 3769 ENDIF 3770C 3771 DO ISYMM = 1, NSYM 3772 ISYMDL = MULD2H(ISYDLM,ISYMM) 3773 DO ISYMD = 1, NSYM 3774 ISYML = MULD2H(ISYMDL,ISYMD) 3775 ISYMDM = MULD2H(ISYMD,ISYMM) 3776 DO M = 1, NRHF(ISYMM) 3777 DO L = 1, NRHF(ISYML) 3778C 3779 KOFF1 = IT2SP(ISYDLM,ISYMID) 3780 * + NCKI(ISYDLM)*(ID-1) 3781 * + ICKI(ISYMDL,ISYMM) 3782 * + NT1AM(ISYMDL)*(M-1) 3783 * + IT1AM(ISYMD,ISYML) 3784 * + NVIR(ISYMD)*(L-1) 3785 * + 1 3786C 3787 KOFF2 = KSCR2 3788 * + ICKI(ISYMDM,ISYML) 3789 * + NT1AM(ISYMDM)*(L-1) 3790 * + IT1AM(ISYMD,ISYMM) 3791 * + NVIR(ISYMD)*(M-1) 3792C 3793 CALL DCOPY(NVIR(ISYMD),T2TP(KOFF1),1, 3794 * WORK(KOFF2),1) 3795C 3796 ENDDO 3797 ENDDO 3798 ENDDO 3799 ENDDO 3800C 3801 DO ISYMDL = 1, NSYM 3802 ISYMM = MULD2H(ISYDLM,ISYMDL) 3803 ISYMIN = MULD2H(ISCKIJ,ISYMDL) 3804C 3805 KOFF1 = ISAIKL(ISYMDL,ISYMIN) 3806 * + 1 3807 KOFF2 = KSCR2 3808 * + ICKI(ISYMDL,ISYMM) 3809 KOFF3 = KSCR1 3810 * + IMAIJK(ISYMIN,ISYMM) 3811C 3812 NTOTDL = MAX(1,NT1AM(ISYMDL)) 3813 NTOTIN = MAX(1,NMATIJ(ISYMIN)) 3814C 3815 CALL DGEMM('T','N',NMATIJ(ISYMIN),NRHF(ISYMM),NT1AM(ISYMDL), 3816 * -ONE,TMAT(KOFF1),NTOTDL,WORK(KOFF2),NTOTDL, 3817 * ONE,WORK(KOFF3),NTOTIN) 3818 ENDDO 3819C 3820 END IF 3821C 3822C------------------------------------------------------- 3823C Sort intermediate and integrals and contract 3824C------------------------------------------------------- 3825C 3826 KEND1 = KSCR2 + NMAIJK(ISYTMP) 3827 LWRK1 = LWORK - KEND1 3828C 3829 IF (LWRK1 .LT. 0) THEN 3830 CALL QUIT('Out of memory CC3_W3_OMEGA1 (OVVO-5)') 3831 ENDIF 3832C 3833C Tt2(inm) = T(mni) 3834C 3835 DO ISYMM = 1, NSYM 3836 ISYMIN = MULD2H(ISYTMP,ISYMM) 3837 DO ISYMN = 1, NSYM 3838 ISYMI = MULD2H(ISYMIN,ISYMN) 3839 ISYMMN = MULD2H(ISYMM,ISYMN) 3840C 3841 DO M = 1, NRHF(ISYMM) 3842 DO N = 1, NRHF(ISYMN) 3843C 3844 KOFF1 = KSCR1 3845 * + IMAIJK(ISYMIN,ISYMM) 3846 * + NMATIJ(ISYMIN)*(M-1) 3847 * + IMATIJ(ISYMI,ISYMN) 3848 * + NRHF(ISYMI)*(N-1) 3849C 3850 KOFF2 = KSCR2 - 1 3851 * + IMAIJK(ISYMMN,ISYMI) 3852 * + IMATIJ(ISYMM,ISYMN) 3853 * + NRHF(ISYMM)*(N-1) 3854 * + M 3855C 3856 CALL DCOPY(NRHF(ISYMI),WORK(KOFF1),1, 3857 * WORK(KOFF2),NMATIJ(ISYMMN)) 3858C 3859 ENDDO 3860 ENDDO 3861 ENDDO 3862 ENDDO 3863C 3864 CALL DCOPY(NMAIJK(ISYTMP),WORK(KSCR2),1,WORK(KSCR1),1) 3865C 3866 ISYAMN = MULD2H(ISYINT,ISYMIB) 3867C 3868C g(mnBa) = xoovv(Bmna) = IB(mna) 3869C 3870 DO ISYMA = 1, NSYM 3871 ISYMMN = MULD2H(ISYMA,ISYAMN) 3872 ISYBMN = MULD2H(ISYINT,ISYMA) 3873 DO ISYMM = 1, NSYM 3874 ISYMN = MULD2H(ISYMMN,ISYMM) 3875 ISYMBN = MULD2H(ISYBMN,ISYMM) 3876C 3877 DO M = 1, NRHF(ISYMM) 3878 DO N = 1, NRHF(ISYMN) 3879C 3880 KOFF1 = IT2SP(ISYBMN,ISYMA) 3881 * + ICKI(ISYMBN,ISYMM) 3882 * + NT1AM(ISYMBN)*(M-1) 3883 * + IT1AM(ISYMIB,ISYMN) 3884 * + NVIR(ISYMIB)*(N-1) 3885 * + IB 3886C 3887 KOFF2 = KSCR2 - 1 3888 * + IMAIJA(ISYMMN,ISYMA) 3889 * + IMATIJ(ISYMN,ISYMM) 3890 * + NRHF(ISYMN)*(M-1) 3891 * + N 3892C 3893 CALL DCOPY(NVIR(ISYMA),XOOVV(KOFF1),NCKI(ISYBMN), 3894 * WORK(KOFF2),NMATIJ(ISYMMN)) 3895C 3896 ENDDO 3897 ENDDO 3898 ENDDO 3899 ENDDO 3900C 3901C omega(ai) = omega(ai) + IB(mna) * T(mni) 3902C 3903 DO ISYMA = 1, NSYM 3904 ISYMI = MULD2H(ISYRES,ISYMA) 3905 ISYMMN = MULD2H(ISYAMN,ISYMA) 3906C 3907 KOFF1 = KSCR2 3908 * + IMAIJA(ISYMMN,ISYMA) 3909 KOFF2 = KSCR1 3910 * + IMAIJK(ISYMMN,ISYMI) 3911 KOFF3 = IT1AM(ISYMA,ISYMI) 3912 * + 1 3913C 3914 NTOTMN = MAX(1,NMATIJ(ISYMMN)) 3915 NTOTA = MAX(1,NVIR(ISYMA)) 3916C 3917C 6.1 + 6.2 omega contribution addomega 3918C 3919 CALL DGEMM('T','N',NVIR(ISYMA),NRHF(ISYMI),NMATIJ(ISYMMN), 3920 * -ONE,WORK(KOFF1),NTOTMN,WORK(KOFF2),NTOTMN, 3921 * ONE,OMEGA1(KOFF3),NTOTA) 3922 ENDDO 3923C----------------------------------------------- 3924C Calculate third contribution 3925C - Wde(finl) * t(dl,em) * g(mnfa) 3926C---------------------------------------------- 3927C 3928C - TMAT(finl) * t2tp(BlmD) xoovv(fmna) 3929C 3930C Tt2(finm) * xoovv(fmna) 3931C tt2(fmni) 3932C 3933C omega(ai) = omega(ai) + xoovv(fmna) * tt2(fmni) 3934C 3935C 3936 ! Calculate only when contracting with first-order multipliers 3937 IF (W3X) THEN 3938C 3939 3940 ISYMBLM = MULD2H(ISYMT2,ISYMID) 3941 ISYMLM = MULD2H(ISYMBLM,ISYMIB) 3942 ISYFINM = MULD2H(ISYMLM,ISCKIJ) 3943C 3944 KSCR1 = 1 3945 KSCR2 = KSCR1 + NMATIJ(ISYMLM) 3946 KEND1 = KSCR2 + NCKIJ(ISYFINM) 3947 LWRK1 = LWORK - KEND1 3948C 3949 IF (LWRK1 .LT. 0) THEN 3950 CALL QUIT('Out of memory CC3_W3_OMEG 6.)') 3951 ENDIF 3952C 3953C 3954 CALL DZERO(WORK(KSCR2),NCKIJ(ISYFINM)) 3955C 3956 DO I = 1, LENGTH 3957 TMAT(I) = WMAT(I) 3958 ENDDO 3959C 3960C ------------------------------- 3961C Sort T^{BD}_{mi} as T_{mi} 3962C ------------------------------- 3963C 3964 ISYMMI = MULD2H(ISYMBD,ISYMT2) 3965 ISYBMI = MULD2H(ISYMMI,ISYMIB) 3966C 3967 DO ISYMI = 1,NSYM 3968 ISYMM = MULD2H(ISYMMI,ISYMI) 3969 ISYMBM = MULD2H(ISYBMI,ISYMI) 3970 DO I = 1,NRHF(ISYMI) 3971 DO M = 1,NRHF(ISYMM) 3972 KOFF1 = IT2SP(ISYBMI,ISYMID) 3973 * + NCKI(ISYBMI)*(ID-1) 3974 * + ICKI(ISYMBM,ISYMI) 3975 * + NT1AM(ISYMBM)*(I-1) 3976 * + IT1AM(ISYMIB,ISYMM) 3977 * + NVIR(ISYMIB)*(M-1) 3978 * + IB 3979 KOFF2 = IMATIJ(ISYMM,ISYMI) 3980 * + NRHF(ISYMM)*(I-1) 3981 * + M 3982 WORK(KSCR1-1+KOFF2) = T2TP(KOFF1) 3983 ENDDO 3984 ENDDO 3985 ENDDO 3986C 3987C - TMAT(finl) * t2tp(BlmD) xoovv(fmna) 3988C 3989 ISYMLM = MULD2H(ISYMBD,ISYMT2) 3990 DO ISYML = 1,NSYM 3991 ISYMFIN = MULD2H(ISCKIJ,ISYML) 3992 ISYMM = MULD2H(ISYMLM,ISYML) 3993C 3994 KOFF1 = ISAIKJ(ISYMFIN,ISYML) + 1 3995 KOFF2 = KSCR1 + IMATIJ(ISYML,ISYMM) 3996 KOFF3 = KSCR2 + ISAIKJ(ISYMFIN,ISYMM) 3997C 3998 NTOTL = MAX(1,NRHF(ISYML)) 3999 NTOTFIN = MAX(1,NCKI(ISYMFIN)) 4000C 4001 CALL DGEMM('N','N',NCKI(ISYMFIN),NRHF(ISYMM),NRHF(ISYML), 4002 * -ONE,TMAT(KOFF1),NTOTFIN,WORK(KOFF2),NTOTL, 4003 * ONE,WORK(KOFF3),NTOTFIN) 4004 END DO 4005C 4006C sort indices of result vedtor FINM as FMNI 4007C 4008 IOPT = 4 4009 CALL CC3_LSORT2(WORK(KSCR2),ISYFINM,WORK(KEND1),LWRK1,IOPT) 4010C 4011C omega(ai) = omega(ai) + xoovv(fmna) * tt2(fmni) 4012C 4013 DO ISYMI = 1,NSYM 4014 ISYMFMN = MULD2H(ISYFINM,ISYMI) 4015 ISYMA = MULD2H(ISYRES,ISYMI) 4016C 4017 KOFF1 = IT2SP(ISYMFMN,ISYMA) + 1 4018 KOFF2 = KSCR2 + ISAIKJ(ISYMFMN,ISYMI) 4019 KOFF3 = IT1AM(ISYMA,ISYMI) + 1 4020C 4021 NTOTFMN = MAX(1,NCKI(ISYMFMN)) 4022 NTOTA = MAX(1,NVIR(ISYMA)) 4023C 4024C 6.3 omega contribution addomega 4025C 4026 CALL DGEMM('T','N',NVIR(ISYMA),NRHF(ISYMI),NCKI(ISYMFMN), 4027 * -ONE,XOOVV(KOFF1),NTOTFMN,WORK(KOFF2),NTOTFMN, 4028 * ONE,OMEGA1(KOFF3),NTOTA) 4029 END DO 4030C 4031 END IF 4032C 4033C------------------- 4034C End. 4035C------------------- 4036C 4037 CALL QEXIT('CC3_W3_OMEGA1') 4038C 4039 RETURN 4040 END 4041C /* Deck cc3_srtvooo */ 4042 SUBROUTINE CC3_SRTVOOO(GVPOOO,GVOOO,ISCKIJ) 4043C 4044C Sort matrix of symmetry ISCKIJ stored in GVOOO(V,O,O,O) 4045C as GVPOOO(V,OOO) 4046C 4047C 4048 IMPLICIT NONE 4049C 4050 4051 INTEGER ISCKIJ, ISYMJ, ISYCKI, ISYMI, ISYMCK, ISYMIJ, ISYMC, ISYMK 4052 INTEGER ISYMKI, ISYKIJ, NKIJ, KVPOOO, KVOOO 4053C 4054#if defined (SYS_CRAY) 4055 REAL GVPOOO(*),GVOOO(*) 4056#else 4057 DOUBLE PRECISION GVPOOO(*),GVOOO(*) 4058#endif 4059C 4060#include "priunit.h" 4061#include "ccorb.h" 4062#include "ccsdsym.h" 4063 4064C 4065 CALL QENTER('CC3_SRTVOOO') 4066C 4067 DO 100 ISYMJ = 1,NSYM 4068C 4069 ISYCKI = MULD2H(ISYMJ,ISCKIJ) 4070C 4071 DO 110 ISYMI = 1,NSYM 4072C 4073 ISYMCK = MULD2H(ISYMI,ISYCKI) 4074 ISYMIJ = MULD2H(ISYMI,ISYMJ) 4075C 4076 DO 120 ISYMC = 1,NSYM 4077C 4078 ISYMK = MULD2H(ISYMCK,ISYMC) 4079 ISYMKI = MULD2H(ISYMK,ISYMI) 4080 ISYKIJ = MULD2H(ISYMIJ,ISYMK) 4081C 4082 DO 130 J = 1,NRHF(ISYMJ) 4083C 4084 DO 140 I = 1,NRHF(ISYMI) 4085C 4086 DO 150 K = 1,NRHF(ISYMK) 4087C 4088 NKIJ = IMAIJK(ISYMKI,ISYMJ) 4089 * + NMATIJ(ISYMKI)*(J-1) 4090 * + IMATIJ(ISYMK,ISYMI) 4091 * + NRHF(ISYMK)*(I - 1) + K 4092C 4093 DO 160 C = 1,NVIR(ISYMC) 4094C 4095 KVOOO = ISAIKJ(ISYCKI,ISYMJ) 4096 * + NCKI(ISYCKI)*(J-1) 4097 * + ISAIK(ISYMCK,ISYMI) 4098 * + NT1AM(ISYMCK)*(I - 1) 4099 * + IT1AM(ISYMC,ISYMK) 4100 * + NVIR(ISYMC)*(K - 1) 4101 * + C 4102C 4103 KVPOOO = I3VOOO(ISYMC,ISYKIJ) 4104 * + NVIR(ISYMC)*(NKIJ - 1) 4105 * + C 4106C 4107 GVPOOO(KVPOOO) = GVOOO(KVOOO) 4108C 4109 160 CONTINUE 4110 150 CONTINUE 4111 140 CONTINUE 4112 130 CONTINUE 4113 120 CONTINUE 4114 110 CONTINUE 4115 100 CONTINUE 4116C 4117 CALL QEXIT('CC3_SRTVOOO') 4118C 4119 RETURN 4120 END 4121 4122C /* Deck cc3_w3_cy2v */ 4123 SUBROUTINE CC3_W3_CY2V(XI2,ISYRES,RBJIA,WMAT,ISWMAT,TMAT,TRVIR, 4124 * TRVIR1,ISYINT,WORK,LWORK,INDSQ,LENSQ, 4125 * ISYMIB,IB,ISYMID,ID,W3X) 4126C 4127C 4128C Modified based on CC3_CY2V routine where TWO*WMAT is replaced by WMAT 4129C and the other terms in TMAT are removed. 4130C 4131C General symmetry: ISWMAT is symmetry of WMAT and TMAT intermediates. 4132C (exclusive isymib,isymid) 4133C ISYINT is symmetry of integrals in TRVIR and TROVIR1. 4134C ISYRES = ISWMAT*ISYINT*ISYMIB*ISYMID 4135C RBJIA intermediate result vector 4136C 4137C 4138C If W3X = .TRUE., then WMAT contains wbar_X and all terms 4139C 1-3 are included 4140C else 4141C we assume that we get tbar_0 in as WMAT and calculate only 4142C term 3. 4143C 4144C Note Wbar_X and Tbar_0 is stored in the same way for fixed IB and ID 4145C 4146C XI2(aibj) = XI2(aibj) 4147C 4148C + sum_{cdk} (2W^{bcd}_{jik}-W^{bcd}_{kij}-W^{bcd}_{jki}) * (ac|kd) 4149C (1): 4150C = sum_{ck} (2W^{DB}_{cijk}-W^{DB}_{cikj}-W^{DB}_{ckji}) * (ac|kD) 4151C (2): 4152C + sum_{dk} (2W^{BC}_{dkij}-W^{BC}_{djik}-W^{BC}_{dikj}) * (aC|kd) 4153C (3): 4154C + sum_{k} (2W^{CD}_{bjki}-W^{CD}_{bkji}-W^{CD}_{bjik}) * (aC|kD) 4155C 4156 IMPLICIT NONE 4157C 4158 LOGICAL W3X 4159C 4160 INTEGER ISYRES, ISWMAT, ISYINT , LENSQ, LWORK 4161 INTEGER INDSQ(LENSQ,6) 4162 INTEGER ISYMIB, IB, ISYMID, ID 4163 INTEGER INDEX ,LENGTH, ISCKIJ, ISYMB, ISYAIJ, KRMAT, KEND, LEND 4164 INTEGER ISYMJ ,ISYBJ, ISYAI, ISYCKI, ISYCK, ISYMA, NTOTCK, NVIRA 4165 INTEGER KOFF1 , KOFF2, KOFF3, ISYMD, ISDKIJ, ISYDKI 4166 INTEGER ISYDK, NTOTDK, ISBJIK, ISYCKA 4167 INTEGER ISYKA, KINT, ISYBJI, KOFFT, KOFFM, KOFFR 4168 INTEGER NTOK, NTOBJI, ISYMK, ISYMC, ISYMI 4169#if defined (SYS_CRAY) 4170 REAL XI2(*), RBJIA(*), WMAT(*), TMAT(*), TRVIR(*), TRVIR1(*) 4171 REAL WORK(*) 4172 REAL ZERO, ONE, TWO 4173#else 4174 DOUBLE PRECISION XI2(*), RBJIA(*), WMAT(*), TMAT(*) 4175 DOUBLE PRECISION TRVIR(*), TRVIR1(*), WORK(*) 4176 DOUBLE PRECISION ZERO, ONE, TWO 4177#endif 4178 4179C 4180 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 4181C 4182#include "priunit.h" 4183#include "ccsdsym.h" 4184#include "ccorb.h" 4185#include "ccsdinp.h" 4186C 4187C 4188C INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 4189C 4190 CALL QENTER('CC3_W3_CY2V') 4191 LENGTH = NCKIJ(ISWMAT) 4192C 4193C------------------------ 4194C (1): 4195C = sum_{ck} (2W^{DB}_{cijk}-W^{DB}_{cikj}-W^{DB}_{ckji}) * (ac|kD) 4196C 4197C------------------------ 4198C 4199C TMAT(ckij) = 2W^{DB}_{cijk}-W^{DB}_{cikj}-W^{DB}_{ckji} 4200C 4201C Use here the symmetry between BJ, DK !! 4202C 4203 ! Calculate only when contracting with first-order multipliers 4204 IF (W3X) THEN 4205C 4206 4207 ISCKIJ = ISWMAT 4208C 4209 ISYMB = ISYMIB 4210 ISYMD = ISYMID 4211 B = IB 4212 D = ID 4213C 4214 ISYAIJ = MULD2H(ISYMB,ISYRES) 4215 KRMAT = 1 4216 KEND = KRMAT + NCKI(ISYAIJ) 4217 LEND = LWORK - KEND 4218 IF (LEND .LT. 0)THEN 4219 CALL QUIT('1. Insufficient core in CC3_W3_CY2V') 4220 ENDIF 4221 CALL DZERO(WORK(KRMAT), NCKI(ISYAIJ)) 4222 DO I = 1,LENGTH 4223C 4224 TMAT(I) = WMAT(INDSQ(I,1)) 4225C 4226 ENDDO 4227C 4228C--------------------------- 4229C (ac|kD) = TRVIR^D_{ck,a} 4230C 4231C RMAT^B_(aij) = TRVIR^D_{ck,a} )^T * TMAT(ckij) 4232C--------------------------- 4233C 4234 DO 200 ISYMJ = 1,NSYM 4235C 4236 ISYBJ = MULD2H(ISYMB,ISYMJ) 4237 ISYAI = MULD2H(ISYBJ,ISYRES) 4238 ISYCKI = MULD2H(ISCKIJ,ISYMJ) 4239C 4240 DO 210 J = 1,NRHF(ISYMJ) 4241C 4242 DO 220 ISYMI = 1,NSYM 4243C 4244 ISYCK = MULD2H(ISYCKI,ISYMI) 4245 ISYMA = MULD2H(ISYAI,ISYMI) 4246C 4247 NTOTCK = MAX(NT1AM(ISYCK),1) 4248 NVIRA = MAX(NVIR(ISYMA),1) 4249C 4250 KOFF1 = ICKATR(ISYCK,ISYMA) + 1 4251 KOFF2 = ISAIKJ(ISYCKI,ISYMJ) 4252 * + NCKI(ISYCKI)*(J - 1) 4253 * + ISAIK(ISYCK,ISYMI) + 1 4254 KOFF3 = KRMAT 4255 * + ISAIK(ISYAI,ISYMJ) 4256 * + NT1AM(ISYAI)*(J - 1) 4257 * + IT1AM(ISYMA,ISYMI) 4258C 4259 CALL DGEMM('T','N',NVIR(ISYMA),NRHF(ISYMI), 4260 * NT1AM(ISYCK), 4261 * -ONE,TRVIR(KOFF1),NTOTCK,TMAT(KOFF2),NTOTCK, 4262 * ONE,WORK(KOFF3),NVIRA) 4263C 4264 220 CONTINUE 4265C 4266 210 CONTINUE 4267C 4268 200 CONTINUE 4269C 4270 ISYAIJ = MULD2H(ISYMB,ISYRES) 4271 IF ( NCKI(ISYAIJ).GT.0 ) THEN 4272 CALL CC3_RACC(XI2,WORK(KRMAT),ISYMB,B,ISYRES) 4273 END IF 4274C 4275C 4276C (2): 4277C 4278C XI2(aibj) = XI2(aibj) 4279C 4280C + sum_{dk} (2W^{BC}_{dkij}-W^{BC}_{djik}-W^{BC}_{dikj}) * (aC|kd) 4281C 4282C TMAT(dkij) = 2W^{BC}_{dkij}-W^{BC}_{djik}-W^{BC}_{dikj} 4283C 4284 ISYMC = ISYMID 4285 ISYMB = ISYMIB 4286 B = IB 4287 C = ID 4288C 4289 ISDKIJ = ISWMAT 4290 ISYAIJ = MULD2H(ISYMB,ISYRES) 4291 KRMAT = 1 4292 KEND = KRMAT + NCKI(ISYAIJ) 4293 LEND = LWORK - KEND 4294 IF (LEND .LT. 0)THEN 4295 CALL QUIT('2. Insufficient core in CC3_W3_CY2V') 4296 ENDIF 4297 CALL DZERO(WORK(KRMAT), NCKI(ISYAIJ)) 4298C 4299 DO I = 1,LENGTH 4300C 4301 TMAT(I) = WMAT(I) 4302C 4303 ENDDO 4304C 4305C (ad|kC) = TRVIR1^C(dka) 4306C 4307C RMAT^B_(aij) = ( TRVIR1^C(dka) )^T * TMAT(dkij) 4308C 4309 DO 300 ISYMJ = 1,NSYM 4310C 4311 ISYBJ = MULD2H(ISYMB,ISYMJ) 4312 ISYAI = MULD2H(ISYBJ,ISYRES) 4313 ISYDKI = MULD2H(ISDKIJ,ISYMJ) 4314C 4315 DO 310 J = 1,NRHF(ISYMJ) 4316C 4317 DO 320 ISYMI = 1,NSYM 4318C 4319 ISYDK = MULD2H(ISYDKI,ISYMI) 4320 ISYMA = MULD2H(ISYAI,ISYMI) 4321C 4322 NTOTDK = MAX(NT1AM(ISYDK),1) 4323 NVIRA = MAX(NVIR(ISYMA),1) 4324C 4325 KOFF1 = ICKATR(ISYDK,ISYMA) + 1 4326 KOFF2 = ISAIKJ(ISYDKI,ISYMJ) 4327 * + NCKI(ISYDKI)*(J - 1) 4328 * + ISAIK(ISYDK,ISYMI) + 1 4329 KOFF3 = KRMAT + ISAIK(ISYAI,ISYMJ) 4330 * + NT1AM(ISYAI)*(J - 1) 4331 * + IT1AM(ISYMA,ISYMI) 4332C 4333 CALL DGEMM('T','N',NVIR(ISYMA),NRHF(ISYMI), 4334 * NT1AM(ISYDK), 4335 * -ONE,TRVIR1(KOFF1),NTOTDK,TMAT(KOFF2),NTOTDK, 4336 * ONE,WORK(KOFF3),NVIRA) 4337C 4338 320 CONTINUE 4339 310 CONTINUE 4340 300 CONTINUE 4341C 4342 ISYAIJ = MULD2H(ISYMB,ISYRES) 4343 IF ( NCKI(ISYAIJ).GT.0 ) THEN 4344 CALL CC3_RACC(XI2,WORK(KRMAT),ISYMB,B,ISYRES) 4345 END IF 4346C 4347 END IF 4348C 4349C----------------------------------- 4350C (3): 4351C + sum_{k} (2W^{CD}_{bjki}-W^{CD}_{bkji}-W^{CD}_{bjik}) * (ac|kD) 4352C 4353C----------------------------------- 4354C 4355 ISYMC = ISYMIB 4356 ISYMD = ISYMID 4357 C = IB 4358 D = ID 4359C 4360C TMAT(bjik) = 2W^{CD}_{bjki}-W^{CD}_{bkji}-W^{CD}_{bjik} 4361C 4362 ISBJIK = ISWMAT 4363C 4364 DO I = 1,LENGTH 4365C 4366 TMAT(I) = WMAT(INDSQ(I,3)) 4367C 4368 ENDDO 4369C 4370C (ac|kD) = TRVIR^D_{Cka} = M^DC_{ka} 4371C 4372c 4373 ISYCKA = MULD2H(ISYINT,ISYMD) 4374 ISYKA = MULD2H(ISYCKA,ISYMC) 4375C 4376 DO ISYMA = 1,NSYM 4377 ISYCK = MULD2H(ISYCKA,ISYMA) 4378 ISYMK = MULD2H(ISYMC,ISYCK) 4379 DO A = 1,NVIR(ISYMA) 4380 DO K = 1,NRHF(ISYMK) 4381 KOFFM = IT1AMT(ISYMK,ISYMA) 4382 * + NRHF(ISYMK)*(A-1) + K 4383 KINT = ICKATR(ISYCK,ISYMA) 4384 * + NT1AM(ISYCK)*(A-1) 4385 * + IT1AM(ISYMC,ISYMK) 4386 * + NVIR(ISYMC)*(K-1) + C 4387 WORK(KOFFM) = TRVIR(KINT) 4388 END DO 4389 END DO 4390 END DO 4391C 4392C RBJIA(bjia) = RBJIA(bjia) + TMAT(bjik) * M^DC_{ka} 4393C 4394 DO ISYMA = 1,NSYM 4395 ISYMK = MULD2H(ISYMA,ISYKA) 4396 ISYBJI = MULD2H(ISBJIK,ISYMK) 4397 KOFFT = ISAIKJ(ISYBJI,ISYMK) + 1 4398 KOFFM = IT1AMT(ISYMK,ISYMA) + 1 4399 KOFFR = IT2SP(ISYBJI,ISYMA) + 1 4400 NTOK = MAX(NRHF(ISYMK),1) 4401 NTOBJI = MAX(NCKI(ISYBJI),1) 4402C 4403 CALL DGEMM('N','N',NCKI(ISYBJI),NVIR(ISYMA), 4404 * NRHF(ISYMK),-ONE,TMAT(KOFFT),NTOBJI, 4405 * WORK(KOFFM),NTOK,ONE,RBJIA(KOFFR),NTOBJI) 4406C 4407 END DO 4408C 4409 CALL QEXIT('CC3_W3_CY2V') 4410 RETURN 4411 END 4412 4413C /* Deck cc3_w3_cy2o */ 4414 SUBROUTINE CC3_W3_CY2O(XI2,ISYRES,WMAT,ISWMAT,TMAT,TROCC, 4415 * TROCC1,ISYINT,WORK,LWORK,INDSQ,LENSQ, 4416 * ISYMIB,IB,ISYMID,ID,W3X) 4417C 4418C 4419C 4420C Modified based on CC3_CY2O routine where TWO*WMAT is replaced by WMAT 4421C and the other terms in TMAT are removed. 4422C 4423C General symmetry: ISWMAT is symmetry of WMAT and TMAT intermediates. 4424C (exclusive isymib,isymid) 4425C ISYINT is symmetry of integrals in TROCC and TROCC1. 4426C ISYRES = ISWMAT*ISYINT*ISYMIB*ISYMID 4427C 4428C If W3X = .TRUE., then WMAT contains wbar_X and all terms 4429C 1-3 are included 4430C else 4431C we assume that we get tbar_0 in as WMAT and calculate only 4432C term 1. 4433C 4434C Note Wbar_X and Tbar_0 is stored in the same way for fixed IB and ID 4435 4436C 4437C XI2(aibj) = XI2(aibj) 4438C 4439C + sum_{ckl} (2W^{abc}_{ikl}-W^{abc}_{lki}-W^{abc}_{ilk}) * (lc|kj) 4440C (1): 4441C = sum_{kl} (2W^{BC}_{ailk}-W^{BC}_{alik}-W^{BC}_{aikl}) * (lc|kj) 4442C (2): 4443C + sum_{kl} (2W^{CA}_{bkil}-W^{CA}_{bkli}-W^{CA}_{blik}) * (lc|kj) 4444C (3): 4445C + sum_{ckl} (2W^{AB}_{clki}-W^{AB}_{cikl}-W^{AB}_{ckli}) * (lc|kj) 4446C 4447 IMPLICIT NONE 4448C 4449 LOGICAL W3X 4450C 4451 INTEGER ISYRES, ISWMAT, ISYINT, LWORK, LENSQ 4452 INTEGER ISYMIB, IB, ISYMID, ID 4453 INTEGER INDSQ(LENSQ,6) 4454 INTEGER ISYMC, ISYMB, ISYBC, JSAIKL, LENGTH, ISYMJ, ISYBJ 4455 INTEGER ISYAI, ISYKL, NTOTAI 4456 INTEGER NTOTKL, KOFF1, KOFF2, KOFF3 4457 INTEGER ISYMA, ISYAB, JSCKLI, JSBIKL, ISYMI 4458 INTEGER ISYCKL, NTOCKL, NRHFI 4459 INTEGER KRMAT1, KRMAT2, KEND, LEND, INDEX 4460 INTEGER ISYCK, ISYCJ, ISYAC, ISYBI, ISYKLJ, NTOTBI 4461 INTEGER ISWB, ISWD, NCKIMX 4462 INTEGER ISYAIJ, ISYBJI, ISYIJ 4463 INTEGER KTMATX 4464C 4465#if defined (SYS_CRAY) 4466 REAL XI2(*), WMAT(*), TMAT(*), TROCC(*), TROCC1(*) 4467 REAL WORK(*) 4468 REAL TWO, ONE, ZERO, XTMAT, XRMAT, DDOT 4469#else 4470 DOUBLE PRECISION XI2(*), WMAT(*), TMAT(*) 4471 DOUBLE PRECISION TROCC(*), TROCC1(*) 4472 DOUBLE PRECISION WORK(*) 4473 DOUBLE PRECISION TWO, ONE, ZERO, XTMAT, XRMAT, DDOT 4474#endif 4475C 4476#include "priunit.h" 4477#include "ccsdsym.h" 4478#include "ccorb.h" 4479#include "ccsdinp.h" 4480C 4481 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 4482C 4483C INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 4484C 4485 CALL QENTER('CC3_W3_CY2O') 4486C 4487 4488 ISWB = MULD2H(ISWMAT,ISYMIB) 4489 ISWD = MULD2H(ISWMAT,ISYMID) 4490 ISYAIJ = MULD2H(ISYMIB,ISYRES) 4491 ISYBJI = MULD2H(ISYMID,ISYRES) 4492 NCKIMX = MAX(NCKI(ISWB),NCKI(ISWD),NCKI(ISYAIJ),NCKI(ISYBJI)) 4493C 4494 KRMAT1 = 1 4495 KRMAT2 = KRMAT1 + NCKIMX 4496 KTMATX = KRMAT2 + NCKIMX 4497 KEND = KTMATX + NCKIJ(ISWMAT) 4498 LEND = LWORK - KEND 4499 IF (LWORK .LT. KEND) THEN 4500 CALL QUIT('Insufficient space in CY2O') 4501 ENDIF 4502C 4503C------------------------- 4504C First occupied term. 4505C------------------------- 4506C (1): 4507C XI2(aibj) = XI2(aibj) 4508C 4509C + sum_{kl} (2W^{BC}_{ailk}-W^{BC}_{alik}-W^{BC}_{aikl}) * (lc|kj) 4510 4511C 4512 B = IB 4513 C = ID 4514C 4515 ISYMB = ISYMIB 4516 ISYMC = ISYMID 4517C 4518 ISYAIJ = MULD2H(ISYMB,ISYRES) 4519 IF (NCKI(ISYAIJ).GT.0 ) THEN 4520 CALL DZERO(WORK(KRMAT1),NCKI(ISYAIJ)) 4521C 4522 ISYBC = MULD2H(ISYMB,ISYMC) 4523 JSAIKL = ISWMAT 4524C 4525 LENGTH = NCKIJ(JSAIKL) 4526C 4527C--------------------------------------- 4528C 4529C TMAT(aikl) = (2W^{BC}_{ailk}-W^{BC}_{alik}-W^{BC}_{aikl}) 4530C 4531C--------------------------------------- 4532C 4533 DO I = 1,LENGTH 4534C 4535 TMAT(I) = WMAT(INDSQ(I,3)) 4536C 4537 ENDDO 4538C 4539C---------------------------------- 4540C Symmetry sorting if symmetry. 4541C---------------------------------- 4542C 4543 4544 IF (NSYM .GT. 1) THEN 4545 CALL CC_GATHER(LENGTH,WORK(KTMATX),TMAT,INDSQ(1,6)) 4546 CALL DCOPY(LENGTH,WORK(KTMATX),1,TMAT,1) 4547 ENDIF 4548C 4549 IF (IPRINT .GT. 55) THEN 4550 XTMAT = DDOT(NCKIJ(JSAIKL),TMAT,1,TMAT,1) 4551 WRITE(LUPRI,*) 'In CC3_W3_CY2O: 1. Norm of TMAT = ',XTMAT 4552 ENDIF 4553C 4554C----------------------- 4555C First contraction. 4556C----------------------- 4557C 4558C TROCC(kljC) = (lc|kj) 4559C 4560C XI2(aibj) = XI2(aibj) 4561C + sum(kl) TMAT(aikl)* TROCC(kljC) 4562C 4563 DO 200 ISYMJ = 1,NSYM 4564C 4565 ISYBJ = MULD2H(ISYMB,ISYMJ) 4566 ISYAI = MULD2H(ISYBJ,ISYRES) 4567 ISYKL = MULD2H(JSAIKL,ISYAI) 4568 ISYKLJ = MULD2H(ISYKL,ISYMJ) 4569C 4570 NTOTAI = MAX(NT1AM(ISYAI),1) 4571 NTOTKL = MAX(NMATIJ(ISYKL),1) 4572C 4573 KOFF1 = ISAIKL(ISYAI,ISYKL) + 1 4574 KOFF2 = ISJIKA(ISYKLJ,ISYMC) 4575 * + NMAJIK(ISYKLJ)*(C - 1) 4576 * + ISJIK(ISYKL,ISYMJ) + 1 4577 KOFF3 = KRMAT1 + ISAIK(ISYAI,ISYMJ) 4578C 4579 CALL DGEMM('N','N',NT1AM(ISYAI),NRHF(ISYMJ),NMATIJ(ISYKL), 4580 * ONE,TMAT(KOFF1),NTOTAI,TROCC(KOFF2),NTOTKL, 4581 * ONE,WORK(KOFF3),NTOTAI) 4582C 4583 200 CONTINUE 4584C 4585 IF (IPRINT .GT. 55) THEN 4586 XRMAT = DDOT(NCKI(ISYAIJ),WORK(KRMAT1),1,WORK(KRMAT1),1) 4587 WRITE(LUPRI,*) '1. In CC3_W3_CY2O: Norm of RMAT1 = ',XRMAT 4588 ENDIF 4589C 4590 CALL CC3_RACC(XI2,WORK(KRMAT1),ISYMB,B,ISYRES) 4591C 4592 END IF 4593C 4594C-------------------------- 4595C Second occupied term. 4596C-------------------------- 4597C 4598C (2): 4599C 4600C XI2(aibj) = XI2(aibj) 4601C 4602C + sum_{kl} (2W^{CA}_{bkil}-W^{CA}_{bkli}-W^{CA}_{blik}) * (lc|kj) 4603C 4604 ! Calculate only when contracting with first-order multipliers 4605 IF (W3X) THEN 4606C 4607 4608 A = ID 4609 C = IB 4610C 4611 ISYMA = ISYMID 4612 ISYMC = ISYMIB 4613C 4614 ISYBJI = MULD2H(ISYMA,ISYRES) 4615 IF (NCKI(ISYBJI).GT.0) THEN 4616 CALL DZERO(WORK(KRMAT1),NCKI(ISYBJI)) 4617C 4618 ISYAC = MULD2H(ISYMA,ISYMC) 4619 JSBIKL = ISWMAT 4620C 4621C--------------------------------------- 4622C 4623C TMAT(bikl) = (2W^{CA}_{bkil}-W^{CA}_{bkli}-W^{CA}_{blik}) 4624C 4625C--------------------------------------- 4626C 4627 LENGTH = NCKIJ(JSBIKL) 4628C 4629 DO I = 1,LENGTH 4630C 4631 TMAT(I) = WMAT(INDSQ(I,1)) 4632C 4633 ENDDO 4634C 4635C---------------------------------- 4636C Symmetry sorting if symmetry. 4637C---------------------------------- 4638C 4639 IF (NSYM .GT. 1) THEN 4640 CALL CC_GATHER(LENGTH,WORK(KTMATX),TMAT,INDSQ(1,6)) 4641 CALL DCOPY(LENGTH,WORK(KTMATX),1,TMAT,1) 4642 ENDIF 4643C 4644 IF (IPRINT .GT. 55) THEN 4645 XTMAT = DDOT(NCKIJ(JSAIKL),TMAT,1,TMAT,1) 4646 WRITE(LUPRI,*) 'In CC3_W3_CY2O: 2. Norm of TMAT = ',XTMAT 4647 ENDIF 4648C 4649C------------------------ 4650C Second contraction. 4651C------------------------ 4652C 4653C TROCC(kljC) = (lc|kj) 4654C 4655C M^A_(bij) = sum(kl) TMAT(bikl)* TTROCC(kljC) 4656C 4657 4658 DO 400 ISYMJ = 1,NSYM 4659C 4660 ISYCJ = MULD2H(ISYMC,ISYMJ) 4661 ISYKL = MULD2H(ISYINT,ISYCJ) 4662 ISYBI = MULD2H(ISYKL,ISWMAT) 4663 ISYKLJ = MULD2H(ISYKL,ISYMJ) 4664C 4665 NTOTBI = MAX(NT1AM(ISYBI),1) 4666 NTOTKL = MAX(NMATIJ(ISYKL),1) 4667C 4668 KOFF1 = ISAIKL(ISYBI,ISYKL) + 1 4669 KOFF2 = ISJIKA(ISYKLJ,ISYMC) 4670 * + NMAJIK(ISYKLJ)*(C - 1) 4671 * + ISJIK(ISYKL,ISYMJ) + 1 4672 KOFF3 = KRMAT1 + ISAIK(ISYBI,ISYMJ) 4673C 4674 4675 CALL DGEMM('N','N',NT1AM(ISYBI),NRHF(ISYMJ),NMATIJ(ISYKL), 4676 * ONE,TMAT(KOFF1),NTOTBI,TROCC(KOFF2),NTOTKL, 4677 * ONE,WORK(KOFF3),NTOTBI) 4678C 4679 400 CONTINUE 4680C 4681 4682 IF (IPRINT .GT. 55) THEN 4683 XRMAT = DDOT(NCKI(ISYBJI),WORK(KRMAT1),1,WORK(KRMAT1),1) 4684 WRITE(LUPRI,*) '2.In CC3_W3_CY2O: Norm of RMAT1 = ',XRMAT 4685 ENDIF 4686C 4687C reorder result vector M^A_(bij) as M^A_(bji) 4688C 4689 CALL CC3_MAJI(WORK(KRMAT1),WORK(KRMAT2),ISYMA,A,ISYRES) 4690C 4691 IF (IPRINT .GT. 55) THEN 4692 XRMAT = DDOT(NCKI(ISYBJI),WORK(KRMAT2),1,WORK(KRMAT2),1) 4693 WRITE(LUPRI,*) 4694 * 'In CC3_W3_CY2O: Reorder :Norm of RMAT2 = ',XRMAT 4695 ENDIF 4696C 4697C add vector to XI2 4698C 4699 CALL CC3_RACC(XI2,WORK(KRMAT2),ISYMA,A,ISYRES) 4700C 4701 END IF 4702C 4703C------------------------- 4704C Third occupied term. 4705C------------------------- 4706C (3): 4707C XI2(aibj) = XI2(aibj) 4708C 4709C + sum_{ckl} (2W^{AB}_{clki}-W^{AB}_{cikl}-W^{AB}_{ckli}) * (lc|kj) 4710C 4711 B = ID 4712 A = IB 4713C 4714 ISYMB = ISYMID 4715 ISYMA = ISYMIB 4716C 4717 ISYAB = MULD2H(ISYMA,ISYMB) 4718 ISYIJ = MULD2H(ISYAB,ISYRES) 4719 IF (NMATIJ(ISYIJ).GT.0) THEN 4720 CALL DZERO(WORK,NMATIJ(ISYIJ)) 4721C 4722 JSCKLI = ISWMAT 4723C 4724 LENGTH = NCKIJ(JSCKLI) 4725C 4726C---------------------------------- 4727C 4728C TMAT(ckli) = 2W^{AB}_{clki}-W^{AB}_{cikl}-W^{AB}_{ckli} 4729C 4730C---------------------------------- 4731C 4732 DO I = 1,LENGTH 4733C 4734 TMAT(I) = WMAT(INDSQ(I,1)) 4735C 4736 ENDDO 4737C 4738C 4739 IF (IPRINT .GT. 55) THEN 4740 XTMAT = DDOT(NCKIJ(JSAIKL),TMAT,1,TMAT,1) 4741 WRITE(LUPRI,*) 'In CC3_W3_CY2O: 3. Norm of TMAT = ',XTMAT 4742 ENDIF 4743C 4744C----------------------- 4745C Third contraction. 4746C----------------------- 4747C 4748C 4749C TROCC1(cklj) = (lc|kj) 4750C 4751C XI2(aibj) = XI2(aibj) 4752C + sum(kl) TMAT(ckli)* TROCC1(cklj) 4753C 4754C 4755 DO 600 ISYMJ = 1,NSYM 4756C 4757 ISYBJ = MULD2H(ISYMB,ISYMJ) 4758 ISYAI = MULD2H(ISYBJ,ISYRES) 4759 ISYMI = MULD2H(ISYAI,ISYMA) 4760 ISYCKL = MULD2H(ISYMI,JSCKLI) 4761C 4762 IF (LWORK .LT. NRHF(ISYMI)*NRHF(ISYMJ)) THEN 4763 CALL QUIT('Insufficient memory in CCSDT_CY2O') 4764 END IF 4765C 4766 NTOCKL = MAX(NCKI(ISYCKL),1) 4767 NRHFI = MAX(NRHF(ISYMI),1) 4768C 4769 KOFF1 = ISAIKJ(ISYCKL,ISYMI) + 1 4770 KOFF2 = ISAIKJ(ISYCKL,ISYMJ) + 1 4771 KOFF3 = IMATIJ(ISYMI,ISYMJ) + 1 4772C 4773 4774 CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),NCKI(ISYCKL), 4775 * ONE,TMAT(KOFF1),NTOCKL,TROCC1(KOFF2),NTOCKL, 4776 * ONE,WORK(KOFF3),NRHFI) 4777C 4778C 4779 4780 600 CONTINUE 4781C 4782 IF (IPRINT .GT. 55) THEN 4783 XRMAT = DDOT(NMATIJ(ISYIJ),WORK,1,WORK,1) 4784 WRITE(LUPRI,*) '3.In CC3_W3_CY2O: Norm of RABIJ = ',XRMAT 4785 ENDIF 4786C 4787 CALL CC3_RABIJ(XI2,WORK,ISYMA,A,ISYMB,B,ISYRES) 4788C 4789 END IF 4790C 4791 END IF 4792C 4793 CALL QEXIT('CC3_W3_CY2O') 4794C 4795 RETURN 4796 END 4797C 4798C /* Deck cc3_eta_2 */ 4799 SUBROUTINE CC3_ETA_2(TMAT,ISTMAT,FOCKYCK,ISYFKY, 4800 * T2TP,ISYMT2,ETA2EFF,ISYRES,RAIJB, 4801 * INDSQ,LENSQ,INDAJLC,ISYAJL, 4802 * ISYMIB,IB,ISYMID,ID,WRK,LWRK) 4803C 4804C Calculate: 4805C 4806C tbar_mu3 <mu_3|[[X,E_aiE_bj],T_2]|HF> = 4807C 4808C - sum_ckdl ( tbar(aiblck) t(ckdl) X(jd) + tbar(aidjck) t(ckdl) X(lb) ) 4809C 4810C Written by P. Jorgensen and F. Pawlowski, Spring 2002. 4811C 4812 IMPLICIT NONE 4813C 4814 INTEGER ISTMAT,ISYFKY,ISYMT2,ISYRES,ISYAJL,ISYMIB,ISYMID 4815 INTEGER ISYMB,ISYMC,ISYMDKL,ISYMKLJ,ISYRMAT,ISYMJ,ISYMD 4816 INTEGER ISYMKL,ISYML,ISYMK,ISYMDK,ISYMLJ,ISYMKJ,ISYMBJ 4817 INTEGER ISYMAI,ISYMDLK,ISYMDC,ISYMKB,ISYMDL,ISYAIJ 4818 INTEGER IB,ID 4819 INTEGER NTOTD,NTOTK,NTOTAI,NTOTKL,NTOTB,NTOTAIJ 4820 INTEGER LENSQ,INDSQ(LENSQ,6),INDAJLC(*) 4821 INTEGER KOFF,KOFF1,KOFF2,KOFF3 4822 INTEGER KDKL,KKJL,KKLJ,KTMAT,KRMAT,KEND1,KT2KL,KKB 4823 INTEGER LWRK,LWRK1 4824 4825#if defined (SYS_CRAY) 4826 REAL TMAT(*),FOCKYCK(*),T2TP(*),ETA2EFF(*),RAIJB(*) 4827 REAL WRK(*) 4828 REAL ONE,XRMAT,DDOT 4829#else 4830 DOUBLE PRECISION TMAT(*),FOCKYCK(*),T2TP(*),ETA2EFF(*),RAIJB(*) 4831 DOUBLE PRECISION WRK(*) 4832 DOUBLE PRECISION ONE,XRMAT,DDOT 4833#endif 4834C 4835#include "priunit.h" 4836#include "ccsdsym.h" 4837#include "ccorb.h" 4838#include "ccsdinp.h" 4839C 4840 PARAMETER (ONE = 1.0D0) 4841C 4842 CALL QENTER('CC3_ETA_2') 4843C 4844C-------------------------- 4845C First contribution 4846C-------------------------- 4847C 4848C - tbar(aiblck) t(ckdl) X(jd) 4849C 4850C TMAT^BC(aikl) T2TP(dlkC) FOCKYCK(dj) 4851C 4852C T^C(dkl) FOCKYCK(dj) 4853C 4854C TF^C(kjl) 4855C 4856C TMAT^BC(aikl) G^C(klj) 4857C 4858 B = IB 4859 ISYMB = ISYMIB 4860 C = ID 4861 ISYMC = ISYMID 4862C 4863 ISYMDKL = MULD2H(ISYMT2,ISYMC) 4864 ISYMKLJ = MULD2H(ISYMDKL,ISYFKY) 4865 ISYRMAT = MULD2H(ISYRES,ISYMB) 4866C 4867 KDKL = 1 4868 KKJL = KDKL + NCKI(ISYMDKL) 4869 KKLJ = KKJL + NMAJIK(ISYMKLJ) 4870 KTMAT = KKLJ + NMAJIK(ISYMKLJ) 4871 KRMAT = KTMAT + NCKIJ(ISTMAT) 4872 KEND1 = KRMAT + NCKI(ISYRMAT) 4873 LWRK1 = LWRK - KEND1 4874C 4875 IF (LWRK1 .LT. 0) THEN 4876 CALL QUIT('Not enough space in cc3_eta_2 (1. contribution)') 4877 END IF 4878C 4879 CALL DZERO(WRK(KKJL),NMAJIK(ISYMKLJ)) 4880 CALL DZERO(WRK(KRMAT),NCKI(ISYRMAT)) 4881C 4882C T2TP(dlkC) put in WORK(dkl) 4883C 4884 KOFF = IT2SP(ISYMDKL,ISYMC) + NCKI(ISYMDKL)*(C - 1) + 1 4885 CALL CC_GATHER(NCKI(ISYMDKL),WRK(KDKL),T2TP(KOFF),INDAJLC) 4886C 4887C TF^C(kjl) = T^C(dkl) FOCKYCK(dj) 4888C 4889 DO ISYMJ = 1,NSYM 4890 ISYMD = MULD2H(ISYMJ,ISYFKY) 4891 ISYMKL = MULD2H(ISYMDKL,ISYMD) 4892 DO ISYML = 1,NSYM 4893 ISYMK = MULD2H(ISYMKL,ISYML) 4894 ISYMDK = MULD2H(ISYMD,ISYMK) 4895 ISYMKJ = MULD2H(ISYMJ,ISYMK) 4896 DO L = 1,NRHF(ISYML) 4897C 4898 KOFF1 = KDKL + ICKI(ISYMDK,ISYML) 4899 * + NT1AM(ISYMDK)*(L-1) 4900 * + IT1AM(ISYMD,ISYMK) 4901 KOFF2 = IT1AM(ISYMD,ISYMJ) + 1 4902 KOFF3 = KKJL + ISJIK(ISYMKJ,ISYML) 4903 * + NMATIJ(ISYMKJ)*(L-1) 4904 * + IMATIJ(ISYMK,ISYMJ) 4905C 4906 NTOTD = MAX(1,NVIR(ISYMD)) 4907 NTOTK = MAX(1,NRHF(ISYMK)) 4908C 4909 CALL DGEMM('T','N',NRHF(ISYMK),NRHF(ISYMJ), 4910 * NVIR(ISYMD),ONE,WRK(KOFF1),NTOTD, 4911 * FOCKYCK(KOFF2),NTOTD,ONE,WRK(KOFF3), 4912 * NTOTK) 4913 END DO 4914 END DO 4915 END DO 4916C 4917C G^C(klj) = TF^C(kjl) 4918C 4919 DO ISYMJ = 1,NSYM 4920 DO ISYML = 1,NSYM 4921 ISYMKJ = MULD2H(ISYMKLJ,ISYML) 4922 ISYMK = MULD2H(ISYMKJ,ISYMJ) 4923 ISYMKL = MULD2H(ISYMK,ISYML) 4924 DO J = 1,NRHF(ISYMJ) 4925 DO L = 1,NRHF(ISYML) 4926 DO K = 1,NRHF(ISYMK) 4927C 4928 KOFF1 = KKJL - 1 4929 * + ISJIK(ISYMKJ,ISYML) 4930 * + NMATIJ(ISYMKJ)*(L-1) 4931 * + IMATIJ(ISYMK,ISYMJ) 4932 * + NRHF(ISYMK)*(J-1) 4933 * + K 4934 KOFF2 = KKLJ - 1 4935 * + ISJIK(ISYMKL,ISYMJ) 4936 * + NMATIJ(ISYMKL)*(J-1) 4937 * + IMATIJ(ISYMK,ISYML) 4938 * + NRHF(ISYMK)*(L-1) 4939 * + K 4940C 4941 WRK(KOFF2) = WRK(KOFF1) 4942C 4943 END DO 4944 END DO 4945 END DO 4946 END DO 4947 END DO 4948C 4949 CALL DCOPY(LENSQ,TMAT,1,WRK(KTMAT),1) 4950C 4951C---------------------------------- 4952C Symmetry sorting if symmetry. 4953C---------------------------------- 4954C 4955 IF (NSYM .GT. 1) THEN 4956 CALL CC_GATHER(LENSQ,WRK(KTMAT),TMAT,INDSQ(1,6)) 4957 ENDIF 4958C 4959C---------------------------------- 4960C ETA2EFF(aibj) = ETA2EFF(aibj) + 4961C 4962C RAIJB1 = TMAT^BC(aikl) G^C(klj) 4963C---------------------------------- 4964C 4965 DO ISYMJ = 1,NSYM 4966C 4967 ISYMBJ = MULD2H(ISYMB,ISYMJ) 4968 ISYMAI = MULD2H(ISYMBJ,ISYRES) 4969 ISYMKL = MULD2H(ISTMAT,ISYMAI) 4970C 4971 NTOTAI = MAX(NT1AM(ISYMAI),1) 4972 NTOTKL = MAX(NMATIJ(ISYMKL),1) 4973 4974 KOFF1 = KTMAT + ISAIKL(ISYMAI,ISYMKL) 4975 KOFF2 = KKLJ + ISJIK(ISYMKL,ISYMJ) 4976 KOFF3 = KRMAT + ISAIK(ISYMAI,ISYMJ) 4977C 4978 CALL DGEMM('N','N',NT1AM(ISYMAI),NRHF(ISYMJ),NMATIJ(ISYMKL), 4979 * ONE,WRK(KOFF1),NTOTAI,WRK(KOFF2),NTOTKL, 4980 * ONE,WRK(KOFF3),NTOTAI) 4981C 4982 END DO 4983C 4984 IF (IPRINT .GT. 55) THEN 4985 XRMAT = DDOT(NCKI(ISYRMAT),WRK(KRMAT),1,WRK(KRMAT),1) 4986 WRITE(LUPRI,*) '1. In CC3_ETA_2: Norm of RMAT = ',XRMAT 4987 ENDIF 4988C 4989 CALL CC3_RACC(ETA2EFF,WRK(KRMAT),ISYMB,B,ISYRES) 4990C 4991C---------------------------------- 4992C Second contribution. 4993C---------------------------------- 4994C 4995C TMAT^DC(aikj) T2TP(DlkC) FOCKYCK(bl) 4996C 4997C T^DC(kl) FOCKYCK(bl) 4998C 4999C WORK^DC(aijk) G^DC(kb) 5000C 5001 D = IB 5002 ISYMD = ISYMIB 5003 C = ID 5004 ISYMC = ISYMID 5005C 5006 ISYMDLK = MULD2H(ISYMT2,ISYMC) 5007 ISYMDC = MULD2H(ISYMD,ISYMC) 5008 ISYMKL = MULD2H(ISYMDC,ISYMT2) 5009 ISYMKB = MULD2H(ISYMKL,ISYFKY) 5010C 5011 KT2KL = 1 5012 KKB = KT2KL + NMATIJ(ISYMKL) 5013 KTMAT = KKB + NT1AM(ISYMKB) 5014 KEND1 = KTMAT + NCKIJ(ISTMAT) 5015 LWRK1 = LWRK - KEND1 5016C 5017 IF (LWRK1 .LT. 0) THEN 5018 CALL QUIT('Not enough space in cc3_eta_2 (2. contribution)') 5019 END IF 5020C 5021 CALL DZERO(WRK(KKB),NT1AM(ISYMKB)) 5022C 5023C ------------------------------- 5024C T2TP(DlkC) sorted as T^{DC}_{kl} equal T_{kl} 5025C ------------------------------- 5026C 5027 DO ISYML = 1,NSYM 5028 ISYMK = MULD2H(ISYMKL,ISYML) 5029 ISYMDL = MULD2H(ISYMDLK,ISYMK) 5030 DO L = 1,NRHF(ISYML) 5031 DO K = 1,NRHF(ISYMK) 5032 KOFF1 = IT2SP(ISYMDLK,ISYMC) 5033 * + NCKI(ISYMDLK)*(C-1) 5034 * + ICKI(ISYMDL,ISYMK) 5035 * + NT1AM(ISYMDL)*(K-1) 5036 * + IT1AM(ISYMD,ISYML) 5037 * + NVIR(ISYMD)*(L-1) 5038 * + D 5039 KOFF2 = IMATIJ(ISYMK,ISYML) 5040 * + NRHF(ISYMK)*(L-1) 5041 * + K 5042 WRK(KT2KL-1+KOFF2) = T2TP(KOFF1) 5043 ENDDO 5044 ENDDO 5045 ENDDO 5046C 5047C TF^DC(kb) = T^DC(kl) FOCKYCK(bl) 5048C 5049 DO ISYML = 1,NSYM 5050 ISYMK = MULD2H(ISYMKL,ISYML) 5051 ISYMB = MULD2H(ISYFKY,ISYML) 5052C 5053 KOFF1 = KT2KL + IMATIJ(ISYMK,ISYML) 5054 KOFF2 = IT1AM(ISYMB,ISYML) + 1 5055 KOFF3 = KKB + IT1AMT(ISYMK,ISYMB) 5056C 5057 NTOTK = MAX(1,NRHF(ISYMK)) 5058 NTOTB = MAX(1,NVIR(ISYMB)) 5059C 5060 CALL DGEMM('N','T',NRHF(ISYMK),NVIR(ISYMB), 5061 * NRHF(ISYML),ONE,WRK(KOFF1),NTOTK, 5062 * FOCKYCK(KOFF2),NTOTB,ONE,WRK(KOFF3), 5063 * NTOTK) 5064C 5065 END DO 5066C 5067C Sort TMAT^DC(aikj) as WORK^DC(aijk) 5068C 5069 DO I = 1, NCKIJ(ISTMAT) 5070 WRK(KTMAT-1+I) = TMAT(INDSQ(I,3)) 5071 ENDDO 5072C 5073C---------------------------------- 5074C 5075C RAIJB1 = WORK^DC(aijk) TF^DC(kb) 5076C 5077C---------------------------------- 5078C 5079 DO ISYMB = 1,NSYM 5080C 5081 ISYMK = MULD2H(ISYMB,ISYMKB) 5082 ISYAIJ = MULD2H(ISTMAT,ISYMK) 5083C 5084 KOFF1 = KTMAT + ISAIKJ(ISYAIJ,ISYMK) 5085 KOFF2 = KKB + IT1AMT(ISYMK,ISYMB) 5086 KOFF3 = IT2SP(ISYAIJ,ISYMB) + 1 5087C 5088 NTOTAIJ = MAX(NCKI(ISYAIJ),1) 5089 NTOTK = MAX(NRHF(ISYMK),1) 5090C 5091 CALL DGEMM('N','N',NCKI(ISYAIJ),NVIR(ISYMB),NRHF(ISYMK), 5092 * ONE,WRK(KOFF1),NTOTAIJ,WRK(KOFF2),NTOTK, 5093 * ONE,RAIJB(KOFF3),NTOTAIJ) 5094C 5095 END DO 5096C 5097 IF (IPRINT .GT. 55) THEN 5098 XRMAT = DDOT(NT2SQ(ISYRES),RAIJB,1,RAIJB,1) 5099 WRITE(LUPRI,*) '2. In CC3_ETA_2: Norm of RAIJB = ',XRMAT 5100 ENDIF 5101C 5102 CALL QEXIT('CC3_ETA_2') 5103C 5104 RETURN 5105 END 5106C 5107C /* Deck cc3_eta_1 */ 5108 SUBROUTINE CC3_ETA_1(FOCKYCK,ISYFKY,ETA1EFF,ISYRES, 5109 * WRK,LWRK) 5110C 5111C--------------------------------------------------------------------- 5112C 5113C Calculate: 5114C 5115C <L3|[[X,tau_ai],T_3]|HF> = 5116C 5117C sum_l X(la) D^0(li) - sum_d X(id) D^0(ad) 5118C 5119C where 5120C D^0(pq) = <L3|[E_pq,T_3]|HF> ; p,q = a,b or i,j 5121C 5122C 5123 IMPLICIT NONE 5124C 5125 CHARACTER FNDPTIJ*5, FNDPTAB*5 5126 PARAMETER(FNDPTIJ = 'DPTIJ', FNDPTAB='DPTAB') 5127C 5128 INTEGER ISYFKY,ISYRES,ISYM0,ISYMI,ISYMA,ISYML,ISYMD 5129 INTEGER KD0AB,KD0IJ,KEND1 5130 INTEGER KOFF1,KOFF2,KOFF3 5131 INTEGER NTOTL,NTOTA,NTOTD 5132 INTEGER LWRK,LWRK1 5133 INTEGER LUPTIJ, LUPTAB 5134C 5135#if defined (SYS_CRAY) 5136 REAL FOCKYCK(*),ETA1EFF(*),WRK(*),ONE 5137#else 5138 DOUBLE PRECISION FOCKYCK(*),ETA1EFF(*),WRK(*),ONE 5139#endif 5140C 5141#include "priunit.h" 5142#include "ccsdsym.h" 5143#include "ccorb.h" 5144C 5145 PARAMETER (ONE = 1.0D0) 5146C 5147 CALL QENTER('CC3_ETA_1') 5148C 5149 ISYM0 = 1 5150C 5151 KD0AB = 1 5152 KD0IJ = KD0AB + NMATAB(ISYM0) 5153 KEND1 = KD0IJ + NMATIJ(ISYM0) 5154 LWRK1 = LWRK - KEND1 5155C 5156 IF (LWRK1 .LT. 0) THEN 5157 CALL QUIT('Not enough space in cc3_eta_1') 5158 END IF 5159C 5160 LUPTIJ = -1 5161 LUPTAB = -1 5162C 5163 ! read intermediates from ground state density from file... 5164C 5165 CALL WOPEN2(LUPTIJ,FNDPTIJ,64,0) 5166 CALL GETWA2(LUPTIJ,FNDPTIJ,WRK(KD0IJ),1,NMATIJ(ISYM0)) 5167 CALL WCLOSE2(LUPTIJ,FNDPTIJ,'KEEP') 5168 5169 CALL WOPEN2(LUPTAB,FNDPTAB,64,0) 5170 CALL GETWA2(LUPTAB,FNDPTAB,WRK(KD0AB),1,NMATAB(ISYM0)) 5171 CALL WCLOSE2(LUPTAB,FNDPTAB,'KEEP') 5172C 5173C ---------------------------------------------------------------- 5174C sum_l X(la) D^0(li) 5175C FOCKYCK(al) 5176C ---------------------------------------------------------------- 5177C 5178 DO ISYMI = 1, NSYM 5179 ISYMA = MULD2H(ISYRES,ISYMI) 5180 ISYML = MULD2H(ISYFKY,ISYMA) 5181C 5182 KOFF1 = IT1AM(ISYMA,ISYML) + 1 5183 KOFF2 = KD0IJ + IMATIJ(ISYML,ISYMI) 5184 KOFF3 = IT1AM(ISYMA,ISYMI) + 1 5185C 5186 NTOTL = MAX(NRHF(ISYML),1) 5187 NTOTA = MAX(NVIR(ISYMA),1) 5188C 5189 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI),NRHF(ISYML), 5190 * ONE,FOCKYCK(KOFF1),NTOTA,WRK(KOFF2),NTOTL, 5191 * ONE,ETA1EFF(KOFF3),NTOTA) 5192C 5193 END DO 5194C 5195C 5196C ---------------------------------------------------------------- 5197C - sum_d D^0(ad) X(id) 5198C FOCKYCK(di) 5199C ---------------------------------------------------------------- 5200C 5201 DO ISYMI = 1, NSYM 5202 ISYMA = MULD2H(ISYRES,ISYMI) 5203 ISYMD = MULD2H(ISYFKY,ISYMI) 5204C 5205 KOFF1 = KD0AB + IMATAB(ISYMA,ISYMD) 5206 KOFF2 = IT1AM(ISYMD,ISYMI) + 1 5207 KOFF3 = IT1AM(ISYMA,ISYMI) + 1 5208C 5209 NTOTA = MAX(NVIR(ISYMA),1) 5210 NTOTD = MAX(NVIR(ISYMD),1) 5211C 5212 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI),NVIR(ISYMD), 5213 * -ONE,WRK(KOFF1),NTOTA,FOCKYCK(KOFF2),NTOTD, 5214 * ONE,ETA1EFF(KOFF3),NTOTA) 5215C 5216 END DO 5217C 5218 CALL QEXIT('CC3_ETA_1') 5219C 5220 RETURN 5221 END 5222C 5223 5224