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 CCSDT_ETA_CONT(LISTL,LISTR,NODDY, 21 & IDOTS,DOTPROD, 22 & NXETRAN,IXETRAN,MXVEC, 23 & NLEFT,IEASTEND, 24 & NDENS,IDXLVEC,IDXRVEC, 25 & WORK,LWORK, 26 & FNDELD,FNCKJD,FNDKBC,FNTOC, 27 & FN3VI,FNDKBC3,FN3FOPX,FN3FOP2X) 28*---------------------------------------------------------------------* 29* 30* Purpose: just a little driver to call the real routines... 31* 32* Written by Christof Haettig, Mai 2002, Friedrichstal 33* 34*=====================================================================* 35 IMPLICIT NONE 36#include "ccsdsym.h" 37#include "priunit.h" 38#include "ccorb.h" 39#include "dummy.h" 40#include "cclists.h" 41#include "ccroper.h" 42 43 LOGICAL LOCDBG 44 PARAMETER (LOCDBG = .FALSE.) 45 46 CHARACTER*12 FNDEFF 47 PARAMETER( FNDEFF = 'TMP-XIETADEN') 48 49 INTEGER ISYM0 50 PARAMETER( ISYM0 = 1 ) 51 52 CHARACTER*3 LISTL, LISTR 53 LOGICAL NODDY 54 INTEGER LWORK, MXVEC, NXETRAN, NLEFT, NDENS 55 INTEGER IDOTS(MXVEC,NXETRAN), IEASTEND(2,NLEFT), 56 & IDXLVEC(NDENS), IDXRVEC(NDENS), 57 & IXETRAN(MXDIM_XEVEC,NXETRAN) 58 CHARACTER*(*) FNDELD, FNCKJD, FNDKBC, FNTOC, FN3VI, FNDKBC3 59 CHARACTER*(*) FN3FOPX, FN3FOP2X 60 61#if defined (SYS_CRAY) 62 REAL DOTPROD(MXVEC,NXETRAN), TRIPCON, CC_XIETA_DENCON 63 REAL WORK(LWORK) 64#else 65 DOUBLE PRECISION DOTPROD(MXVEC,NXETRAN), TRIPCON, CC_XIETA_DENCON 66 DOUBLE PRECISION WORK(LWORK) 67#endif 68 69 CHARACTER MODEL*(10), LABELH*(8) 70 LOGICAL SKIPXI, SKIPETA 71 INTEGER LUDELD, LUCKJD, LUDKBC, LUTOC, LU3VI, LUDKBC3, LU3FOPX, 72 * LU3FOP2X, LUDEFF, LUFOCK, 73 * KLAMP0, KLAMH0, KFOCK0, KEND0, LWRK0, KT1AMP0, KEND1, 74 * LWRK1, IOPT, ILEFT, IDLSTL, IDENS, IDLSTR, ISYMR, ISYML, 75 * ISYDEN, KDAB, KDIJ, KDIA, IADRIJ, IADRAB, IADRIA, IOPER, 76 * ISYHOP, ISYCTR, ISYETA, IVEC, IDX, M2BST, ITRAN, ISYM 77C 78 79 INTEGER ILSTSYM 80 81*---------------------------------------------------------------------* 82* begin 83*---------------------------------------------------------------------* 84 CALL QENTER('CCSDT_ETA_CONT') 85 86 IF (LOCDBG) THEN 87 WRITE(LUPRI,*) 'entered CCSDT_ETA_CONT... NODDY=',NODDY 88 WRITE(LUPRI,*) 'LISTL,LISTR:',LISTL,LISTR 89 END IF 90 91 LUDEFF = -1 92 CALL WOPEN2(LUDEFF,FNDEFF,64,0) 93 94 M2BST = 0 95 DO ISYM = 1, NSYM 96 M2BST = MAX(M2BST,N2BST(ISYM)) 97 END DO 98*---------------------------------------------------------------------* 99* initialize 0.th-order Lambda and Fock matrices: 100*---------------------------------------------------------------------* 101 KLAMP0 = 1 102 KLAMH0 = KLAMP0 + NLAMDT 103 KFOCK0 = KLAMH0 + NLAMDT 104 KEND0 = KFOCK0 + N2BST(ISYM0) 105 LWRK0 = LWORK - KEND0 106 107 KT1AMP0 = KEND0 108 KEND1 = KT1AMP0 + NT1AMX 109 LWRK1 = LWORK - KEND1 110 111 IF (LWRK1.LT.0) CALL QUIT('Out of memory in CCSDT_ETA_CONT') 112 113 IOPT = 1 114 CALL CC_RDRSP('R0',0,ISYM0,IOPT,MODEL,WORK(KT1AMP0),DUMMY) 115 116 CALL LAMMAT(WORK(KLAMP0),WORK(KLAMH0),WORK(KT1AMP0), 117 & WORK(KEND1),LWRK1) 118 119* read zeroth-order AO Fock matrix from file: 120 LUFOCK = -1 121 CALL GPOPEN(LUFOCK,'CC_FCKH','OLD',' ','UNFORMATTED',IDUMMY, 122 & .FALSE.) 123 REWIND(LUFOCK) 124 READ(LUFOCK) (WORK(KFOCK0-1+I),I=1,N2BST(ISYM0)) 125 CALL GPCLOSE(LUFOCK,'KEEP') 126 127 CALL CC_FCKMO(WORK(KFOCK0),WORK(KLAMP0),WORK(KLAMH0), 128 & WORK(KEND1),LWRK1,1,1,1) 129 130*---------------------------------------------------------------------* 131* precompute effective densities and store on file: 132*---------------------------------------------------------------------* 133 LUDELD = -1 134 LUCKJD = -1 135 LUDKBC = -1 136 LUTOC = -1 137 LU3VI = -1 138 LUDKBC3 = -1 139 LU3FOPX = -1 140 LU3FOP2X= -1 141 142 CALL WOPEN2(LUDELD,FNDELD,64,0) 143 CALL WOPEN2(LUCKJD,FNCKJD,64,0) 144 CALL WOPEN2(LUDKBC,FNDKBC,64,0) 145 CALL WOPEN2(LUTOC,FNTOC,64,0) 146 CALL WOPEN2(LU3VI,FN3VI,64,0) 147 CALL WOPEN2(LUDKBC3,FNDKBC3,64,0) 148 CALL WOPEN2(LU3FOPX,FN3FOPX,64,0) 149 CALL WOPEN2(LU3FOP2X,FN3FOP2X,64,0) 150 151 DO ILEFT = 1, NLEFT 152 IDLSTL = IDXLVEC(IEASTEND(1,ILEFT)) 153 154 DO IDENS = IEASTEND(1,ILEFT), IEASTEND(2,ILEFT) 155 156 IF (IDXLVEC(IDENS).NE.IDLSTL) 157 * CALL QUIT('Loop error in CCSDT_ETA_CONT.') 158 159 IDLSTR = IDXRVEC(IDENS) 160 ISYMR = ILSTSYM(LISTR,IDLSTR) 161 ISYML = ILSTSYM(LISTL,IDLSTL) 162 ISYDEN = MULD2H(ISYMR,ISYML) 163 164 KDAB = KEND0 165 KDIJ = KDAB + NMATAB(ISYDEN) 166 KDIA = KDIJ + NMATIJ(ISYDEN) 167 KEND1 = KDIA + NT1AM(ISYDEN) 168 LWRK1 = LWORK - KEND1 169 170 IF (LWRK1.LT.0) CALL QUIT('Out of memory in CCSDT_ETA_CONT') 171 172 IF (LISTL.EQ.'L0 ') THEN 173C -------------------- 174C Compute the density: 175C -------------------- 176 CALL CCSDT_ETA_DEN(LISTL,IDLSTL,LISTR,IDLSTR, 177 * WORK(KLAMP0),WORK(KLAMH0),WORK(KFOCK0), 178 * WORK(KDIJ),WORK(KDAB),WORK(KDIA), 179 * WORK(KEND1),LWRK1, 180 * LUDELD,FNDELD,LUCKJD,FNCKJD,LUDKBC, 181 * FNDKBC,LUTOC,FNTOC,LU3VI,FN3VI, 182 * LUDKBC3,FNDKBC3,LU3FOPX,FN3FOPX, 183 * LU3FOP2X,FN3FOP2X) 184 ELSE 185 186 IF (NODDY) THEN 187 CALL CCSDT_ADEN_NODDY(LISTL,IDLSTL,LISTR,IDLSTR, 188 * WORK(KLAMP0),WORK(KLAMH0),WORK(KFOCK0), 189 * WORK(KDIJ),WORK(KDAB),.TRUE.,WORK(KDIA), 190 * .FALSE.,DUMMY,WORK(KEND1),LWRK1, 191 * LUDELD,FNDELD,LUCKJD,FNCKJD,LUDKBC, 192 * FNDKBC,LUTOC,FNTOC,LU3VI,FN3VI, 193 * LUDKBC3,FNDKBC3,LU3FOPX,FN3FOPX, 194 * LU3FOP2X,FN3FOP2X) 195 ELSE 196 IF (LISTR(1:3).EQ.'R1 '.OR.LISTR(1:3).EQ.'RE ') THEN 197 CALL CC3_ADEN(LISTL,IDLSTL,LISTR,IDLSTR, 198 * WORK(KLAMP0),WORK(KLAMH0),WORK(KFOCK0), 199 * WORK(KDIJ),WORK(KDAB), 200 * .TRUE.,WORK(KDIA),ISYDEN, 201 * .FALSE.,DUMMY, 202 * WORK(KEND1),LWRK1, 203 * LUDELD,FNDELD,LUCKJD,FNCKJD,LUDKBC, 204 * FNDKBC,LUTOC,FNTOC,LU3VI,FN3VI, 205 * LUDKBC3,FNDKBC3,LU3FOPX,FN3FOPX, 206 * LU3FOP2X,FN3FOP2X) 207 ELSEIF (LISTR(1:3).EQ.'R2 '.OR.LISTR(1:3).EQ.'ER1') THEN 208 CALL CC3_ADEN_CUB(LISTL,IDLSTL,LISTR,IDLSTR, 209 * WORK(KLAMP0),WORK(KLAMH0),WORK(KFOCK0), 210 * WORK(KDIJ),WORK(KDAB),WORK(KDIA),ISYDEN, 211 * WORK(KEND1),LWRK1, 212 * LUDELD,FNDELD,LUCKJD,FNCKJD,LUDKBC, 213 * FNDKBC,LUTOC,FNTOC,LU3VI,FN3VI, 214 * LUDKBC3,FNDKBC3,LU3FOPX,FN3FOPX, 215 * LU3FOP2X,FN3FOP2X) 216 ELSE 217 CALL QUIT('Unknown LISTR in CCSDT_ETA_CONT.') 218 ENDIF 219 220 END IF 221 222 END IF 223 224 225 IADRIJ = 1 + M2BST*(IDENS-1) 226 CALL PUTWA2(LUDEFF,FNDEFF,WORK(KDIJ),IADRIJ,NMATIJ(ISYDEN)) 227 228 IADRAB = IADRIJ + NMATIJ(ISYDEN) 229 CALL PUTWA2(LUDEFF,FNDEFF,WORK(KDAB),IADRAB,NMATAB(ISYDEN)) 230 231 IADRIA = IADRAB + NMATAB(ISYDEN) 232 CALL PUTWA2(LUDEFF,FNDEFF,WORK(KDIA),IADRIA,NT1AM(ISYDEN)) 233 234 END DO ! IDENS 235 236 END DO ! ILEFT 237 238 CALL WCLOSE2(LUDELD,FNDELD,'KEEP') 239 CALL WCLOSE2(LUCKJD,FNCKJD,'KEEP') 240 CALL WCLOSE2(LUTOC,FNTOC,'KEEP') 241 CALL WCLOSE2(LUDKBC,FNDKBC,'KEEP') 242 CALL WCLOSE2(LU3VI,FN3VI,'KEEP') 243 CALL WCLOSE2(LUDKBC3,FNDKBC3,'KEEP') 244 CALL WCLOSE2(LU3FOPX,FN3FOPX,'KEEP') 245 CALL WCLOSE2(LU3FOP2X,FN3FOP2X,'KEEP') 246 247*---------------------------------------------------------------------* 248* compute from the densities the polarizability contributions: 249*---------------------------------------------------------------------* 250 DO ITRAN = 1, NXETRAN 251 IOPER = IXETRAN(1,ITRAN) ! operator index 252 IDLSTL = IXETRAN(2,ITRAN) ! Left vector index 253 254 ISYHOP = ISYOPR(IOPER) ! symmetry of O operator 255 LABELH = LBLOPR(IOPER) ! operator label 256 257 SKIPXI = ( IXETRAN(3,ITRAN) .EQ. -1 ) 258 SKIPETA = ( IXETRAN(4,ITRAN) .EQ. -1 ) 259 260 IF (.NOT.SKIPETA) THEN 261 ISYCTR = ILSTSYM(LISTL,IDLSTL) ! sym. of Left (ZETA) vector 262 ISYETA = MULD2H(ISYHOP,ISYCTR) ! sym. of ETA result vect. 263 264 IVEC = 1 265 DO WHILE(IDOTS(IVEC,ITRAN).NE.0 .AND. IVEC.LE.MXVEC) 266 IDLSTR = IDOTS(IVEC,ITRAN) 267 ISYMR = ILSTSYM(LISTR,IDLSTR) 268 269 IF (ISYMR.EQ.ISYETA) THEN 270 ! find index of density 271 IDENS = 0 272 DO IDX = 1, NDENS 273 IF (IDLSTL.EQ.IDXLVEC(IDX) .AND. 274 & IDLSTR.EQ.IDXRVEC(IDX) ) THEN 275 IDENS = IDX 276 END IF 277 END DO 278 279 IF (IDENS.EQ.0) 280 & CALL QUIT('Density not found in ccsdt_eta_cont.') 281 282 TRIPCON = CC_XIETA_DENCON(IDENS,LABELH,ISYHOP, 283 & WORK(KLAMP0),WORK(KLAMH0), 284 & LUDEFF,FNDEFF,M2BST,WORK(KEND0),LWRK0) 285 286 IF (LOCDBG) THEN 287 WRITE(LUPRI,*) 'IVEC,ITRAN,TRIPCON:', 288 & IVEC,ITRAN,TRIPCON 289 WRITE(LUPRI,*) 'DOTPROD before,after:', 290 & DOTPROD(IVEC,ITRAN),DOTPROD(IVEC,ITRAN)+TRIPCON 291 END IF 292 293 DOTPROD(IVEC,ITRAN) = DOTPROD(IVEC,ITRAN) + TRIPCON 294 295 END IF 296 297 IVEC = IVEC + 1 298 END DO 299 300 END IF 301 END DO ! ITRAN 302 303 304 CALL WCLOSE2(LUDEFF,FNDEFF,'DELETE') 305 306 CALL QEXIT('CCSDT_ETA_CONT') 307 RETURN 308 END 309*=====================================================================* 310*=====================================================================* 311 SUBROUTINE CCSDT_ETA_DEN(LISTL,IDLSTL,LISTR,IDLSTR, 312 * XLAMDP,XLAMDH,FOCK0, 313 * DENIJ,DENAB,DENIA, 314 * WORK,LWORK, 315 * LUDELD,FNDELD,LUCKJD, FNCKJD, 316 * LUDKBC,FNDKBC,LUTOC, FNTOC, 317 * LU3VI, FN3VI, LUDKBC3,FNDKBC3, 318 * LU3FOP,FN3FOP,LU3FOP2,FN3FOP2) 319*---------------------------------------------------------------------* 320* 321* Purpose: compute triples contribution to a contraction of 322* an Eta vector with some 'right' vector 323* 324* Calculation is done via a Eta-density: 325* 326* Eta(X) x T^Y = sum_pq D^eta_pq X_pq 327* 328* Written by Christof Haettig, Mai 2002, Aarhus 329* 330*********************************************************************** 331* 332* Spring 2004, Filip Pawlowski, Aarhus: Modified to treat also the 333* Cauchy moments. 334* 335*=====================================================================* 336 IMPLICIT NONE 337#include "priunit.h" 338#include "dummy.h" 339#include "ccsdsym.h" 340#include "ccorb.h" 341#include "ccr1rsp.h" 342#include "ccexci.h" 343#include "ccrc1rsp.h" 344C 345 LOGICAL LOCDBG 346 LOGICAL LCAUCHY 347 PARAMETER(LOCDBG = .FALSE.) 348C 349 INTEGER ISYM0 350 PARAMETER(ISYM0 = 1) 351C 352 CHARACTER*10 FNT2Y 353 PARAMETER(FNT2Y = 'CCT2Y__TMP' ) 354C 355 CHARACTER FNDPTIJ*5, FNDPTAB*5 356 PARAMETER(FNDPTIJ = 'DPTIJ', FNDPTAB='DPTAB') 357C 358 CHARACTER LISTL*3, LISTR*3 359 CHARACTER*(*) FNDELD, FNCKJD, FNDKBC, FNTOC, FN3VI, FNDKBC3 360 CHARACTER*(*) FN3FOP, FN3FOP2 361 362 INTEGER LUDELD, LUCKJD, LUDKBC, LUTOC, LU3VI, LUDKBC3, LU3FOP, 363 * LU3FOP2, IDLSTR, IDLSTL 364 365 INTEGER ISYCTR, ISTAMY, KT2Y, ISYDEN, KDAB, KDIJ, KDIA1, 366 * KDIA2, KDIA3, KDIA4, KMMAT, KYMMAT, KT2TP, KL1AM, KL2TP 367 INTEGER LUPTIJ, LUPTAB, LUT2Y 368 369 CHARACTER LABELY*8, MODEL*10, CDUMMY*1 370 INTEGER ISYM, ILEN, ISYMIM, NIMFN(8), IOPT, KEND0, LWRK0, 371 & KD0AB, KD0IJ, KT1AMY, KEND2, LWRK2, IOPTT2, ISYMA, 372 & ISYMI, ISYMB, ISYMJ, KOFFDBA, KOFFDIJ, KOFFDIA, KOFFTBI, 373 & KOFFTAJ, NVIRA, NVIRB, NRHFI, ISYMFN, LWORK 374C 375 INTEGER NCAU 376 377#if defined (SYS_CRAY) 378 REAL XLAMDP(*), XLAMDH(*), FOCK0(*) 379 REAL PRECISION FREQY, WORK(*), DDOT 380 REAL ZERO, ONE, TWO, HALF 381 REAL DENIJ(*),DENAB(*),DENIA(*) 382#else 383 DOUBLE PRECISION XLAMDP(*), XLAMDH(*), FOCK0(*) 384 DOUBLE PRECISION FREQY, WORK(*), DDOT 385 DOUBLE PRECISION ZERO, ONE, TWO, HALF 386 DOUBLE PRECISION DENIJ(*),DENAB(*),DENIA(*) 387#endif 388C 389 PARAMETER (ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, TWO = 2.0D0) 390C 391 INTEGER ILSTSYM 392C 393C------------------------------------------------------------ 394C some initializations: 395C------------------------------------------------------------ 396C 397 CALL QENTER('CTETADEN') 398 399 IF (LISTL(1:3).EQ.'L0 ') THEN 400 ISYCTR = ISYM0 401 ELSE 402 ! ups, probably higher-order response, not yet implemented here 403 CALL QUIT('Unknown type of left vector in CCSDT_ETA_DEN.') 404 END IF 405C 406 !initialize LCAUCHY and NCAU 407 LCAUCHY = .FALSE. 408 NCAU = 0 409 410 IF (LISTR(1:3).EQ.'R1 ') THEN 411 ISTAMY = ISYLRT(IDLSTR) 412 FREQY = FRQLRT(IDLSTR) 413 LABELY = LRTLBL(IDLSTR) 414 415 IF (LORXLRT(IDLSTR)) 416 & CALL QUIT('NO ORBITAL RELAXED PERTURBATION IN CCSDT_ETA_DEN') 417 418 ELSE IF (LISTR(1:3).EQ.'RE ') THEN 419 ISTAMY = ILSTSYM(LISTR,IDLSTR) 420 FREQY = EIGVAL(IDLSTR) 421 LABELY = '- none -' 422 423 ELSE IF (LISTR(1:3).EQ.'RC ') THEN 424 ISTAMY = ISYLRC(IDLSTR) 425 FREQY = ZERO 426 LABELY = LRCLBL(IDLSTR) 427 NCAU = ILRCAU(IDLSTR) 428 LCAUCHY = .TRUE. 429 430 ELSE 431 ! ups, probably higher-order response, not yet implemented here 432 WRITE(LUPRI,*)'LISTR = ',LISTR 433 CALL QUIT('Unknown type of right vector in CCSDT_ETA_DEN.') 434 END IF 435C 436C------------------------------------------------------------ 437C compute dimensions of M matrices NIMFN: 438C------------------------------------------------------------ 439 440 DO ISYM = 1, NSYM 441 ILEN = 0 442 DO ISYMFN = 1, NSYM 443 ISYMIM = MULD2H(ISYM,ISYMFN) 444 ILEN = ILEN + NMATIJ(ISYMIM)*NT1AM(ISYMFN) 445 END DO 446 NIMFN(ISYM) = ILEN 447 END DO 448C 449C------------------------------------------------------------------ 450C resort response amplitudes and dump them on a temporary file: 451C------------------------------------------------------------------ 452C 453 KT2Y = 1 454 KEND0 = KT2Y + NT2SQ(ISTAMY) 455 LWRK0 = LWORK - KEND0 456 457 IF (LWRK0 .LT. NT2AM(ISTAMY)) THEN 458 CALL QUIT('Not enough memory for squared T2Y (CCSDT_ETA_DEN)') 459 ENDIF 460 461 IOPT = 2 462 CALL CC_RDRSP(LISTR,IDLSTR,ISTAMY,IOPT,MODEL,DUMMY,WORK(KEND0)) 463 464 Call CCLR_DIASCL(WORK(KEND0),TWO,ISTAMY) 465 466 CALL CC_T2SQ(WORK(KEND0),WORK(KT2Y),ISTAMY) 467 468 LUT2Y = -1 469 CALL WOPEN2(LUT2Y,FNT2Y,64,0) 470 CALL PUTWA2(LUT2Y,FNT2Y,WORK(KT2Y),1,NT2SQ(ISTAMY)) 471C 472C------------------------------------------------------- 473C initial allocations, prepare T2 and L2 amplitudes: 474C------------------------------------------------------- 475C 476 ISYDEN = MULD2H(ISYCTR,ISTAMY) 477 478 KEND0 = 1 479 480 KDAB = KEND0 481 KDIJ = KDAB + NMATAB(ISYDEN) 482 KDIA1 = KDIJ + NMATIJ(ISYDEN) 483 KDIA2 = KDIA1 + NT1AM(ISYDEN) 484 KDIA3 = KDIA2 + NT1AM(ISYDEN) 485 KDIA4 = KDIA3 + NT1AM(ISYDEN) 486 KEND0 = KDIA4 + NT1AM(ISYDEN) 487 488 KMMAT = KEND0 489 KYMMAT = KMMAT + NIMFN(ISYCTR) 490 KEND0 = KYMMAT + NIMFN(ISYDEN) 491 492 KT2TP = KEND0 493 KL1AM = KT2TP + NT2SQ(ISYM0) 494 KL2TP = KL1AM + NT1AM(ISYCTR) 495 KEND0 = KL2TP + NT2SQ(ISYCTR) 496 497 LWRK0 = LWORK - KEND0 498 499 IF (LWRK0.LT.NT2SQ(ISYM0) .OR. LWRK0.LT.NT2SQ(ISYCTR)) THEN 500 CALL QUIT('Not enough memory to construct T2TP (CCSDT_ETA_DEN)') 501 ENDIF 502 503 504 IOPT = 2 505 CALL CC_RDRSP('R0',0,ISYM0,IOPT,MODEL,DUMMY,WORK(KT2TP)) 506 507 CALL CC_T2SQ(WORK(KT2TP),WORK(KEND0),ISYM0) 508 509 CALL CC3_T2TP(WORK(KT2TP),WORK(KEND0),ISYM0) 510 511 IF (LOCDBG) WRITE(LUPRI,*) 'Norm of T2TP ', 512 * DDOT(NT2SQ(ISYM0),WORK(KT2TP),1,WORK(KT2TP),1) 513 514 IOPT = 3 515 CALL CC_RDRSP(LISTL,IDLSTL,ISYCTR,IOPT,MODEL, 516 * WORK(KL1AM),WORK(KL2TP)) 517 518 CALL CC_T2SQ(WORK(KL2TP),WORK(KEND0),ISYCTR) 519 520 CALL CC3_T2TP(WORK(KL2TP),WORK(KEND0),ISYCTR) 521 522 IF (LOCDBG) WRITE(LUPRI,*) 'Norm of L2TP ', 523 * DDOT(NT2SQ(ISYCTR),WORK(KL2TP),1,WORK(KL2TP),1) 524 525 ! initialize density matrices 526 CALL DZERO(WORK(KDAB), NMATAB(ISYDEN)) 527 CALL DZERO(WORK(KDIJ), NMATIJ(ISYDEN)) 528 CALL DZERO(WORK(KDIA1),NT1AM(ISYDEN)) 529 CALL DZERO(WORK(KDIA2),NT1AM(ISYDEN)) 530 CALL DZERO(WORK(KDIA3),NT1AM(ISYDEN)) 531 CALL DZERO(WORK(KDIA4),NT1AM(ISYDEN)) 532 533 ! initialize M-matrix intermediates 534 CALL DZERO(WORK(KMMAT), NIMFN(ISYCTR)) 535 CALL DZERO(WORK(KYMMAT),NIMFN(ISYDEN)) 536 537 !frequency of LISTL is handled inside 538 CALL CCSDT_ETA_FA_DEN(LISTL,IDLSTL,ISYCTR, 539 * LISTR,IDLSTR,ISTAMY,FREQY,LABELY, 540 * LCAUCHY,NCAU, 541 * XLAMDP,XLAMDH,FOCK0, 542 * WORK(KT2TP),WORK(KL2TP),WORK(KL1AM), 543 * ISYDEN,WORK(KDAB),WORK(KDIJ), 544 * .TRUE.,WORK(KDIA3), 545 * .TRUE.,WORK(KMMAT),.TRUE.,WORK(KYMMAT), 546 * WORK(KEND0),LWRK0, 547 * LUDELD,FNDELD,LUCKJD, FNCKJD, 548 * LUDKBC,FNDKBC,LUTOC, FNTOC, 549 * LU3VI, FN3VI, LUDKBC3,FNDKBC3, 550 * LU3FOP,FN3FOP,LU3FOP2,FN3FOP2, 551 * LUT2Y,FNT2Y) 552 553C --------------------------------------------- 554C D(ia) <-- D(ia) + sum_fnm y^t(am,fn) M(imfn): 555C --------------------------------------------- 556 IOPTT2 = 1 557 CALL CCSDT_ETA_TM2(WORK(KDIA2),ISYDEN,WORK(KMMAT),ISYCTR, 558 & DUMMY,LUT2Y,FNT2Y,ISTAMY,IOPTT2, 559 & WORK(KEND0),LWRK0) 560 561C --------------------------------------------- 562C D(ia) <-- D(ia) + sum_fnm t(am,fn) y^M(imfn): 563C --------------------------------------------- 564 IOPTT2 = 0 565 CALL CCSDT_ETA_TM2(WORK(KDIA4),ISYDEN,WORK(KYMMAT),ISYDEN, 566 & WORK(KT2TP),IDUMMY,CDUMMY,ISYM0,IOPTT2, 567 & WORK(KEND0),LWRK0) 568 569C 570C ---------------------------------------------------------------- 571C D(ia) <-- D(ia) - sum_b y^t(bi) D^0(ba) + sum_j y^t(aj) D^0(ij): 572C ---------------------------------------------------------------- 573 ! read D^0(ab) 574 KD0AB = KEND0 575 KD0IJ = KD0AB + NMATAB(ISYM0) 576 KT1AMY = KD0IJ + NMATIJ(ISYM0) 577 KEND2 = KT1AMY + NT1AM(ISTAMY) 578 LWRK2 = LWORK - KEND2 579 580 IF (LWRK2 .LT. 0) THEN 581 CALL QUIT('Out of memory in CCSDT_ETA_DEN.') 582 ENDIF 583 584 IOPT = 1 585 CALL CC_RDRSP(LISTR,IDLSTR,ISTAMY,IOPT,MODEL,WORK(KT1AMY),DUMMY) 586 587 LUPTIJ = -1 588 LUPTAB = -1 589 590 ! read intermediates from ground state density from file... 591 CALL WOPEN2( LUPTIJ,FNDPTIJ,64,0) 592 CALL GETWA2( LUPTIJ,FNDPTIJ,WORK(KD0IJ),1,NMATIJ(ISYM0)) 593 CALL WCLOSE2(LUPTIJ,FNDPTIJ,'KEEP') 594 595 CALL WOPEN2( LUPTAB,FNDPTAB,64,0) 596 CALL GETWA2( LUPTAB,FNDPTAB,WORK(KD0AB),1,NMATAB(ISYM0)) 597 CALL WCLOSE2(LUPTAB,FNDPTAB,'KEEP') 598 599 600 DO ISYMA = 1, NSYM 601 ISYMI = MULD2H(ISYDEN,ISYMA) 602 ISYMB = MULD2H(ISYCTR,ISYMA) 603 ISYMJ = MULD2H(ISYCTR,ISYMI) 604 605 KOFFDBA = KD0AB + IMATAB(ISYMB,ISYMA) 606 KOFFDIJ = KD0IJ + IMATIJ(ISYMI,ISYMJ) 607 KOFFDIA = KDIA1 + IT1AM(ISYMA,ISYMI) 608 KOFFTBI = KT1AMY+ IT1AM(ISYMB,ISYMI) 609 KOFFTAJ = KT1AMY+ IT1AM(ISYMA,ISYMJ) 610 611 NVIRA = MAX(NVIR(ISYMA),1) 612 NVIRB = MAX(NVIR(ISYMB),1) 613 NRHFI = MAX(NRHF(ISYMI),1) 614 615 CALL DGEMM('T','N',NVIR(ISYMA),NRHF(ISYMI),NVIR(ISYMB), 616 & -ONE,WORK(KOFFDBA),NVIRB,WORK(KOFFTBI),NVIRB, 617 & ONE,WORK(KOFFDIA),NVIRA) 618 619 CALL DGEMM('N','T',NVIR(ISYMA),NRHF(ISYMI),NRHF(ISYMJ), 620 & ONE,WORK(KOFFTAJ),NVIRA,WORK(KOFFDIJ),NRHFI, 621 & ONE,WORK(KOFFDIA),NVIRA) 622 623 END DO 624C 625C----------------------------------------- 626C copy/add densities to result arrays: 627C----------------------------------------- 628C 629 CALL DCOPY(NMATIJ(ISYDEN),WORK(KDIJ),1,DENIJ,1) 630 CALL DCOPY(NMATAB(ISYDEN),WORK(KDAB),1,DENAB,1) 631 632 CALL DCOPY(NT1AM(ISYDEN),WORK(KDIA1),1,DENIA,1) 633 CALL DAXPY(NT1AM(ISYDEN),+1.0D0,WORK(KDIA2),1,DENIA,1) 634 CALL DAXPY(NT1AM(ISYDEN),-1.0D0,WORK(KDIA3),1,DENIA,1) 635 CALL DAXPY(NT1AM(ISYDEN),+1.0D0,WORK(KDIA4),1,DENIA,1) 636 637 638C----------------------------------------- 639C close/deleta scratch files: 640C----------------------------------------- 641 CALL WCLOSE2(LUT2Y,FNT2Y,'DELETE') 642C 643 CALL QEXIT('CTETADEN') 644 RETURN 645 END 646*=====================================================================* 647*=====================================================================* 648 SUBROUTINE CCSDT_ETA_FA_DEN(LISTL,IDLSTL,ISYMBL, 649 * LISTR,IDLSTR,ISTAMY,FREQY,LABELY, 650 * LCAUCHY,NCAU, 651 * XLAMDP,XLAMDH,FOCK0, 652 * T2TP,ZL2TP,ZL1AM, 653 * ISYDEN,DAB,DIJ,DO_DIA,DIA, 654 * DO_MMAT,XMMAT,DO_YMMAT,YMMAT, 655 * WORK,LWORK, 656 * LUDELD,FNDELD,LUCKJD, FNCKJD, 657 * LUDKBC,FNDKBC,LUTOC, FNTOC, 658 * LU3VI, FN3VI, LUDKBC3,FNDKBC3, 659 * LU3FOP,FN3FOP,LU3FOP2,FN3FOP2, 660 * LUT2Y,FNT2Y) 661*---------------------------------------------------------------------* 662* 663* Purpose: compute contractions needed for triples contributions 664* to the Eta and F{A} densities 665* 666* 667* Written by Christof Haettig, Mai 2002, Aarhus 668* 669*********************************************************************** 670* 671* Spring 2004, Filip Pawlowski, Aarhus: Modified to treat also the 672* Cauchy moments. 673* 674*=====================================================================* 675 IMPLICIT NONE 676C 677#include "priunit.h" 678#include "dummy.h" 679#include "iratdef.h" 680#include "ccsdsym.h" 681#include "ccorb.h" 682#include "ccsdinp.h" 683#include "ccinftap.h" 684#include "inftap.h" 685#include "ccexci.h" 686C 687 LOGICAL LOCDBG, TOT_T3Y 688 LOGICAL LCAUCHY 689 PARAMETER(LOCDBG = .FALSE., TOT_T3Y = .TRUE.) 690C 691 INTEGER ISYM0 692 PARAMETER(ISYM0 = 1) 693C 694 CHARACTER*10 FNTBAR, FNWMAT 695 PARAMETER(FNTBAR = 'CCTBAR_TMP', FNWMAT = 'CCWMAT_TMP') 696C 697 CHARACTER LISTL*3, LISTR*3 698 CHARACTER*12 FN3SRTR, FNDELDR 699 CHARACTER*16 FNCKJDR, FNDKBCR 700 CHARACTER*(*) FNDELD, FNCKJD, FNDKBC, FNTOC, FN3VI, FNDKBC3 701 CHARACTER*(*) FN3FOP, FN3FOP2, FNT2Y 702C 703 PARAMETER(FN3SRTR = 'CCSDT_ETAFD1', FNDELDR = 'CCSDT_ETAFD3') 704 705 CHARACTER MODEL*10, LABELY*8, CDUMMY*1 706 707 LOGICAL DO_DIA, DO_MMAT, DO_YMMAT 708 709 INTEGER LWORK, IDLSTL, IDLSTR 710 711 INTEGER LU3SRTR, LUCKJDR, LUDELDR, LUDKBCR 712 INTEGER KFOCKD, KEND0, LWRK0, 713 * KTROC0, KXIAJB, KEND1, LWRK1, KINTOC, KEND2, LWRK2, 714 * LENGTH, IOFF, ISYMD, ISYCKB, ISTAMY, 715 * KTRVI1, KTRVI2, KTRVI0, 716 * KTRVI3, KEND3, LWRK3, KINTVI, KEND4, LWRK4, 717 * ISALJB, ISYMBD, ISCKIJ, KBLSMATBD, KSMATBD, 718 * KDIAG, KDIAGW, ISYMC, ISYMK, KOFF1, KOFF2, KOFF3, 719 * KINDSQ, KINDEXB, KTMAT, LENSQ, KINDSQW, 720 * KFCKBA, KTRVI4, KTRVI5, KTRVI6, 721 * KUMATBD, KBLUMATBD, 722 * KTROC02, KTRVI7, KTRVI8, KTRVI9, KTRVI10, 723 * KTRVI11, KTRVI12, KTRVI13, 724 * ISYCKD, KSMATDB, KUMATDB, KEND5, LWRK5, 725 * ISALJD, KINDEXD, KBLSMATDB, KBLUMATDB, 726 * KTRVI14, KTRVI15, KTRVI16, KTRVI17, KTRVI18, KTRVI19, 727 * KTRVI20, ISYTMP, KTROC01, KTROC21, KTROC03, KTROC23, 728 * LUDELD, LUCKJD, LUDKBC, LUTOC, LU3VI, LUDKBC3, LU3FOP, 729 * LU3FOP2, LUTBAR, LUWMAT, LUT2Y, 730 * ISYMAI, ISYML, ISYAIL, ISAIBJ, ISYMDL, ISYMBJ, 731 * KOFFA, IOPT, NBJ, IADR, ISYMJ, ISYMB, ISYMI, NRHFI, 732 * ISYMBL, ISAIB, ISYM, ILEN, ISYMFN, ISYMIM, ISYDEN, 733 * ISWMAT, KFOCKY, IRREP, NVIRB, 734 * IERR, KWMAT, LENSQW, ISTBAR, NDL, ISYMN, ISYEMF, KTBAR, 735 * ISYMBN, ISYMEM, ISYMF, KT2Y, KOFFD, NTOTEM, NNVIRF, 736 * ISYMFI, NNEMF, ISYMM, ISYME, ISYMEI, ISYEIL, 737 * KFN, KOFFM, NVIRE, NRFHI, ISYMDN, KOFFT, KBN, KOFFW, 738 * ISYEMN, ISYMEN, ISYENL, ISYEML, 739 * KT1B, KT2B, KTROC0Y, KTRVI0Y, 740 * IMAT, ISCKBY, KTRVI8Y, 741 * ISCKDY 742 INTEGER ISWTL(8,8), ISTLN(8,8), NIMFN(8) 743 744 INTEGER ISBLALJB,ISBLALJD,ISBLCKIJ,KBLINDEXB,KBLINDEXD 745 INTEGER KBLDIAG,LENSQBL,KBLINDSQ 746C 747 INTEGER NCAU,MCAU,IDLSTRM,KOFFOCC,KOFFINTD,KOFFINTB 748C 749 !functions: 750 INTEGER ILSTSYM,ILRCAMP 751C 752 integer kx3am 753 754#if defined (SYS_CRAY) 755 REAL XLAMDP(*), XLAMDH(*), FOCK0(*) 756 REAL WORK(LWORK) 757 REAL FREQY, DDOT 758 REAL ZERO, ONE, TWO, HALF 759 REAL DAB(*),DIJ(*),DIA(*) 760 REAL XMMAT(*),YMMAT(*) 761 REAL T2TP(*),ZL2TP(*),ZL1AM(*) 762 REAL FREQL 763#else 764 DOUBLE PRECISION XLAMDP(*), XLAMDH(*), FOCK0(*) 765 DOUBLE PRECISION WORK(LWORK) 766 DOUBLE PRECISION FREQY, DDOT 767 DOUBLE PRECISION ZERO, ONE, TWO, HALF 768 DOUBLE PRECISION DAB(*),DIJ(*),DIA(*) 769 DOUBLE PRECISION XMMAT(*),YMMAT(*) 770 DOUBLE PRECISION T2TP(*),ZL2TP(*),ZL1AM(*) 771 DOUBLE PRECISION FREQL 772#endif 773C 774 PARAMETER (ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, TWO = 2.0D0) 775C 776C------------------------------------------------------------ 777C some initializations: 778C set offset arrays ISWTL and ISTLN and dimensions NIMFN: 779C------------------------------------------------------------ 780 CALL QENTER('CTETAFAD') 781C 782 !read frequency for listl 783 IF (LISTL(1:3) .EQ. 'L0 ') THEN 784 FREQL = 0.0D0 785 ELSE IF (LISTL(1:3) .EQ. 'LE ') THEN 786 FREQL = EIGVAL(IDLSTL) 787 ELSE 788 WRITE(LUPRI,*) 'LISTL = ', LISTL(1:3) 789 WRITE(LUPRI,*) 'CCSDT_ETA_FA_DEN handles only L0 or LE as LISTL' 790 CALL QUIT('Unknown left list in CCSDT_ETA_FA_DEN') 791 END IF 792C 793 DO ISYMDL = 1, NSYM 794 IOFF = 0 795 DO ISYML = 1, NSYM 796 ISWMAT = MULD2H(ISYMDL,ISYML) 797 ISWTL(ISWMAT,ISYML) = IOFF 798 IOFF = IOFF + NT2SQ(ISWMAT)*NRHF(ISYML) 799 END DO 800 END DO 801 802 DO ISAIBJ = 1, NSYM 803 IOFF = 0 804 DO ISYMJ = 1, NSYM 805 ISAIB = MULD2H(ISAIBJ,ISYMJ) 806 ISTLN(ISAIB,ISYMJ) = IOFF 807 IOFF = IOFF + NCKATR(ISAIB)*NRHF(ISYMJ) 808 END DO 809 END DO 810 811 DO ISYM = 1, NSYM 812 ILEN = 0 813 DO ISYMFN = 1, NSYM 814 ISYMIM = MULD2H(ISYM,ISYMFN) 815 ILEN = ILEN + NMATIJ(ISYMIM)*NT1AM(ISYMFN) 816 END DO 817 NIMFN(ISYM) = ILEN 818 END DO 819C 820C------------------------------------------------------- 821C initial allocations, prepare T2B amplitudes: 822C------------------------------------------------------- 823C 824 KEND0 = 1 825 826 IF ( (LISTR(1:3).EQ.'R1 ') .OR. (LISTR(1:3).EQ.'RC ') ) THEN 827 KFOCKY = KEND0 828 KEND0 = KFOCKY + N2BST(ISTAMY) 829 END IF 830 831 KFOCKD = KEND0 832 KFCKBA = KFOCKD + NORBTS 833 KEND0 = KFCKBA + NT1AMX 834 835 IF (TOT_T3Y) THEN 836 KT1B = KEND0 837 KT2B = KT1B + (NCAU+1)*NT1AM(ISTAMY) 838 KEND0 = KT2B + (NCAU+1)*NT2SQ(ISTAMY) 839 END IF 840 841 LWRK0 = LWORK - KEND0 842 843 IF (LWRK0.LT.0) THEN 844 CALL QUIT('Not enough memory in CCSDT_ETA_FA_DEN') 845 ENDIF 846 847C 848C --------------------------------------------------------------------- 849C Get KT1B and KT2B vectors 850C 851C For non-Cacuhy calculations (LCAUCHY = .FALSE., NCAU = 0) the 852C loop below is dummy. 853C 854C For Cacuhy calculations (LCAUCHY = .TRUE., NCAU >= 0) 855C KT1B and KT2B correspond to singles and doubles part of Cauchy 856C vectors. 857C --------------------------------------------------------------------- 858C 859 IF (TOT_T3Y) THEN 860C 861 IF (LWRK0.LT.NT2SQ(ISTAMY)) THEN 862 CALL QUIT('Not enough memory for T2Y (CCSDT_ETA_FA_DEN)') 863 ENDIF 864C 865 DO MCAU = 0, NCAU 866C 867 IF (LCAUCHY) THEN 868 !Get the list for current MCAU 869 IDLSTRM = ILRCAMP(LABELY,MCAU,ISTAMY) 870 !Consistency check 871 IF (MCAU.EQ.NCAU .AND. IDLSTRM.NE.IDLSTR) THEN 872 CALL QUIT('Sth wrong in Cauchy loop in CCSDT_ETA_FA_DEN') 873 END IF 874 ELSE 875 IDLSTRM = IDLSTR 876 END IF 877C 878 IOPT = 3 879 KOFF1 = KT1B + MCAU*NT1AM(ISTAMY) 880 KOFF2 = KT2B + MCAU*NT2SQ(ISTAMY) 881 CALL CC_RDRSP(LISTR,IDLSTRM,ISTAMY,IOPT,MODEL, 882 * WORK(KOFF1),WORK(KOFF2)) 883 CALL CCLR_DIASCL(WORK(KOFF2),TWO,ISTAMY) 884 CALL CC_T2SQ(WORK(KOFF2),WORK(KEND0),ISTAMY) 885 CALL CC3_T2TP(WORK(KOFF2),WORK(KEND0),ISTAMY) 886C 887 END DO !MCAU 888C 889 END IF 890C 891C ---------------------------------------------------------------------- 892C Calculate (ck|de)-tilde(Y) and (ck|lm)-tilde(Y) needed in construction 893C of W3 ( <mu3|[[H^0,T1B],T2^0]|HF> term). 894C ---------------------------------------------------------------------- 895C 896 IF (TOT_T3Y) THEN 897C 898 DO MCAU = 0, NCAU 899C 900 !Open temporary files 901 LU3SRTR = -1 902 LUDELDR = -1 903 CALL WOPEN2(LU3SRTR,FN3SRTR,64,0) 904 CALL WOPEN2(LUDELDR,FNDELDR,64,0) 905C 906 !Generate file names for T1-transformed integrals 907 IF (LCAUCHY) THEN 908 CALL GET_FILENAME(FNCKJDR,FNDKBCR,MCAU) 909 ELSE 910 FNCKJDR = 'CCSDT_____ETAFD2' 911 FNDKBCR = 'CCSDT_____ETAFD4' 912 END IF 913C 914 !Open the files for T1-transformed integrals 915 LUCKJDR = -1 916 LUDKBCR = -1 917 CALL WOPEN2(LUCKJDR,FNCKJDR,64,0) 918 CALL WOPEN2(LUDKBCR,FNDKBCR,64,0) 919C 920 !Get offset for the T1Y amplitudes 921 KOFF1 = KT1B + MCAU*NT1AM(ISTAMY) 922C 923 !Construct T1-transformed integrals and store on file 924 CALL CC3_BARINT(WORK(KOFF1),ISTAMY,XLAMDP, 925 * XLAMDH,WORK(KEND0),LWRK0, 926 * LU3SRTR,FN3SRTR,LUCKJDR,FNCKJDR) 927C 928 CALL CC3_SORT1(WORK(KEND0),LWRK0,2,ISTAMY,LU3SRTR,FN3SRTR, 929 * LUDELDR,FNDELDR,IDUMMY,CDUMMY,IDUMMY,CDUMMY, 930 * IDUMMY,CDUMMY) 931C 932 CALL CC3_SINT(XLAMDH,WORK(KEND0),LWRK0,ISTAMY, 933 * LUDELDR,FNDELDR,LUDKBCR,FNDKBCR) 934C 935 !Close the files for T1-transformed integrals 936 CALL WCLOSE2(LUCKJDR,FNCKJDR,'KEEP') 937 CALL WCLOSE2(LUDKBCR,FNDKBCR,'KEEP') 938C 939 !Close and delete temporary files 940 CALL WCLOSE2(LU3SRTR,FN3SRTR,'DELETE') 941 CALL WCLOSE2(LUDELDR,FNDELDR,'DELETE') 942C 943 END DO !MCAU 944C 945 ENDIF 946c 947 if (.false.) then 948 kx3am = kend0 949 kend0 = kx3am + nt1amx*nt1amx*nt1amx 950 call dzero(work(kx3am),nt1amx*nt1amx*nt1amx) 951 lwrk0 = lwork - kend0 952 if (lwrk0 .lt. 0) then 953 write(lupri,*) 'Memory available : ',lwork 954 write(lupri,*) 'Memory needed : ',kend0 955 call quit('Insufficient space (T3) in CCSDT_ETA_FA_DEN') 956 END IF 957 endif 958C 959C------------------------------------- 960C Read canonical orbital energies: 961C------------------------------------- 962C 963 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY, 964 & .FALSE.) 965 REWIND LUSIFC 966 967 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 968 READ (LUSIFC) 969 READ (LUSIFC) (WORK(KFOCKD+I-1), I=1,NORBTS) 970 971 CALL GPCLOSE(LUSIFC,'KEEP') 972C 973C--------------------------------------------- 974C Delete frozen orbitals in Fock diagonal. 975C--------------------------------------------- 976C 977 IF (FROIMP .OR. FROEXP) 978 * CALL CCSD_DELFRO(WORK(KFOCKD),WORK(KEND0),LWRK0) 979C 980C----------------------------------------------------- 981C Sort the fock matrix 982C----------------------------------------------------- 983C 984 DO ISYMC = 1,NSYM 985 986 ISYMK = MULD2H(ISYMC,ISYM0) 987 988 DO K = 1,NRHF(ISYMK) 989 990 DO C = 1,NVIR(ISYMC) 991 992 KOFF1 = IFCVIR(ISYMK,ISYMC) + 993 * NORB(ISYMK)*(C - 1) + K 994 KOFF2 = IT1AM(ISYMC,ISYMK) 995 * + NVIR(ISYMC)*(K - 1) + C 996 997 WORK(KFCKBA-1+KOFF2) = FOCK0(KOFF1) 998 999 ENDDO 1000 ENDDO 1001 ENDDO 1002 1003 IF (LOCDBG) THEN 1004 CALL AROUND('In CCSDT_ETA_FA_DEN: Triples Fock MO matrix (sort)') 1005 CALL CC_PRFCKMO(WORK(KFCKBA),ISYM0) 1006 ENDIF 1007C 1008C----------------------------------------------------------- 1009C Prepare one-electron operators needed to compute the 1010C amplitude response vectors: 1011C----------------------------------------------------------- 1012C 1013 IF ( (LISTR(1:3).EQ.'R1 ') .OR. (LISTR(1:3).EQ.'RC ') ) THEN 1014 ! read property integrals from file: 1015 CALL CCPRPAO(LABELY,.TRUE.,WORK(KFOCKY),IRREP,IMAT,IERR, 1016 * WORK(KEND0),LWRK0) 1017 IF ((IERR.GT.0) .OR. (IERR.EQ.0 .AND. IRREP.NE.ISTAMY)) THEN 1018 WRITE(LUPRI,*) 'ISTAMY:',ISTAMY 1019 WRITE(LUPRI,*) 'IRREP :',IRREP 1020 WRITE(LUPRI,*) 'IERR :',IERR 1021 CALL QUIT('CCSDT_ETA_FA_DEN: error reading operator '//LABELY) 1022 ELSE IF (IERR.LT.0) THEN 1023 CALL DZERO(WORK(KFOCKY),N2BST(ISTAMY)) 1024 END IF 1025 1026 ! transform property integrals to Lambda-MO basis 1027 CALL CC_FCKMO(WORK(KFOCKY),XLAMDP,XLAMDH,WORK(KEND0),LWRK0, 1028 * ISTAMY,1,1) 1029 END IF 1030 1031C 1032C----------------------------- 1033C Memory allocation. 1034C----------------------------- 1035C 1036 KTROC0 = KEND0 1037 KTROC02 = KTROC0 + NTRAOC(ISYM0) 1038 KXIAJB = KTROC02 + NTRAOC(ISYM0) 1039 KEND1 = KXIAJB + NT2AM(ISYM0) 1040 IF (TOT_T3Y) THEN 1041 KTROC0Y = KEND1 1042 KEND1 = KTROC0Y + (NCAU+1)*NTRAOC(ISTAMY) 1043 ENDIF 1044 1045 KTROC01 = KEND1 1046 KTROC21 = KTROC01 + NTRAOC(ISYM0) 1047 KTROC03 = KTROC21 + NTRAOC(ISYM0) 1048 KTROC23 = KTROC03 + NTRAOC(ISYM0) 1049 KEND1 = KTROC23 + NTRAOC(ISYM0) 1050 LWRK1 = LWORK - KEND1 1051 1052 KINTOC = KEND1 1053 KEND2 = KINTOC + MAX(NTOTOC(ISYM0),NTOTOC(ISTAMY)) 1054 LWRK2 = LWORK - KEND2 1055 1056 IF (LWRK2 .LT. 0) THEN 1057 WRITE(LUPRI,*) 'Memory available : ',LWORK 1058 WRITE(LUPRI,*) 'Memory needed : ',KEND2 1059 CALL QUIT('Insufficient space in CCSDT_ETA_FA_DEN') 1060 END IF 1061C 1062C------------------------ 1063C Construct L(ia,jb). 1064C------------------------ 1065C 1066 LENGTH = IRAT*NT2AM(ISYM0) 1067 1068 REWIND(LUIAJB) 1069 CALL READI(LUIAJB,LENGTH,WORK(KXIAJB)) 1070 1071 CALL CCSD_TCMEPK(WORK(KXIAJB),1.0D0,ISYM0,1) 1072C 1073C------------------------ 1074C Occupied integrals. 1075C------------------------ 1076C 1077 IOFF = 1 1078 IF (NTOTOC(ISYM0) .GT. 0) THEN 1079 CALL GETWA2(LUCKJD,FNCKJD,WORK(KINTOC),IOFF,NTOTOC(ISYM0)) 1080 ENDIF 1081 1082 IF (LOCDBG) WRITE(LUPRI,*) 'Norm of OCC-INT ', 1083 * DDOT(NTOTOC(ISYM0),WORK(KINTOC),1,WORK(KINTOC),1) 1084C 1085C---------------------------------------------------------------------- 1086C Transform (ia|j delta) integrals to (ia|j k) and sort as (ij,k,a) 1087C---------------------------------------------------------------------- 1088C 1089 CALL CC3_TROCC(WORK(KINTOC),WORK(KTROC0),XLAMDP, 1090 * WORK(KEND2),LWRK2,ISYM0) 1091 1092C 1093C-------------------------------------------------------------------- 1094C Read in T1Y-transformed occupied integrals used to construct W3 1095C-------------------------------------------------------------------- 1096C 1097 IF (TOT_T3Y) THEN 1098C 1099 DO MCAU = 0,NCAU 1100C 1101 !Generate file name for T1Y-transformed integrals 1102 IF (LCAUCHY) THEN 1103 !(FNDKBCR is not needed here) 1104 CALL GET_FILENAME(FNCKJDR,FNDKBCR,MCAU) 1105 ELSE 1106 FNCKJDR = 'CCSDT_____ETAFD2' 1107 END IF 1108C 1109 !Open the file for T1Y-transformed integrals 1110 LUCKJDR = -1 1111 CALL WOPEN2(LUCKJDR,FNCKJDR,64,0) 1112C 1113 !Get the offset for the integrals depending on MCAU 1114 KOFF1 = KTROC0Y + MCAU*NTRAOC(ISTAMY) 1115C 1116 IOFF = 1 1117 IF (NTOTOC(ISTAMY) .GT. 0) THEN 1118 CALL GETWA2(LUCKJDR,FNCKJDR,WORK(KINTOC),IOFF, 1119 * NTOTOC(ISTAMY)) 1120 ENDIF 1121C 1122 IF (LOCDBG) 1123 * WRITE(LUPRI,*) 'Norm of CCSDT_OC-INT (Y transformed) ', 1124 * DDOT(NTOTOC(ISTAMY),WORK(KINTOC),1,WORK(KINTOC),1) 1125C 1126 CALL CC3_TROCC(WORK(KINTOC),WORK(KOFF1),XLAMDP, 1127 * WORK(KEND2),LWRK2,ISTAMY) 1128C 1129 !Close and delete the file for T1Y-transformed occupied integrals 1130 CALL WCLOSE2(LUCKJDR,FNCKJDR,'DELETE') 1131C 1132 END DO !MCAU 1133C 1134 ENDIF 1135C 1136C----------------------------------------------------------- 1137C Construct 2*C-E of the integrals. 1138C Have integral for both (ij,k,a) and (a,k,j,i) 1139C----------------------------------------------------------- 1140C 1141 IOFF = 1 1142 IF (NTOTOC(ISYM0) .GT. 0) THEN 1143 CALL GETWA2(LUTOC,FNTOC,WORK(KINTOC),IOFF,NTOTOC(ISYM0)) 1144 ENDIF 1145 1146 CALL CC3_TROCC(WORK(KINTOC),WORK(KTROC01),XLAMDH, 1147 * WORK(KEND2),LWRK2,ISYM0) 1148 IF (LOCDBG) WRITE(LUPRI,*) 'Norm of CKJDEL-INT ', 1149 * DDOT(NTOTOC(ISYM0),WORK(KINTOC),1,WORK(KINTOC),1) 1150 1151 CALL CCSDT_TCMEOCC(WORK(KTROC01),WORK(KTROC21),ISYM0) 1152 1153 CALL CCFOP_SORT(WORK(KTROC01),WORK(KTROC03),ISYM0,1) 1154 1155 CALL CCFOP_SORT(WORK(KTROC21),WORK(KTROC23),ISYM0,1) 1156 1157 CALL CCFOP_SORT(WORK(KTROC0),WORK(KTROC02),ISYM0,1) 1158 IF (LOCDBG) WRITE(LUPRI,*) 'Norm of TROC0 ', 1159 * DDOT(NTRAOC(ISYM0),WORK(KTROC0),1,WORK(KTROC0),1) 1160C 1161C--------------------------------------------- 1162C Open files for Tbar and W intermediates: 1163C--------------------------------------------- 1164C 1165 LUTBAR = -1 1166 LUWMAT = -1 1167 1168 CALL WOPEN2(LUTBAR,FNTBAR,64,0) 1169 CALL WOPEN2(LUWMAT,FNWMAT,64,0) 1170 1171C 1172C---------------------------- 1173C Loop over D 1174C---------------------------- 1175C 1176 DO ISYMD = 1,NSYM 1177 1178 ISYCKB = MULD2H(ISYMD,ISYM0) 1179 ISCKBY = MULD2H(ISTAMY,ISYMD) 1180 IF (LOCDBG) 1181 * WRITE(LUPRI,*) 'In CCSDT_ETA_FA_DEN: ISYCKB :',ISYCKB 1182 1183 DO D = 1,NVIR(ISYMD) 1184C 1185C ------------------ 1186C Memory allocation. 1187C ------------------ 1188 KTRVI0 = KEND1 1189 KTRVI1 = KTRVI0 + NCKATR(ISYCKB) 1190 KTRVI2 = KTRVI1 + NCKATR(ISYCKB) 1191 KTRVI3 = KTRVI2 + NCKATR(ISYCKB) 1192 KTRVI4 = KTRVI3 + NCKATR(ISYCKB) 1193 KTRVI5 = KTRVI4 + NCKATR(ISYCKB) 1194 KTRVI6 = KTRVI5 + NCKATR(ISYCKB) 1195 KTRVI7 = KTRVI6 + NCKATR(ISYCKB) 1196 KEND3 = KTRVI7 + NCKATR(ISYCKB) 1197 LWRK3 = LWORK - KEND3 1198 1199 KTRVI14 = KEND3 1200 KTRVI15 = KTRVI14 + NCKATR(ISYCKB) 1201 KTRVI18 = KTRVI15 + NCKATR(ISYCKB) 1202 KTRVI19 = KTRVI18 + NCKATR(ISYCKB) 1203 KEND3 = KTRVI19 + NCKATR(ISYCKB) 1204 LWRK3 = LWORK - KEND3 1205 1206 IF (TOT_T3Y) THEN 1207 KTRVI0Y = KEND3 1208 KEND3 = KTRVI0Y + (NCAU+1)*NCKATR(ISCKBY) 1209 LWRK3 = LWORK - KEND3 1210 ENDIF 1211 1212 KINTVI = KEND3 1213 KEND4 = KINTVI + MAX(NCKA(ISYMD),NCKA(ISYCKB)) 1214 LWRK4 = LWORK - KEND4 1215 1216 IF (LWRK4 .LT. 0) THEN 1217 WRITE(LUPRI,*) 'Memory available : ',LWORK 1218 WRITE(LUPRI,*) 'Memory needed : ',KEND4 1219 CALL QUIT('Insufficient space in CCSDT_ETA_FA_DEN(99)') 1220 END IF 1221C 1222C ------------------------------------ 1223C Integrals used in s3am. 1224C ------------------------------------ 1225C 1226 IOFF = ICKBD(ISYCKB,ISYMD) + NCKATR(ISYCKB)*(D - 1) + 1 1227 IF (NCKATR(ISYCKB) .GT. 0) THEN 1228 CALL GETWA2(LUDKBC,FNDKBC,WORK(KTRVI0),IOFF, 1229 & NCKATR(ISYCKB)) 1230 ENDIF 1231 1232C 1233C ------------------------------------------------------------- 1234C Read T1Y-transformed virt ints (fixed D) used to construct W3 1235C ------------------------------------------------------------- 1236C 1237 IF (TOT_T3Y) THEN 1238C 1239 DO MCAU = 0, NCAU 1240C 1241 !Generate file names for T1Y-transformed integrals 1242 IF (LCAUCHY) THEN 1243 !(FNCKJDR is not needed here) 1244 CALL GET_FILENAME(FNCKJDR,FNDKBCR,MCAU) 1245 ELSE 1246 FNDKBCR = 'CCSDT_____ETAFD4' 1247 END IF 1248C 1249 !Open the files for T1Y-transformed integrals 1250 LUDKBCR = -1 1251 CALL WOPEN2(LUDKBCR,FNDKBCR,64,0) 1252C 1253 !Get the offset for a given MCAU 1254 KOFF1 = KTRVI0Y + MCAU*NCKATR(ISCKBY) 1255C 1256 !Read in the integrals 1257 IOFF = ICKBD(ISCKBY,ISYMD) + NCKATR(ISCKBY)*(D-1) + 1 1258 IF (NCKATR(ISCKBY) .GT. 0) THEN 1259 CALL GETWA2(LUDKBCR,FNDKBCR,WORK(KOFF1),IOFF, 1260 & NCKATR(ISCKBY)) 1261 ENDIF 1262C 1263 !Close the file for T1B-transformed virtual integrals 1264 !(and keep it: it will be needed in B-loop) 1265 CALL WCLOSE2(LUDKBCR,FNDKBCR,'KEEP') 1266C 1267 END DO !MCAU 1268C 1269 ENDIF 1270C 1271C ------------------------------------------- 1272C Read 2*C-E of integral used for t3-bar 1273C ------------------------------------------- 1274C 1275 IOFF = ICKBD(ISYCKB,ISYMD) + NCKATR(ISYCKB)*(D - 1) + 1 1276 IF (NCKATR(ISYCKB) .GT. 0) THEN 1277 CALL GETWA2(LU3FOP2,FN3FOP2,WORK(KTRVI4),IOFF, 1278 & NCKATR(ISYCKB)) 1279 ENDIF 1280C 1281C ------------------------------------------------- 1282C Integrals used for t3-bar for cc3 1283C ------------------------------------------------- 1284C 1285 IOFF = ICKBD(ISYCKB,ISYMD) + NCKATR(ISYCKB)*(D - 1) + 1 1286 IF (NCKATR(ISYCKB) .GT. 0) THEN 1287 CALL GETWA2(LUDKBC3,FNDKBC3,WORK(KTRVI14),IOFF, 1288 & NCKATR(ISYCKB)) 1289 ENDIF 1290 CALL CCSDT_SRVIR3(WORK(KTRVI14),WORK(KEND4), 1291 * ISYMD,D,ISYM0) 1292 CALL CCSDT_SRTVIR(WORK(KTRVI14),WORK(KTRVI15),WORK(KEND4) 1293 * ,LWRK4,ISYMD,ISYM0) 1294C 1295C ------------------------------------------------ 1296C Sort the integrals for s3am and for t3-bar 1297C ------------------------------------------------ 1298C 1299 CALL CCSDT_SRTVIR(WORK(KTRVI0),WORK(KTRVI2),WORK(KEND4), 1300 * LWRK4,ISYMD,ISYM0) 1301 1302 CALL CCSDT_SRTVIR(WORK(KTRVI4),WORK(KTRVI5),WORK(KEND4), 1303 * LWRK4,ISYMD,ISYM0) 1304C 1305C ------------------------------------------- 1306C Read virtual integrals used in contraction. 1307C ------------------------------------------- 1308C 1309 IOFF = ICKAD(ISYCKB,ISYMD) + NCKA(ISYCKB)*(D - 1) + 1 1310 IF (NCKA(ISYCKB) .GT. 0) THEN 1311 CALL GETWA2(LUDELD,FNDELD,WORK(KINTVI),IOFF, 1312 * NCKA(ISYCKB)) 1313 ENDIF 1314 1315 CALL CCSDT_TRVIR(WORK(KINTVI),WORK(KTRVI1),XLAMDH, 1316 * ISYMD,D,ISYM0,WORK(KEND4),LWRK4) 1317C 1318C --------------------------------------------- 1319C Calculate virtual integrals used in q3am. 1320C --------------------------------------------- 1321C 1322 CALL DCOPY(NCKATR(ISYCKB),WORK(KTRVI1),1,WORK(KTRVI3),1) 1323 1324 IF (LWRK3 .LT. NCKATR(ISYCKB)) THEN 1325 CALL QUIT('Insufficient space for allocation in '// 1326 & 'CCSDT_ETA_FA_DEN (1)') 1327 END IF 1328 1329 CALL CCSDT_SRVIR3(WORK(KTRVI3),WORK(KEND4),ISYMD,D,ISYM0) 1330C 1331C ---------------------------------------------------- 1332C Read virtual integrals used in q3am/u3am for t3-bar. 1333C ---------------------------------------------------- 1334C 1335 IOFF = ICKAD(ISYCKB,ISYMD) + NCKA(ISYCKB)*(D - 1) + 1 1336 IF (NCKA(ISYCKB) .GT. 0) THEN 1337 CALL GETWA2(LU3VI,FN3VI,WORK(KINTVI),IOFF, 1338 * NCKA(ISYCKB)) 1339 ENDIF 1340 1341 CALL CCSDT_TRVIR(WORK(KINTVI),WORK(KTRVI19),XLAMDP, 1342 * ISYMD,D,ISYM0,WORK(KEND4),LWRK4) 1343 1344 IF (LWRK4 .LT. NCKATR(ISYCKB)) THEN 1345 CALL QUIT('Insufficient space for allocation in '// 1346 * 'CCSDT_ETA_FA_DEN (CC3 TRVI)') 1347 END IF 1348 1349 CALL CCSDT_SRTVIR(WORK(KTRVI19),WORK(KTRVI18),WORK(KEND4) 1350 * ,LWRK4,ISYMD,ISYM0) 1351 1352 IOFF = ICKBD(ISYCKB,ISYMD) + NCKATR(ISYCKB)*(D - 1) + 1 1353 IF (NCKATR(ISYCKB) .GT. 0) THEN 1354 CALL GETWA2(LU3FOP,FN3FOP,WORK(KTRVI6),IOFF, 1355 * NCKATR(ISYCKB)) 1356 ENDIF 1357 1358 IF (LWRK3 .LT. NCKATR(ISYCKB)) THEN 1359 CALL QUIT('Insufficient space for allocation in '// 1360 & 'CCSDT_ETA_FA_DEN (2)') 1361 END IF 1362 1363 CALL DCOPY(NCKATR(ISYCKB),WORK(KTRVI6),1,WORK(KTRVI7),1) 1364 1365 CALL CCSDT_SRVIR3(WORK(KTRVI6),WORK(KEND4),ISYMD,D,ISYM0) 1366C 1367C 1368 DO ISYMB = 1,NSYM 1369 1370 ISYMBD = MULD2H(ISYMB,ISYMD) 1371 1372 ISALJB = MULD2H(ISYMB,ISYM0) 1373 ISALJD = MULD2H(ISYMD,ISYM0) 1374 ISBLALJB = MULD2H(ISYMB,ISYMBL) 1375 ISBLALJD = MULD2H(ISYMD,ISYMBL) 1376C 1377 ISCKIJ = MULD2H(ISYMBD,ISYM0) 1378 ISBLCKIJ = MULD2H(ISYMBD,ISYMBL) 1379C 1380 ISYCKD = MULD2H(ISYM0,ISYMB) 1381 ISCKDY = MULD2H(ISTAMY,ISYMB) 1382 ISWMAT = MULD2H(ISCKIJ,ISTAMY) 1383 1384 IF (LOCDBG) THEN 1385 WRITE(LUPRI,*) 'In CCSDT_ETA_FA_DEN: ISYMD :',ISYMD 1386 WRITE(LUPRI,*) 'In CCSDT_ETA_FA_DEN: ISYMB :',ISYMB 1387 WRITE(LUPRI,*) 'In CCSDT_ETA_FA_DEN: ISALJB:',ISALJB 1388 WRITE(LUPRI,*) 'In CCSDT_ETA_FA_DEN: ISYMBD:',ISYMBD 1389 WRITE(LUPRI,*) 'In CCSDT_ETA_FA_DEN: ISCKIJ:',ISCKIJ 1390 ENDIF 1391C 1392C Can use kend3 since we do not need the integrals anymore. 1393 KSMATBD = KEND3 1394 KBLSMATBD = KSMATBD + NCKIJ(ISCKIJ) 1395 KSMATDB = KBLSMATBD + NCKIJ(ISBLCKIJ) 1396 KUMATBD = KSMATDB + NCKIJ(ISCKIJ) 1397 KBLUMATBD = KUMATBD + NCKIJ(ISCKIJ) 1398 KUMATDB = KBLUMATBD + NCKIJ(ISBLCKIJ) 1399 KDIAG = KUMATDB + NCKIJ(ISCKIJ) 1400 KDIAGW = KDIAG + NCKIJ(ISCKIJ) 1401 KINDSQW = KDIAGW + NCKIJ(ISWMAT) 1402 KINDSQ = KINDSQW + (6*NCKIJ(ISWMAT) - 1)/IRAT + 1 1403 KINDEXB = KINDSQ + (6*NCKIJ(ISCKIJ) - 1)/IRAT + 1 1404 KINDEXD = KINDEXB + (NCKI(ISALJB) - 1)/IRAT + 1 1405 KTMAT = KINDEXD + (NCKI(ISALJD) - 1)/IRAT + 1 1406 KTRVI8 = KTMAT + NCKIJMAX 1407 KTRVI9 = KTRVI8 + NCKATR(ISYCKD) 1408 KTRVI10 = KTRVI9 + NCKATR(ISYCKD) 1409 KEND4 = KTRVI10 + NCKATR(ISYCKD) 1410 LWRK4 = LWORK - KEND4 1411C 1412 KBLINDEXB = KEND4 1413 KBLINDEXD = KBLINDEXB + (NCKI(ISBLALJB) - 1)/IRAT + 1 1414 KEND4 = KBLINDEXD + (NCKI(ISBLALJD) - 1)/IRAT + 1 1415 LWRK4 = LWORK - KEND4 1416C 1417 KBLINDSQ = KEND4 1418 KEND4 = KBLINDSQ + (6*NCKIJ(ISBLCKIJ) - 1)/IRAT + 1 1419 LWRK4 = LWORK - KEND4 1420C 1421 KBLDIAG = KEND4 1422 KEND4 = KBLDIAG + NCKIJ(ISBLCKIJ) 1423 LWRK4 = LWORK - KEND4 1424C 1425 KTRVI16 = KEND4 1426 KTRVI17 = KTRVI16 + NCKATR(ISYCKD) 1427 KTRVI20 = KTRVI17 + NCKATR(ISYCKD) 1428 KEND4 = KTRVI20 + NCKATR(ISYCKD) 1429 LWRK4 = LWORK - KEND4 1430 IF (TOT_T3Y) THEN 1431 KTRVI8Y = KEND4 1432 KEND4 = KTRVI8Y + (NCAU+1)*NCKATR(ISCKDY) 1433 ENDIF 1434 1435 KBLSMATDB = KEND4 1436 KBLUMATDB = KBLSMATDB + NCKIJ(ISBLCKIJ) 1437 KTRVI11 = KBLUMATDB + NCKIJ(ISBLCKIJ) 1438 KTRVI12 = KTRVI11 + NCKATR(ISYCKD) 1439 KTRVI13 = KTRVI12 + NCKATR(ISYCKD) 1440 KEND4 = KTRVI13 + NCKATR(ISYCKD) 1441 LWRK4 = LWORK - KEND4 1442 1443 KWMAT = KEND4 1444 KEND4 = KWMAT + NCKIJ(ISWMAT) 1445 LWRK4 = LWORK - KEND4 1446 1447 KINTVI = KEND4 1448 KEND5 = KINTVI + MAX(NCKA(ISYMB),NCKA(ISYCKD)) 1449 LWRK5 = LWORK - KEND5 1450 1451 IF (LWRK5 .LT. 0) THEN 1452 WRITE(LUPRI,*) 'Memory available : ',LWORK 1453 WRITE(LUPRI,*) 'Memory needed : ',KEND5 1454 CALL QUIT('Insufficient space in CCSDT_ETA_FA_DEN') 1455 END IF 1456C 1457C ------------------------------- 1458C Construct part of the diagonal. 1459C ------------------------------- 1460C 1461 CALL CC3_DIAG(WORK(KDIAG),WORK(KFOCKD),ISCKIJ) 1462 CALL CC3_DIAG(WORK(KBLDIAG),WORK(KFOCKD),ISBLCKIJ) 1463C 1464 CALL CC3_DIAG(WORK(KDIAGW),WORK(KFOCKD),ISWMAT) 1465 1466 IF (LOCDBG) THEN 1467 WRITE(LUPRI,*) 'Norm of DIA ', 1468 * DDOT(NCKIJ(ISCKIJ),WORK(KDIAG),1,WORK(KDIAG),1) 1469 WRITE(LUPRI,*) 'Norm of DIA_W', 1470 * DDOT(NCKIJ(ISWMAT),WORK(KDIAGW),1,WORK(KDIAGW),1) 1471 END IF 1472C 1473C ----------------------- 1474C Construct index arrays. 1475C ----------------------- 1476C 1477 LENSQ = NCKIJ(ISCKIJ) 1478 CALL CC3_INDSQ(WORK(KINDSQ),LENSQ,ISCKIJ) 1479 LENSQBL = NCKIJ(ISBLCKIJ) 1480 CALL CC3_INDSQ(WORK(KBLINDSQ),LENSQBL,ISBLCKIJ) 1481C 1482 LENSQW = NCKIJ(ISWMAT) 1483 CALL CC3_INDSQ(WORK(KINDSQW),LENSQW,ISWMAT) 1484C 1485 CALL CC3_INDEX(WORK(KBLINDEXB),ISBLALJB) 1486 CALL CC3_INDEX(WORK(KBLINDEXD),ISBLALJD) 1487 CALL CC3_INDEX(WORK(KINDEXB),ISALJB) 1488 CALL CC3_INDEX(WORK(KINDEXD),ISALJD) 1489 1490 DO B = 1,NVIR(ISYMB) 1491 IF (LOCDBG) write(lupri,*) 'For B, D : ',B,D 1492C 1493C ---------------------------- 1494C Read and transform integrals 1495C ---------------------------- 1496 IOFF = ICKBD(ISYCKD,ISYMB) 1497 * + NCKATR(ISYCKD)*(B - 1) + 1 1498 IF (NCKATR(ISYCKD) .GT. 0) THEN 1499 CALL GETWA2(LUDKBC3,FNDKBC3,WORK(KTRVI16),IOFF, 1500 * NCKATR(ISYCKD)) 1501 ENDIF 1502 CALL CCSDT_SRVIR3(WORK(KTRVI16),WORK(KEND5), 1503 * ISYMB,B,ISYM0) 1504 CALL CCSDT_SRTVIR(WORK(KTRVI16),WORK(KTRVI17), 1505 * WORK(KEND4),LWRK4,ISYMB,ISYM0) 1506C 1507C ------------------------------------ 1508C Read and transform integrals used in 1509C second S-bar and U-bar 1510C ------------------------------------ 1511 IOFF = ICKBD(ISYCKD,ISYMB) 1512 * + NCKATR(ISYCKD)*(B-1) + 1 1513 IF (NCKATR(ISYCKD) .GT. 0) THEN 1514 CALL GETWA2(LU3FOP2,FN3FOP2,WORK(KTRVI11), 1515 * IOFF,NCKATR(ISYCKD)) 1516 ENDIF 1517 1518 CALL CCSDT_SRTVIR(WORK(KTRVI11),WORK(KTRVI12), 1519 * WORK(KEND5),LWRK5,ISYMB, 1520 * ISYM0) 1521 1522 IOFF = ICKBD(ISYCKD,ISYMB) 1523 * + NCKATR(ISYCKD)*(B - 1) + 1 1524 IF (NCKATR(ISYCKD) .GT. 0) THEN 1525 CALL GETWA2(LU3FOP,FN3FOP,WORK(KTRVI13),IOFF, 1526 * NCKATR(ISYCKD)) 1527 ENDIF 1528 1529 IOFF = ICKAD(ISYCKD,ISYMB) 1530 * + NCKA(ISYCKD)*(B - 1) + 1 1531 IF (NCKA(ISYCKD) .GT. 0) THEN 1532 CALL GETWA2(LU3VI,FN3VI,WORK(KINTVI),IOFF, 1533 * NCKA(ISYCKD)) 1534 ENDIF 1535 1536 CALL CCSDT_TRVIR(WORK(KINTVI),WORK(KTRVI20), 1537 * XLAMDP,ISYMB,B,ISYM0, 1538 * WORK(KEND4),LWRK4) 1539C 1540C ---------------------------------------------------- 1541C Calculate the S(ci,bk,dj) matrix for B,D for T3-BAR: 1542C ---------------------------------------------------- 1543 CALL DZERO(WORK(KBLSMATBD),NCKIJ(ISBLCKIJ)) 1544 1545 CALL CCFOP_SMAT(FREQL,ZL1AM,ISYMBL,ZL2TP, 1546 * ISYMBL,WORK(KTMAT), 1547 * WORK(KFCKBA),WORK(KXIAJB),ISYM0, 1548 * WORK(KTRVI14),WORK(KTRVI15), 1549 * WORK(KTRVI4),WORK(KTRVI5), 1550 * WORK(KTROC01),WORK(KTROC21), 1551 * ISYM0,WORK(KFOCKD),WORK(KBLDIAG), 1552 * WORK(KBLSMATBD),WORK(KEND4),LWRK4, 1553 * WORK(KBLINDEXB), 1554 * WORK(KBLINDSQ),LENSQBL, 1555 * ISYMB,B,ISYMD,D) 1556 1557 CALL DSCAL(NCKIJ(ISBLCKIJ),HALF,WORK(KBLSMATBD),1) 1558 1559 CALL T3_FORBIDDEN(WORK(KBLSMATBD),ISYMBL,ISYMB,B, 1560 * ISYMD,D) 1561C 1562C ---------------------------------------------------- 1563C Calculate the S(ci,bk,dj) matrix for D,B for T3-BAR: 1564C ---------------------------------------------------- 1565 CALL DZERO(WORK(KBLSMATDB),NCKIJ(ISBLCKIJ)) 1566 1567 CALL CCFOP_SMAT(FREQL,ZL1AM,ISYMBL,ZL2TP, 1568 * ISYMBL,WORK(KTMAT),WORK(KFCKBA), 1569 * WORK(KXIAJB),ISYM0, 1570 * WORK(KTRVI16),WORK(KTRVI17), 1571 * WORK(KTRVI11),WORK(KTRVI12), 1572 * WORK(KTROC01),WORK(KTROC21), 1573 * ISYM0,WORK(KFOCKD),WORK(KBLDIAG), 1574 * WORK(KBLSMATDB),WORK(KEND4),LWRK4, 1575 * WORK(KBLINDEXD),WORK(KBLINDSQ), 1576 * LENSQBL,ISYMD,D,ISYMB,B) 1577 1578 CALL DSCAL(NCKIJ(ISBLCKIJ),HALF,WORK(KBLSMATDB),1) 1579C 1580C 1581 CALL T3_FORBIDDEN(WORK(KBLSMATDB),ISYMBL,ISYMD,D, 1582 * ISYMB,B) 1583C 1584C ------------------------------------------------ 1585C Calculate U(ci,jk) for fixed b,d for t3-bar. 1586C ------------------------------------------------ 1587 CALL DZERO(WORK(KBLUMATBD),NCKIJ(ISBLCKIJ)) 1588 1589 CALL CCFOP_UMAT(FREQL,ZL1AM,ISYMBL,ZL2TP, 1590 * ISYMBL, 1591 * WORK(KXIAJB),ISYM0,WORK(KFCKBA), 1592 * WORK(KTRVI19),WORK(KTRVI7), 1593 * WORK(KTROC03),WORK(KTROC23),ISYM0, 1594 * WORK(KFOCKD),WORK(KBLDIAG), 1595 * WORK(KBLUMATBD), 1596 * WORK(KTMAT),WORK(KEND4),LWRK4, 1597 * WORK(KBLINDSQ),LENSQBL, 1598 * ISYMB,B,ISYMD,D) 1599 1600 CALL DSCAL(NCKIJ(ISBLCKIJ),HALF,WORK(KBLUMATBD),1) 1601C 1602C 1603 CALL T3_FORBIDDEN(WORK(KBLUMATBD),ISYMBL,ISYMB,B, 1604 * ISYMD,D) 1605C 1606C ------------------------------------------------ 1607C Calculate U(ci,jk) for fixed d,b for t3-bar. 1608C ------------------------------------------------ 1609 CALL DZERO(WORK(KBLUMATDB),NCKIJ(ISBLCKIJ)) 1610 1611 CALL CCFOP_UMAT(FREQL,ZL1AM,ISYMBL,ZL2TP, 1612 * ISYMBL,WORK(KXIAJB),ISYM0, 1613 * WORK(KFCKBA),WORK(KTRVI20), 1614 * WORK(KTRVI13),WORK(KTROC03), 1615 * WORK(KTROC23),ISYM0, 1616 * WORK(KFOCKD),WORK(KBLDIAG), 1617 * WORK(KBLUMATDB),WORK(KTMAT), 1618 * WORK(KEND4),LWRK4,WORK(KBLINDSQ), 1619 * LENSQBL,ISYMD,D,ISYMB,B) 1620 1621 CALL DSCAL(NCKIJ(ISBLCKIJ),HALF,WORK(KBLUMATDB),1) 1622 1623 CALL T3_FORBIDDEN(WORK(KBLUMATDB),ISYMBL,ISYMD,D, 1624 * ISYMB,B) 1625C 1626C ------------------------------------------- 1627C Sum up the S-bar and U-bar to get a real T3-bar 1628C ------------------------------------------- 1629 CALL CC3_T3BD(ISBLCKIJ,WORK(KBLSMATBD), 1630 * WORK(KBLSMATDB), 1631 * WORK(KBLUMATBD),WORK(KBLUMATDB), 1632 * WORK(KTMAT),WORK(KBLINDSQ), 1633 * LENSQBL) 1634C 1635C ------------------------------------- 1636C Write T3-bar as T3^D(ai,bj,l) to disc 1637C ------------------------------------- 1638 DO ISYML = 1, NSYM 1639 ISYMDL = MULD2H(ISYMD,ISYML) 1640 ISAIBJ = MULD2H(ISYMBL,ISYMDL) 1641 DO L = 1, NRHF(ISYML) 1642 DO ISYMJ = 1, NSYM 1643 ISYMBJ = MULD2H(ISYMJ,ISYMB) 1644 ISYMAI = MULD2H(ISAIBJ,ISYMBJ) 1645 ISYAIL = MULD2H(ISYMAI,ISYML) 1646 ISAIB = MULD2H(ISYMAI,ISYMB) 1647 DO J = 1, NRHF(ISYMJ) 1648 1649 KOFFA = ISAIKJ(ISYAIL,ISYMJ)+ NCKI(ISYAIL)*(J-1) 1650 * + ISAIK(ISYMAI, ISYML)+NT1AM(ISYMAI)*(L-1) 1651 1652 IADR = ISWTL(ISAIBJ,ISYML) + NT2SQ(ISAIBJ)*(L-1)+ 1653 * ISTLN(ISAIB,ISYMJ) + NCKATR(ISAIB)*(J-1)+ 1654 * ICKATR(ISYMAI,ISYMB)+ NT1AM(ISYMAI)*(B-1)+1 1655 1656 CALL PUTWA2(LUTBAR,FNTBAR,WORK(KTMAT+KOFFA), 1657 * IADR,NT1AM(ISYMAI)) 1658 END DO 1659 END DO 1660 END DO 1661 END DO 1662C 1663C --------------------------------------------- 1664C Initialize WMAT: 1665C --------------------------------------------- 1666 CALL DZERO(WORK(KWMAT),NCKIJ(ISWMAT)) 1667 1668C --------------------------------------------- 1669C Read and transform integrals used in second S 1670C --------------------------------------------- 1671 IOFF = ICKBD(ISYCKD,ISYMB) 1672 * + NCKATR(ISYCKD)*(B - 1) + 1 1673 IF (NCKATR(ISYCKD) .GT. 0) THEN 1674 CALL GETWA2(LUDKBC,FNDKBC,WORK(KTRVI8),IOFF, 1675 * NCKATR(ISYCKD)) 1676 ENDIF 1677 1678 1679 CALL CCSDT_SRTVIR(WORK(KTRVI8),WORK(KTRVI9), 1680 * WORK(KEND4),LWRK4,ISYMB,ISYM0) 1681C 1682C -------------------------------------------------- 1683C Read T1Y-transformed vir int (fixed B) used for W3 1684C -------------------------------------------------- 1685C 1686 IF (TOT_T3Y) THEN 1687C 1688 DO MCAU = 0, NCAU 1689C 1690 !Generate file names for T1Y-transformed integrals 1691 IF (LCAUCHY) THEN 1692 !(FNCKJDR is not needed here) 1693 CALL GET_FILENAME(FNCKJDR,FNDKBCR,MCAU) 1694 ELSE 1695 FNDKBCR = 'CCSDT_____ETAFD4' 1696 END IF 1697C 1698 !Open the files for T1Y-transformed integrals 1699 LUDKBCR = -1 1700 CALL WOPEN2(LUDKBCR,FNDKBCR,64,0) 1701C 1702 !Get the offset for a given MCAU 1703 KOFF1 = KTRVI8Y + MCAU*NCKATR(ISCKDY) 1704C 1705 !Read in the integrals 1706 IOFF = ICKBD(ISCKDY,ISYMB) + 1707 & NCKATR(ISCKDY)*(B - 1) + 1 1708 IF (NCKATR(ISCKDY) .GT. 0) THEN 1709 CALL GETWA2(LUDKBCR,FNDKBCR,WORK(KOFF1),IOFF, 1710 & NCKATR(ISCKDY)) 1711 ENDIF 1712C 1713 !Close and keep the file for T1-transformed virt int 1714 !...or delte it if you don't need it anymore 1715 1716 IF ( (ISYMD.EQ.NSYM) .AND. (ISYMB.EQ.NSYM) 1717 * .AND. (D.EQ.NVIR(ISYMD)) 1718 * .AND. (B.EQ.NVIR(ISYMB))) THEN 1719C 1720 CALL WCLOSE2(LUDKBCR,FNDKBCR,'DELETE') 1721 ELSE 1722 CALL WCLOSE2(LUDKBCR,FNDKBCR,'KEEP') 1723 END IF 1724C 1725 END DO !MCAU 1726C 1727 ENDIF 1728C 1729C ----------------------------------------- 1730C Read virtual integrals used in second U 1731C ----------------------------------------- 1732 IOFF = ICKAD(ISYCKD,ISYMB) 1733 * + NCKA(ISYCKD)*(B - 1) + 1 1734 IF (NCKA(ISYCKD) .GT. 0) THEN 1735 CALL GETWA2(LUDELD,FNDELD,WORK(KINTVI),IOFF, 1736 * NCKA(ISYCKD)) 1737 ENDIF 1738 1739 CALL CCSDT_TRVIR(WORK(KINTVI),WORK(KTRVI10), 1740 * XLAMDH,ISYMB,B,ISYM0, 1741 * WORK(KEND5),LWRK5) 1742C 1743C -------------------------------------------------- 1744C Calculate the S(ci,bk,dj) matrix for T3 for B,D. 1745C -------------------------------------------------- 1746 CALL CC3_SMAT(0.0D0,T2TP,ISYM0,WORK(KTMAT), 1747 * WORK(KTRVI0),WORK(KTRVI2), 1748 * WORK(KTROC0),ISYM0, 1749 * WORK(KFOCKD),WORK(KDIAG), 1750 * WORK(KSMATBD),WORK(KEND4),LWRK4, 1751 * WORK(KINDEXB),WORK(KINDSQ),LENSQ, 1752 * ISYMB,B,ISYMD,D) 1753 1754 CALL T3_FORBIDDEN(WORK(KSMATBD),ISYM0,ISYMB,B,ISYMD,D) 1755C 1756 IF (LOCDBG) WRITE(LUPRI,*) 'Norm of SMAT :', 1757 * DDOT(NCKIJ(ISCKIJ),WORK(KSMATBD),1,WORK(KSMATBD),1) 1758C 1759C -------------------------------------------------- 1760C Calculate the S(ci,bk,dj) matrix for T3 for D,B. 1761C -------------------------------------------------- 1762 CALL CC3_SMAT(0.0D0,T2TP,ISYM0,WORK(KTMAT), 1763 * WORK(KTRVI8),WORK(KTRVI9), 1764 * WORK(KTROC0),ISYM0, 1765 * WORK(KFOCKD),WORK(KDIAG), 1766 * WORK(KSMATDB),WORK(KEND4),LWRK4, 1767 * WORK(KINDEXD),WORK(KINDSQ),LENSQ, 1768 * ISYMD,D,ISYMB,B) 1769 1770 CALL T3_FORBIDDEN(WORK(KSMATDB),ISYM0,ISYMD,D,ISYMB,B) 1771C 1772 IF (LOCDBG) WRITE(LUPRI,*) 'Norm of SMAT3:', 1773 * DDOT(NCKIJ(ISCKIJ),WORK(KSMATDB),1,WORK(KSMATDB),1) 1774C 1775C --------------------------------- 1776C Calculate U(ci,jk) for fixed b,d. 1777C --------------------------------- 1778 CALL CC3_UMAT(0.0D0,T2TP,ISYM0,WORK(KTRVI1), 1779 * WORK(KTROC02),ISYM0,WORK(KFOCKD), 1780 * WORK(KDIAG),WORK(KUMATBD),WORK(KTMAT), 1781 * WORK(KEND4),LWRK4,WORK(KINDSQ),LENSQ, 1782 * ISYMB,B,ISYMD,D) 1783 1784 CALL T3_FORBIDDEN(WORK(KUMATBD),ISYM0,ISYMB,B,ISYMD,D) 1785C 1786 IF (LOCDBG) WRITE(LUPRI,*) 'Norm of UMAT :', 1787 * DDOT(NCKIJ(ISCKIJ),WORK(KUMATBD),1,WORK(KUMATBD),1) 1788C 1789C --------------------------------- 1790C Calculate U(ci,jk) for fixed d,b. 1791C --------------------------------- 1792 CALL CC3_UMAT(0.0D0,T2TP,ISYM0,WORK(KTRVI10), 1793 * WORK(KTROC02),ISYM0,WORK(KFOCKD), 1794 * WORK(KDIAG),WORK(KUMATDB),WORK(KTMAT), 1795 * WORK(KEND4),LWRK4,WORK(KINDSQ),LENSQ, 1796 * ISYMD,D,ISYMB,B) 1797 1798 CALL T3_FORBIDDEN(WORK(KUMATDB),ISYM0,ISYMD,D,ISYMB,B) 1799 1800 IF (LOCDBG) WRITE(LUPRI,*) 'Norm of UMAT3:', 1801 * DDOT(NCKIJ(ISCKIJ),WORK(KUMATDB),1,WORK(KUMATDB),1) 1802C 1803C ----------------------------------- 1804C Sum up the S and U to get T30 1805C ----------------------------------- 1806 CALL CC3_T3BD(ISCKIJ,WORK(KSMATBD),WORK(KSMATDB), 1807 * WORK(KUMATBD),WORK(KUMATDB), 1808 * WORK(KTMAT),WORK(KINDSQ),LENSQ) 1809C 1810C ----------------------------------- 1811C Use the T30 in TMAT to calculate W3 1812C ----------------------------------- 1813 IF ( (LISTR(1:3).EQ.'R1 ') .OR. 1814 * (LISTR(1:3).EQ.'RC ') ) THEN 1815C 1816 IF (LOCDBG) WRITE(LUPRI,*) 'Norm of TMAT:', 1817 * DDOT(NCKIJ(ISCKIJ),WORK(KTMAT),1,WORK(KTMAT),1) 1818C 1819C -------------------------------------------------------- 1820C Calculate the term <mu3|[Y,T3]|HF> virtual contribution 1821C added in W^BD(KWMAT) 1822C -------------------------------------------------------- 1823C 1824 CALL WBD_V(WORK(KTMAT),ISCKIJ,WORK(KFOCKY),ISTAMY, 1825 * WORK(KWMAT),ISWMAT,WORK(KEND4),LWRK4) 1826 1827 IF (LOCDBG) WRITE(LUPRI,*) 'Norm of WMAT (WBD_V):', 1828 * DDOT(NCKIJ(ISWMAT),WORK(KWMAT),1,WORK(KWMAT),1) 1829 1830C 1831C --------------------------------------------------------- 1832C Calculate the term <mu3|[Y,T3]|HF> occupied contribution 1833C added in W^BD(KWMAT) 1834C --------------------------------------------------------- 1835C 1836 CALL WBD_O(WORK(KTMAT),ISCKIJ,WORK(KFOCKY),ISTAMY, 1837 * WORK(KWMAT),ISWMAT,WORK(KEND4),LWRK4) 1838 1839 IF (LOCDBG) WRITE(LUPRI,*) 'Norm of WMAT (WBD_O)', 1840 * DDOT(NCKIJ(ISWMAT),WORK(KWMAT),1,WORK(KWMAT),1) 1841 1842C 1843C --------------------------------------------------------- 1844C Calculate the term <mu3|[[Y,T2],T2]|HF> 1845C added in W^BD(KWMAT) 1846C --------------------------------------------------------- 1847C 1848 CALL WBD_T2(.FALSE.,B,ISYMB,D,ISYMD, 1849 * T2TP,ISYM0,T2TP,ISYM0, 1850 * WORK(KFOCKY),ISTAMY, 1851 * WORK(KINDSQW),LENSQW, 1852 * WORK(KWMAT),ISWMAT,WORK(KEND4),LWRK4) 1853 1854 IF (LOCDBG) WRITE(LUPRI,*) 'Norm of WMAT (WBD_T2)', 1855 * DDOT(NCKIJ(ISWMAT),WORK(KWMAT),1,WORK(KWMAT),1) 1856c 1857* call sum_pt3(work(kwmat),isymb,b,isymd,d, 1858* * 1,work(kx3am),4) 1859C 1860 END IF 1861C 1862C -------------------------------------------------------------- 1863C Calculate the terms <mu3|[H^0,T2Y]|HF> and <mu3|[H^Y,T2^0]|HF> 1864C added in W^BD(KWMAT) 1865C -------------------------------------------------------------- 1866C 1867 DO MCAU = 0, NCAU !for non-Cauchy calc this is dummy loop 1868C 1869 IF (TOT_T3Y) THEN 1870C 1871 !Get offset for KT2B vec for a given MCAU 1872 KOFF2 = KT2B + MCAU*NT2SQ(ISTAMY) 1873 1874 !Calculate the term <mu3|[H^0,T2Y]|HF> 1875 1876 1877 CALL WBD_GROUND(WORK(KOFF2),ISTAMY,WORK(KTMAT), 1878 * WORK(KTRVI0),WORK(KTRVI8), 1879 * WORK(KTROC0),1, 1880 * WORK(KWMAT), 1881 * WORK(KEND4),LWRK4, 1882 * WORK(KINDSQW),LENSQW, 1883 * ISYMB,B,ISYMD,D) 1884C 1885 !Get offset for T1Y-trans occ int for a given MCAU 1886 KOFFOCC = KTROC0Y + MCAU*NTRAOC(ISTAMY) 1887C 1888 !Get offset for T1Y-trans virt int with fixed D for MCAU 1889 KOFFINTD = KTRVI0Y + MCAU*NCKATR(ISCKBY) 1890 !Get offset for T1Y-trans virt int with fixed B for MCAU 1891 KOFFINTB = KTRVI8Y + MCAU*NCKATR(ISCKDY) 1892C 1893 !Calculate the term <mu3|[H^Y,T2^0]|HF> 1894 CALL WBD_GROUND(T2TP,ISYM0,WORK(KTMAT), 1895 * WORK(KOFFINTD),WORK(KOFFINTB), 1896 * WORK(KOFFOCC),ISTAMY, 1897 * WORK(KWMAT), 1898 * WORK(KEND4),LWRK4, 1899 * WORK(KINDSQW),LENSQW, 1900 * ISYMB,B,ISYMD,D) 1901 ENDIF 1902 1903 CALL WBD_DIA(B,ISYMB,D,ISYMD,FREQY, 1904 * ISWMAT,WORK(KWMAT), 1905 * WORK(KDIAGW),WORK(KFOCKD)) 1906 1907 CALL T3_FORBIDDEN(WORK(KWMAT),ISTAMY, 1908 * ISYMB,B,ISYMD,D) 1909C 1910 END DO !MCAU 1911C 1912 !AFTER ALL THE SIGN SHOULD NOT BE TURNED !!!!! 1913 1914* IF (LCAUCHY) THEN 1915* !trun the sign C_T <- -C_T 1916* CALL DSCAL(NCKIJ(ISTAMY),-1.0D0,WORK(KWMAT),1) 1917* END IF !LCAUCHY 1918 1919* call sum_pt3(work(kwmat),isymb,b,isymd,d, 1920* * 1,work(kx3am),4) 1921 1922 IF (LOCDBG) WRITE(LUPRI,*) 'Norm of WMAT (WBD_DIA)', 1923 * DDOT(NCKIJ(ISWMAT),WORK(KWMAT),1,WORK(KWMAT),1) 1924C 1925C ------------------------------------- 1926C Write WMAT as WMAT^D(ai,bj,l) to disc 1927C ------------------------------------- 1928 DO ISYML = 1, NSYM 1929 ISYMDL = MULD2H(ISYMD,ISYML) 1930 ISAIBJ = MULD2H(ISTAMY,ISYMDL) 1931 DO L = 1, NRHF(ISYML) 1932 DO ISYMJ = 1, NSYM 1933 ISYMBJ = MULD2H(ISYMJ,ISYMB) 1934 ISYMAI = MULD2H(ISAIBJ,ISYMBJ) 1935 ISYAIL = MULD2H(ISYMAI,ISYML) 1936 DO J = 1, NRHF(ISYMJ) 1937 1938 KOFFA = ISAIKJ(ISYAIL,ISYMJ)+ NCKI(ISYAIL)*(J-1) 1939 * + ISAIK(ISYMAI, ISYML)+NT1AM(ISYMAI)*(L-1) 1940 1941 NBJ = IT1AM(ISYMB,ISYMJ)+NVIR(ISYMB)*(J-1)+B 1942 IADR = ISWTL(ISAIBJ,ISYML)+NT2SQ(ISAIBJ)*(L-1)+ 1943 * IT2SQ(ISYMAI,ISYMBJ)+NT1AM(ISYMAI)*(NBJ-1)+1 1944 1945 CALL PUTWA2(LUWMAT,FNWMAT,WORK(KWMAT+KOFFA), 1946 * IADR,NT1AM(ISYMAI)) 1947 END DO 1948 END DO 1949 END DO 1950 END DO 1951 1952 ENDDO ! B 1953 ENDDO ! ISYMB 1954 1955 DO ISYML = 1, NSYM 1956 1957 ISYMDL = MULD2H(ISYMD,ISYML) 1958 ISWMAT = MULD2H(ISYMDL,ISTAMY) 1959 ISTBAR = MULD2H(ISYMBL,ISYMDL) 1960 1961 KWMAT = KEND1 1962 KT2Y = KWMAT + NT2SQ(ISWMAT) 1963 KEND2 = KT2Y + NT1AM(ISWMAT) 1964 LWRK2 = LWORK - KEND2 1965 1966 IF ( LWRK2 .LT. 0 ) THEN 1967 CALL QUIT('Out of memory in CCSDT_ETA_FA_DEN (x)') 1968 ENDIF 1969 1970 DO L = 1, NRHF(ISYML) 1971 1972C -------------------------------------------- 1973C Read WMAT from file and generate WMAT-tilde: 1974C -------------------------------------------- 1975 IADR = ISWTL(ISWMAT,ISYML) + NT2SQ(ISWMAT)*(L-1) + 1 1976 CALL GETWA2(LUWMAT,FNWMAT,WORK(KWMAT),IADR,NT2SQ(ISWMAT)) 1977 1978 CALL CC_T2MOD(WORK(KWMAT),ISWMAT,HALF) 1979 1980C -------------------------------------------- 1981C Read response doubles amplitudes T2^y,DL: 1982C -------------------------------------------- 1983 NDL = IT1AM(ISYMD,ISYML) + NVIR(ISYMD)*(L-1) + D 1984 IADR = IT2SQ(ISWMAT,ISYMDL) + NT1AM(ISWMAT)*(NDL-1) + 1 1985 1986 CALL GETWA2(LUT2Y,FNT2Y,WORK(KT2Y),IADR,NT1AM(ISWMAT)) 1987 1988 1989C ----------------------------------- 1990C Loop over outermost occupied index: 1991C ----------------------------------- 1992 DO ISYMN = 1, NSYM 1993 ISYEMF = MULD2H(ISTBAR,ISYMN) 1994 1995 KTBAR = KEND2 1996 KEND3 = KTBAR + NCKATR(ISYEMF) 1997 LWRK3 = LWORK - KEND2 1998 1999 IF ( LWRK3 .LT. 0 ) THEN 2000 CALL QUIT('Out of memory in CCSDT_ETA_FA_DEN (x2)') 2001 ENDIF 2002 2003 DO N = 1, NRHF(ISYMN) 2004 2005C ---------------------- 2006C Read T3-BAR from file: 2007C ---------------------- 2008 IADR = ISWTL(ISTBAR,ISYML) + NT2SQ(ISTBAR)*(L-1) + 2009 * ISTLN(ISYEMF,ISYMN) + NCKATR(ISYEMF)*(N-1)+1 2010 CALL GETWA2(LUTBAR,FNTBAR,WORK(KTBAR),IADR, 2011 * NCKATR(ISYEMF)) 2012 2013C ------------------------------------------------------- 2014C D(fb) <- D(fb)+ sum_em tbar^DLN(em,f) Wtilde^DL(em,bN): 2015C ------------------------------------------------------- 2016 DO ISYMBN = 1, NSYM 2017 ISYMEM = MULD2H(ISWMAT,ISYMBN) 2018 ISYMB = MULD2H(ISYMBN,ISYMN) 2019 ISYMF = MULD2H(ISYEMF,ISYMEM) 2020 2021 KOFFT = KTBAR + ICKATR(ISYMEM,ISYMF) 2022 KBN = IT1AM(ISYMB,ISYMN)+NVIR(ISYMB)*(N-1)+1 2023 KOFFW = KWMAT + IT2SQ(ISYMEM,ISYMBN) + 2024 * NT1AM(ISYMEM)*(KBN-1) 2025 KOFFD = 1 + IMATAB(ISYMF,ISYMB) 2026 2027 NTOTEM = MAX(NT1AM(ISYMEM),1) 2028 NNVIRF = MAX(NVIR(ISYMF),1) 2029 2030 CALL DGEMM('T','N',NVIR(ISYMF),NVIR(ISYMB), 2031 * NT1AM(ISYMEM), 2032 * ONE,WORK(KOFFT),NTOTEM,WORK(KOFFW),NTOTEM, 2033 * ONE,DAB(KOFFD),NNVIRF) 2034 2035 END DO 2036 2037C ------------------------------------------------------- 2038C D(iN) <- D(iN)- sum_em tbar^DLN(em,f) Wtilde^DL(em,fi): 2039C ------------------------------------------------------- 2040 DO ISYMFI = 1, NSYM 2041 ISYMEM = MULD2H(ISWMAT,ISYMFI) 2042 ISYMF = MULD2H(ISYEMF,ISYMEM) 2043 ISYMI = MULD2H(ISYMFI,ISYMF) 2044 2045 KOFFT = KTBAR + ICKATR(ISYMEM,ISYMF) 2046 KOFFW = KWMAT + IT2SQ(ISYMEM,ISYMFI) + 2047 * NT1AM(ISYMEM)*IT1AM(ISYMF,ISYMI) 2048 KOFFD = 1 + IMATIJ(ISYMI,ISYMN) + 2049 * NRHF(ISYMI)*(N-1) 2050 2051 NNEMF = MAX(NT1AM(ISYMEM)*NVIR(ISYMF),1) 2052 2053 CALL DGEMV('T',NT1AM(ISYMEM)*NVIR(ISYMF),NRHF(ISYMI), 2054 * -ONE,WORK(KOFFW),NNEMF,WORK(KOFFT),1, 2055 * ONE,DIJ(KOFFD),1) 2056 2057 END DO 2058 2059C ------------------------------------------------------- 2060C M(imfN) <-- M(imfN) + sum_em t^DL(ei) tbar^DLN(em,f): 2061C ------------------------------------------------------- 2062 IF (DO_MMAT) THEN 2063 DO ISYMF = 1, NSYM 2064 ISYMFN = MULD2H(ISYMF,ISYMN) 2065 ISYMEM = MULD2H(ISYEMF,ISYMF) 2066 DO ISYMM = 1, NSYM 2067 ISYME = MULD2H(ISYMEM,ISYMM) 2068 ISYMI = MULD2H(ISYMDL,ISYME) 2069 ISYMIM = MULD2H(ISYMM, ISYMI) 2070 ISYMEI = MULD2H(ISYME, ISYMI) 2071 ISYEIL = MULD2H(ISYMEI,ISYML) 2072 2073 DO F = 1, NVIR(ISYMF) 2074 2075 KOFFA = 1 + IT2SP(ISYEIL,ISYMD) + 2076 * NCKI(ISYEIL) *(D-1) + ISAIK(ISYMEI,ISYML) + 2077 * NT1AM(ISYMEI)*(L-1) + IT1AM(ISYME,ISYMI) 2078 2079 KOFFT = KTBAR + ICKATR(ISYMEM,ISYMF) + 2080 * NT1AM(ISYMEM)*(F-1) + IT1AM(ISYME,ISYMM) 2081 2082 KFN = IT1AM(ISYMF,ISYMN) + NVIR(ISYMF)*(N-1) + F 2083 KOFFM = 1 + ISAIKL(ISYMFN,ISYMIM) + 2084 * NMATIJ(ISYMIM)*(KFN-1) + IMATIJ(ISYMI,ISYMM) 2085 2086 NVIRE = MAX(NVIR(ISYME),1) 2087 NRHFI = MAX(NRHF(ISYMI),1) 2088 2089 CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMM), 2090 * NVIR(ISYME), 2091 * ONE,T2TP(KOFFA),NVIRE,WORK(KOFFT),NVIRE, 2092 * ONE,XMMAT(KOFFM),NRHFI) 2093 2094 END DO 2095 END DO 2096 END DO 2097 END IF 2098 2099C ------------------------------------------------------- 2100C y^M(imfN) <- M(imfN)+ sum_em y^t^DL(ei) tbar^DLN(em,f): 2101C ------------------------------------------------------- 2102 IF (DO_YMMAT) THEN 2103 DO ISYMF = 1, NSYM 2104 ISYMFN = MULD2H(ISYMF,ISYMN) 2105 ISYMEM = MULD2H(ISYEMF,ISYMF) 2106 DO ISYMM = 1, NSYM 2107 ISYME = MULD2H(ISYMEM,ISYMM) 2108 ISYMI = MULD2H(MULD2H(ISYMDL,ISYME),ISTAMY) 2109 ISYMIM = MULD2H(ISYMM, ISYMI) 2110 ISYMEI = MULD2H(ISYME, ISYMI) 2111 ISYEIL = MULD2H(ISYMEI,ISYML) 2112 2113 DO F = 1, NVIR(ISYMF) 2114 2115 KOFFA = KT2Y + IT1AM(ISYME,ISYMI) 2116 2117 KOFFT = KTBAR + ICKATR(ISYMEM,ISYMF) + 2118 * NT1AM(ISYMEM)*(F-1) + IT1AM(ISYME,ISYMM) 2119 2120 KFN = IT1AM(ISYMF,ISYMN) + NVIR(ISYMF)*(N-1) + F 2121 KOFFM = 1 + ISAIKL(ISYMFN,ISYMIM) + 2122 * NMATIJ(ISYMIM)*(KFN-1) + IMATIJ(ISYMI,ISYMM) 2123 2124 NVIRE = MAX(NVIR(ISYME),1) 2125 NRHFI = MAX(NRHF(ISYMI),1) 2126 2127 CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMM), 2128 * NVIR(ISYME), 2129 * ONE,WORK(KOFFA),NVIRE,WORK(KOFFT),NVIRE, 2130 * ONE,YMMAT(KOFFM),NRHFI) 2131 2132 END DO 2133 END DO 2134 END DO 2135 END IF 2136 2137 IF (DO_DIA) THEN 2138C ------------------------------------------------------ 2139C D(Lb) <-- D(Lb) - sum_em tbar^DN(em) Wtilde^DL(em,bN): 2140C ------------------------------------------------------ 2141 ISYMDN = MULD2H(ISYMD,ISYMN) 2142 ISYMEM = MULD2H(ISYMBL,ISYMDN) 2143 ISYMBN = MULD2H(ISWMAT,ISYMEM) 2144 ISYMB = MULD2H(ISYMBN,ISYMN) 2145 ISYEMN = MULD2H(ISYMEM,ISYMN) 2146 2147 KOFFT = 1 + IT2SP(ISYEMN,ISYMD) + 2148 * NCKI(ISYEMN) * (D-1) + ISAIK(ISYMEM,ISYMN)+ 2149 * NT1AM(ISYMEM)* (N-1) 2150 2151 KBN = IT1AM(ISYMB,ISYMN)+NVIR(ISYMB)*(N-1)+1 2152 KOFFW = KWMAT + IT2SQ(ISYMEM,ISYMBN) + 2153 * NT1AM(ISYMEM)*(KBN-1) 2154 2155 KOFFD = 1 + IT1AM(ISYMB,ISYML)+NVIR(ISYMB)*(L-1) 2156 2157 NTOTEM = MAX(NT1AM(ISYMEM),1) 2158 2159 CALL DGEMV('T',NT1AM(ISYMEM),NVIR(ISYMB), 2160 * -ONE,WORK(KOFFW),NTOTEM,ZL2TP(KOFFT),1, 2161 * ONE,DIA(KOFFD),1) 2162 2163 2164C ----------------------------------------------------- 2165C D(mb) <-- D(mb) - sum_e tbar^DL(eN) Wtilde^DL(em,bN): 2166C ----------------------------------------------------- 2167 DO ISYMBN = 1, NSYM 2168 ISYMEM = MULD2H(ISWMAT,ISYMBN) 2169 ISYMB = MULD2H(ISYMBN,ISYMN) 2170 ISYMEN = MULD2H(ISYMBL,ISYMDL) 2171 ISYME = MULD2H(ISYMEN,ISYMN) 2172 ISYMM = MULD2H(ISYMEM,ISYME) 2173 ISYENL = MULD2H(ISYMEN,ISYML) 2174 2175 NVIRE = MAX(NVIR(ISYME),1) 2176 NVIRB = MAX(NVIR(ISYMB),1) 2177 2178 DO B = 1, NVIR(ISYMB) 2179 2180 KOFFT = 1 + IT2SP(ISYENL,ISYMD) + 2181 * NCKI(ISYENL) * (D-1)+ISAIK(ISYMEN,ISYML)+ 2182 * NT1AM(ISYMEN)* (L-1)+IT1AM(ISYME,ISYMN) + 2183 * NVIR(ISYME) * (N-1) 2184 2185 KBN = IT1AM(ISYMB,ISYMN)+NVIR(ISYMB)*(N-1)+B 2186 KOFFW = KWMAT + IT2SQ(ISYMEM,ISYMBN) + 2187 * NT1AM(ISYMEM)*(KBN-1)+IT1AM(ISYME,ISYMM) 2188 2189 KOFFD = IT1AM(ISYMB,ISYMM) + B 2190 2191 CALL DGEMV('T',NVIR(ISYME),NRHF(ISYMM), 2192 * -ONE,WORK(KOFFW),NVIRE,ZL2TP(KOFFT),1, 2193 * ONE,DIA(KOFFD),NVIRB) 2194 END DO 2195 2196 END DO 2197 END IF ! DO_DIA 2198 2199 END DO ! N 2200 END DO ! ISYMN 2201 2202 IF (DO_DIA) THEN 2203C -------------------------------------------------------- 2204C D(nb) <-- D(nb) + 2 sum_em tbar^DL(em) Wtilde^DL(em,bN): 2205C -------------------------------------------------------- 2206 ISYMEM = MULD2H(ISYMBL,ISYMDL) 2207 ISYMBN = MULD2H(ISWMAT,ISYMEM) 2208 ISYEML = MULD2H(ISYMEM,ISYML) 2209 2210 KOFFT = 1 + IT2SP(ISYEML,ISYMD) + 2211 * NCKI(ISYEML) * (D-1) + ISAIK(ISYMEM,ISYML) + 2212 * NT1AM(ISYMEM)* (L-1) 2213 2214 KOFFW = KWMAT + IT2SQ(ISYMEM,ISYMBN) 2215 2216 NTOTEM = MAX(NT1AM(ISYMEM),1) 2217 2218 CALL DGEMV('T',NT1AM(ISYMEM),NT1AM(ISYMBN), 2219 * TWO,WORK(KOFFW),NTOTEM,ZL2TP(KOFFT),1, 2220 * ONE,DIA,1) 2221 END IF ! DO_DIA 2222 2223 2224 END DO ! L 2225 END DO ! ISYML 2226 2227 ENDDO ! D 2228 ENDDO ! ISYMD 2229 2230C 2231C------------- 2232C End 2233C------------- 2234C 2235 CALL WCLOSE2(LUTBAR,FNTBAR,'DELETE') 2236 CALL WCLOSE2(LUWMAT,FNWMAT,'DELETE') 2237C 2238C--------------------------------------------------------------------- 2239C It might have happened that 'CC3_CAUINT_V*' or 'CCSDT_____ETAFD4' 2240C files have not been deleted up there. Do it now! 2241C--------------------------------------------------------------------- 2242C 2243 IF (LCAUCHY) THEN 2244 CALL DELETE_FILES('CC3_CAUINT_V*') 2245 ELSE 2246 CALL DELETE_FILES('CCSDT_____ETAFD4') 2247 END IF 2248 2249C 2250* call print_pt3(work(kx3am),1,4) 2251C 2252 CALL QEXIT('CTETAFAD') 2253C 2254 RETURN 2255 END 2256*=====================================================================* 2257 SUBROUTINE CC_T2MOD(T2AMP,ISYAMP,FAC) 2258*---------------------------------------------------------------------* 2259* 2260* Purpose: construct for a squared vector the following 2261* combination in place: 2262* 2263* T2-new(ai,bj) = T2-old(ai,bj) + FAC * T2-old(bj,ai) 2264* 2265* Written by Christof Haettig, Mai 2002, Aarhus 2266* 2267*=====================================================================* 2268 IMPLICIT NONE 2269#include "ccsdsym.h" 2270#include "ccorb.h" 2271 2272#if defined (SYS_CRAY) 2273 REAL T2AMP(*), XAIBJ, XBJAI, FAC 2274#else 2275 DOUBLE PRECISION T2AMP(*), XAIBJ, XBJAI, FAC 2276#endif 2277 2278 INTEGER ISYMAI, ISYMBJ, NAI, NBJ, NAIBJ, NBJAI, NAIAI, ISYAMP 2279 2280 CALL QENTER('CC_T2MOD') 2281 2282 DO ISYMAI = 1, NSYM 2283 ISYMBJ = MULD2H(ISYMAI,ISYAMP) 2284 2285 IF (ISYMAI.GT.ISYMBJ) THEN 2286 2287 DO NAI = 1, NT1AM(ISYMAI) 2288 DO NBJ = 1, NT1AM(ISYMBJ) 2289 2290 NAIBJ = IT2SQ(ISYMAI,ISYMBJ) + NT1AM(ISYMAI)*(NBJ-1)+NAI 2291 NBJAI = IT2SQ(ISYMBJ,ISYMAI) + NT1AM(ISYMBJ)*(NAI-1)+NBJ 2292 2293 XAIBJ = T2AMP(NAIBJ) 2294 XBJAI = T2AMP(NBJAI) 2295 2296 T2AMP(NAIBJ) = XAIBJ + FAC * XBJAI 2297 T2AMP(NBJAI) = XBJAI + FAC * XAIBJ 2298 2299 END DO 2300 END DO 2301 2302 ELSE IF (ISYMAI.EQ.ISYMBJ) THEN 2303 2304 DO NAI = 1, NT1AM(ISYMAI) 2305 DO NBJ = 1, NAI-1 2306 2307 NAIBJ = IT2SQ(ISYMAI,ISYMBJ) + NT1AM(ISYMAI)*(NBJ-1)+NAI 2308 NBJAI = IT2SQ(ISYMBJ,ISYMAI) + NT1AM(ISYMBJ)*(NAI-1)+NBJ 2309 2310 XAIBJ = T2AMP(NAIBJ) 2311 XBJAI = T2AMP(NBJAI) 2312 2313 T2AMP(NAIBJ) = XAIBJ + FAC * XBJAI 2314 T2AMP(NBJAI) = XBJAI + FAC * XAIBJ 2315 2316 END DO 2317 2318 NAIAI = IT2SQ(ISYMAI,ISYMAI) + NT1AM(ISYMAI)*(NAI-1)+NAI 2319 T2AMP(NAIAI) = (1.0D0 + FAC) * T2AMP(NAIAI) 2320 2321 END DO 2322 2323 END IF 2324 2325 END DO 2326 2327 CALL QEXIT('CC_T2MOD') 2328 RETURN 2329 END 2330*=====================================================================* 2331*=====================================================================* 2332 SUBROUTINE CCSDT_ETA_TM2(DIA,ISYDIA,XMMAT,ISMMAT,T2TP, 2333 & LUT2,FNT2,ISYMT2,IOPTT2,WORK,LWORK) 2334*---------------------------------------------------------------------* 2335* 2336* Purpose: compute density matrix contribution 2337* 2338* D(ia) <-- D(ia) + sum_fnm T(am,fn) M(imfn): 2339* 2340* Note: D(ia) is stored as D(ai) using NT1AM, IT1AM 2341* 2342* IOPTT2 = 0 : take T2 from T2TP 2343* = 1 : read T2 from LUT2 2344* 2345* Written by Christof Haettig, Mai 2002, Aarhus 2346* 2347*=====================================================================* 2348 IMPLICIT NONE 2349#include "priunit.h" 2350#include "ccsdsym.h" 2351#include "ccorb.h" 2352 2353#if defined (SYS_CRAY) 2354 REAL DIA(*), XMMAT(*), T2TP(*), WORK(*) 2355 REAL ONE 2356#else 2357 DOUBLE PRECISION DIA(*), XMMAT(*), T2TP(*), WORK(*) 2358 DOUBLE PRECISION ONE 2359#endif 2360 PARAMETER(ONE=1.0D0) 2361 2362 CHARACTER*(*) FNT2 2363 INTEGER LUT2, ISYDIA, ISMMAT, ISYMT2, LWORK, IOPTT2 2364 2365 INTEGER ISYMFN, ISYMAM, ISYMIM, KTAM, KEND1, LWRK1, ISYMN, 2366 * ISYMF, ISYAMN, NFN, IADR, ISYMM, ISYMI, ISYMA, KOFFM, 2367 * KOFFT, KOFFW, NVIRA, NRHFI, KOFFD 2368 2369 CALL QENTER('CCSDT_ETA_TM') 2370 2371 IF (MULD2H(ISYMT2,ISMMAT).NE.ISYDIA) 2372 * CALL QUIT('Symmetry mismatch in CCSDT_ETA_TM2') 2373 2374 DO ISYMFN = 1, NSYM 2375 ISYMAM = MULD2H(ISYMFN,ISYMT2) 2376 ISYMIM = MULD2H(ISYMFN,ISMMAT) 2377 2378 KTAM = 1 2379 KEND1 = KTAM + NT1AM(ISYMAM) 2380 LWRK1 = LWORK - KEND1 2381 2382 IF (LWRK1 .LT. 0) THEN 2383 CALL QUIT('Out of memory in CCSDT_ETA_TM2') 2384 ENDIF 2385 2386 DO ISYMN = 1, NSYM 2387 ISYMF = MULD2H(ISYMFN,ISYMN) 2388 ISYAMN = MULD2H(ISYMAM,ISYMN) 2389 2390 DO N = 1, NRHF(ISYMN) 2391 DO F = 1, NVIR(ISYMF) 2392 NFN = IT1AM(ISYMF,ISYMN) + NVIR(ISYMF) * (N-1) + F 2393 2394 IF (IOPTT2.EQ.0) THEN 2395 IADR = 1 + IT2SP(ISYAMN,ISYMF) + NCKI(ISYAMN) *(F-1) + 2396 * ISAIK(ISYMAM,ISYMN) + NT1AM(ISYMAM)*(N-1) 2397 CALL DCOPY(NT1AM(ISYMAM),T2TP(IADR),1,WORK(KTAM),1) 2398 ELSE IF (IOPTT2.EQ.1) THEN 2399 IADR = IT2SQ(ISYMAM,ISYMFN) + NT1AM(ISYMAM)*(NFN-1) + 1 2400 CALL GETWA2(LUT2,FNT2,WORK(KTAM),IADR,NT1AM(ISYMAM)) 2401 ELSE 2402 CALL QUIT('Illegal case IOPTT2 in CCSDT_ETA_TM') 2403 END IF 2404 2405 DO ISYMM = 1, NSYM 2406 ISYMI = MULD2H(ISYMIM,ISYMM) 2407 ISYMA = MULD2H(ISYMAM,ISYMM) 2408 2409 KOFFM = 1 + ISAIKL(ISYMFN,ISYMIM) + 2410 * NMATIJ(ISYMIM)*(NFN-1) + IMATIJ(ISYMI,ISYMM) 2411 2412 KOFFT = KTAM + IT1AM(ISYMA,ISYMM) 2413 2414 KOFFD = 1 + IT1AM(ISYMA,ISYMI) 2415 2416 NVIRA = MAX(NVIR(ISYMA),1) 2417 NRHFI = MAX(NRHF(ISYMI),1) 2418 2419 CALL DGEMM('N','T',NVIR(ISYMA),NRHF(ISYMI),NRHF(ISYMM), 2420 * ONE,WORK(KOFFT),NVIRA,XMMAT(KOFFM),NRHFI, 2421 * ONE,DIA(KOFFD),NVIRA) 2422 2423 END DO 2424 2425 END DO ! N 2426 END DO ! F 2427 2428 END DO ! ISYMN 2429 END DO ! ISYMFN 2430 2431 CALL QEXIT('CCSDT_ETA_TM2') 2432 RETURN 2433 END 2434*=====================================================================* 2435*=====================================================================* 2436 SUBROUTINE CC_FOCK_RESORT(FIJ,LOO,FIA,LOV,FAI,LVO,FAB,LVV, 2437 & FOCK,ISYFCK) 2438*---------------------------------------------------------------------* 2439* 2440* Purpose: grep occupied/occupied, occupied/virtual etc. blocks 2441* out of the complete MO Fock matrix 2442* 2443* if LOO.eq.true --> return occupied/occupied block FIJ 2444* if LOV.eq.true --> return occupied/virtual block FIA 2445* if LVO.eq.true --> return virtual/occupied block FAI 2446* if LVV.eq.true --> return virtual/virtual block FAB 2447* 2448* 2449* Written by Christof Haettig, Mai 2002, Aarhus 2450* 2451*=====================================================================* 2452 IMPLICIT NONE 2453#include "ccsdsym.h" 2454#include "ccorb.h" 2455 2456#if defined (SYS_CRAY) 2457 REAL FIA(*), FIJ(*), FAI(*), FAB(*), FOCK(*) 2458#else 2459 DOUBLE PRECISION FIA(*), FIJ(*), FAI(*), FAB(*), FOCK(*) 2460#endif 2461 2462 LOGICAL LOV, LOO, LVV, LVO 2463 INTEGER ISYFCK, ISYMK, ISYMC, ISYMI, ISYMJ, ISYMA, ISYMB, 2464 & KOFF1, KOFF2 2465 2466 CALL QENTER('CC_FOCK_RESORT') 2467 2468 2469 IF (LOV) THEN 2470 DO ISYMC = 1,NSYM 2471 ISYMK = MULD2H(ISYMC,ISYFCK) 2472 DO K = 1,NRHF(ISYMK) 2473 DO C = 1,NVIR(ISYMC) 2474 KOFF1 = IFCVIR(ISYMK,ISYMC) + NORB(ISYMK)*(C - 1) + K 2475 KOFF2 = IT1AM(ISYMC,ISYMK) + NVIR(ISYMC)*(K - 1) + C 2476 FIA(KOFF2) = FOCK(KOFF1) 2477 ENDDO 2478 ENDDO 2479 ENDDO 2480 END IF 2481 2482 IF (LVO) THEN 2483 DO ISYMC = 1,NSYM 2484 ISYMK = MULD2H(ISYMC,ISYFCK) 2485 DO K = 1,NRHF(ISYMK) 2486 KOFF1 = 1+IFCRHF(ISYMK,ISYMC)+NORB(ISYMC)*(K-1)+NRHF(ISYMK) 2487 KOFF2 = 1+IT1AM(ISYMC,ISYMK) +NVIR(ISYMC)*(K-1) 2488 CALL DCOPY(NVIR(ISYMC),FOCK(KOFF1),1,FAI(KOFF2),1) 2489 ENDDO 2490 ENDDO 2491 END IF 2492 2493 2494 IF (LVV) THEN 2495 DO ISYMB = 1,NSYM 2496 ISYMA = MULD2H(ISYMB,ISYFCK) 2497 DO B = 1,NVIR(ISYMB) 2498 KOFF1 = 1+IFCVIR(ISYMA,ISYMB)+NORB(ISYMA)*(B-1)+NRHF(ISYMA) 2499 KOFF2 = 1+IMATAB(ISYMA,ISYMB)+NVIR(ISYMA)*(B-1) 2500 CALL DCOPY(NVIR(ISYMA),FOCK(KOFF1),1,FAB(KOFF2),1) 2501 END DO 2502 END DO 2503 END IF 2504 2505 IF (LOO) THEN 2506 DO ISYMJ = 1,NSYM 2507 ISYMI = MULD2H(ISYMJ,ISYFCK) 2508 DO J = 1,NRHF(ISYMJ) 2509 KOFF1 = 1 + IFCRHF(ISYMI,ISYMJ) + NORB(ISYMI)*(J - 1) 2510 KOFF2 = 1 + IMATIJ(ISYMI,ISYMJ) + NRHF(ISYMI)*(J - 1) 2511 CALL DCOPY(NRHF(ISYMI),FOCK(KOFF1),1,FIJ(KOFF2),1) 2512 END DO 2513 END DO 2514 END IF 2515 2516 CALL QEXIT('CC_FOCK_RESORT') 2517 RETURN 2518 END 2519*=====================================================================* 2520*=====================================================================* 2521 SUBROUTINE CC_XIETA_DENPREP(IXETRAN,NXETRAN,MXVEC, 2522 & IXDOTS,LISTDPX, 2523 & IEDOTS,LISTDPE,LISTL, 2524 & IDXL_XIDEN,NXIDEN,MXXIDEN, 2525 & IDXL_EADEN,IDXR_EADEN,NEADEN,MXEADEN, 2526 & IEASTEND,NEALEFT,MXEALEFT) 2527*---------------------------------------------------------------------* 2528* 2529* Purpose: set up loops for the calculation of effective Xi/Eta 2530* densities 2531* 2532* Written by Christof Haettig, Mai 2002, Aarhus 2533* 2534*=====================================================================* 2535 IMPLICIT NONE 2536#include "ccsdsym.h" 2537#include "ccorb.h" 2538#include "ccroper.h" 2539#include "cclists.h" 2540 2541* input: 2542 CHARACTER*3 LISTL, LISTDPE, LISTDPX 2543 INTEGER NXETRAN, IORDER, MXVEC 2544 INTEGER IXETRAN(MXDIM_XEVEC,NXETRAN) 2545 INTEGER IXDOTS(MXVEC,NXETRAN), IEDOTS(MXVEC,NXETRAN) 2546 2547* output for Xi densities: 2548 INTEGER NXIDEN, MXXIDEN 2549 INTEGER IDXL_XIDEN(MXXIDEN) 2550 2551* output for Eta densities: 2552 INTEGER NEADEN, MXEADEN, NEALEFT, MXEALEFT 2553 INTEGER IDXL_EADEN(MXEADEN) 2554 INTEGER IDXR_EADEN(MXEADEN) 2555 INTEGER IEASTEND(2,MXEALEFT) 2556 2557* local variables: 2558 CHARACTER*8 LABELH 2559 LOGICAL LTWOEL, LRELAX, SKIPXI, SKIPETA, CHANGES 2560 INTEGER IVEC, ITRAN, IOPER, IDLSTL, IRELAX, ISYHOP, ISYCTR, 2561 & ISYETA, IDLSTR, ISYMR, ISYML, IDX, IEADEN, IXIDEN 2562 2563* external functions: 2564 INTEGER ILSTSYM 2565 2566 CALL QENTER('CXIETADP') 2567 2568 2569 NEADEN = 0 2570 NXIDEN = 0 2571 2572 DO ITRAN = 1, NXETRAN 2573 IOPER = IXETRAN(1,ITRAN) ! operator index 2574 IDLSTL = IXETRAN(2,ITRAN) ! Left vector index 2575 IRELAX = IXETRAN(5,ITRAN) ! flag for relax. contrib. 2576 2577 ISYHOP = ISYOPR(IOPER) ! symmetry of O operator 2578 LABELH = LBLOPR(IOPER) ! operator label 2579 LTWOEL = LPDBSOP(IOPER) ! two-electron contribution to O? 2580 2581 SKIPXI = ( IXETRAN(3,ITRAN) .EQ. -1 ) 2582 SKIPETA = ( IXETRAN(4,ITRAN) .EQ. -1 ) 2583 2584 LRELAX = LTWOEL .OR. (IRELAX.GE.1) 2585 IF (LRELAX) CALL QUIT( 2586 & 'Relaxed perturbations not yet implemented CC_XIETA_DENPREP') 2587 2588 2589C --------------------------------------------------------------- 2590C set up list of effecitive Eta{O} x R or L x A{O} x R densities: 2591C --------------------------------------------------------------- 2592 IF (.NOT.SKIPETA) THEN 2593 ISYCTR = ILSTSYM(LISTL,IDLSTL) ! sym. of Left (ZETA) vector 2594 ISYETA = MULD2H(ISYHOP,ISYCTR) ! sym. of ETA result vect. 2595 2596 IVEC = 1 2597 DO WHILE(IEDOTS(IVEC,ITRAN).NE.0 .AND. IVEC.LE.MXVEC) 2598 IDLSTR = IEDOTS(IVEC,ITRAN) 2599 ISYMR = ILSTSYM(LISTDPE,IDLSTR) 2600 2601 IF (ISYMR.EQ.ISYETA) THEN 2602 2603 ! check if respective density already on our to-do list 2604 IEADEN = 0 2605 DO IDX = 1, NEADEN 2606 IF (IDLSTL.EQ.IDXL_EADEN(IDX) .AND. 2607 & IDLSTR.EQ.IDXR_EADEN(IDX) ) THEN 2608 IEADEN = IDX 2609 END IF 2610 END DO 2611 2612 ! not found, put on to-do list for densities 2613 IF (IEADEN.EQ.0) THEN 2614 NEADEN = NEADEN + 1 2615 IF (NEADEN.GT.MXEADEN) CALL QUIT('NEADEN out of range') 2616 IDXL_EADEN(NEADEN) = IDLSTL 2617 IDXR_EADEN(NEADEN) = IDLSTR 2618 END IF 2619 2620 END IF 2621 2622 IVEC = IVEC + 1 2623 END DO 2624 END IF 2625 2626C -------------------------------------------- 2627C set up list of effecitive L Xi{O} densities: 2628C -------------------------------------------- 2629 IF (.NOT.SKIPXI) THEN 2630 IVEC = 1 2631 DO WHILE(IXDOTS(IVEC,ITRAN).NE.0 .AND. IVEC.LE.MXVEC) 2632 IDLSTL = IXDOTS(IVEC,ITRAN) 2633 ISYML = ILSTSYM(LISTDPX,IDLSTL) 2634 2635 IF (ISYML.EQ.ISYHOP) THEN 2636 2637 ! check if respective density already on our to-do list 2638 IXIDEN = 0 2639 DO IDX = 1, NXIDEN 2640 IF (IDLSTL.EQ.IDXL_XIDEN(IDX)) THEN 2641 IXIDEN = IDX 2642 END IF 2643 END DO 2644 2645 ! not found, put on to-do list for densities 2646 IF (IXIDEN.EQ.0) THEN 2647 NXIDEN = NXIDEN + 1 2648 IF (NXIDEN.GT.MXXIDEN) CALL QUIT('NXIDEN out of range') 2649 IDXL_XIDEN(NXIDEN) = IDLSTL 2650 END IF 2651 2652 END IF 2653 2654 IVEC = IVEC + 1 2655 END DO 2656 END IF 2657 2658 END DO 2659 2660C ------------------------------------------------------------- 2661C sort list of effecitive Eta{O} x R or L x A{O} x R densities: 2662C ------------------------------------------------------------- 2663 CHANGES = .TRUE. 2664 DO WHILE(CHANGES) 2665 CHANGES = .FALSE. 2666 DO IEADEN = 2, NEADEN 2667 IF (IDXL_EADEN(IEADEN-1).GT.IDXL_EADEN(IEADEN)) THEN 2668 CHANGES = .TRUE. 2669 IDLSTL = IDXL_EADEN(IEADEN) 2670 IDLSTR = IDXR_EADEN(IEADEN) 2671 IDXL_EADEN(IEADEN) = IDXL_EADEN(IEADEN-1) 2672 IDXR_EADEN(IEADEN) = IDXR_EADEN(IEADEN-1) 2673 IDXL_EADEN(IEADEN-1) = IDLSTL 2674 IDXR_EADEN(IEADEN-1) = IDLSTR 2675 END IF 2676 END DO 2677 END DO 2678 2679C ------------------------------------------------------- 2680C set start and end indeces for loop over left vectors in 2681C calculation of effective Eta densities: 2682C ------------------------------------------------------- 2683 IF (NEADEN.GE.1) THEN 2684 NEALEFT = 1 2685 IEASTEND(1,1) = 1 2686 IEASTEND(2,1) = 1 2687 DO IEADEN = 2, NEADEN 2688 IF (IDXL_EADEN(IEADEN-1).NE.IDXL_EADEN(IEADEN)) THEN 2689 NEALEFT = NEALEFT + 1 2690 IEASTEND(1,NEALEFT) = IEADEN 2691 END IF 2692 IEASTEND(2,NEALEFT) = IEADEN 2693 END DO 2694 ELSE 2695 NEALEFT = 0 2696 IEASTEND(1,1) = 0 2697 IEASTEND(2,1) = 0 2698 END IF 2699 2700 CALL QEXIT('CXIETADP') 2701 RETURN 2702 END 2703*=====================================================================* 2704*=====================================================================* 2705#if defined (SYS_CRAY) 2706 REAL 2707#else 2708 DOUBLE PRECISION 2709#endif 2710 & FUNCTION CC_XIETA_DENCON(IDENS,LABELH,ISYHOP,XLAMDP,XLAMDH, 2711 & LUDEFF,FNDEFF,M2BST,WORK,LWORK) 2712*---------------------------------------------------------------------* 2713* 2714* Purpose: contract effective density with one-electron integrals 2715* 2716* Written by Christof Haettig, Mai 2002, Friedrichstal 2717* 2718*=====================================================================* 2719 IMPLICIT NONE 2720#include "ccsdsym.h" 2721#include "ccorb.h" 2722#include "priunit.h" 2723#include "dummy.h" 2724 2725 LOGICAL LOCDBG 2726 PARAMETER (LOCDBG = .FALSE.) 2727 2728 CHARACTER LABELH*(*), FNDEFF*(*) 2729 INTEGER IDENS, ISYHOP, LUDEFF, M2BST, LWORK 2730 2731#if defined (SYS_CRAY) 2732 REAL 2733#else 2734 DOUBLE PRECISION 2735#endif 2736 & WORK(LWORK), XLAMDP(*), XLAMDH(*), DDOT 2737 2738 INTEGER KDAB,KDIJ,KDIA,KEND1,KFOCK,KFOCKIJ,KFOCKIA,KFOCKAB, 2739 & LWRK1, IADRIJ, IADRAB, IADRIA, IRREP, IMAT, IERR 2740 2741 CALL QENTER('CXIETDC') 2742 2743 KDAB = 1 2744 KDIJ = KDAB + NMATAB(ISYHOP) 2745 KDIA = KDIJ + NMATIJ(ISYHOP) 2746 KEND1 = KDIA + NT1AM(ISYHOP) 2747 2748 KFOCK = KEND1 2749 KEND1 = KFOCK + N2BST(ISYHOP) 2750 2751 KFOCKIJ = KEND1 2752 KFOCKIA = KFOCKIJ + NMATIJ(ISYHOP) 2753 KFOCKAB = KFOCKIA + NT1AM(ISYHOP) 2754 KEND1 = KFOCKAB + NMATAB(ISYHOP) 2755 2756 LWRK1 = LWORK - KEND1 2757 IF (LWRK1.LT.0) CALL QUIT('Out of memory in CC_XIETA_DENCON') 2758 2759 IADRIJ = 1 + M2BST*(IDENS-1) 2760 CALL GETWA2(LUDEFF,FNDEFF,WORK(KDIJ),IADRIJ,NMATIJ(ISYHOP)) 2761 2762 IADRAB = IADRIJ + NMATIJ(ISYHOP) 2763 CALL GETWA2(LUDEFF,FNDEFF,WORK(KDAB),IADRAB,NMATAB(ISYHOP)) 2764 2765 IADRIA = IADRAB + NMATAB(ISYHOP) 2766 CALL GETWA2(LUDEFF,FNDEFF,WORK(KDIA),IADRIA,NT1AM(ISYHOP)) 2767 2768 CALL CCPRPAO(LABELH,.TRUE.,WORK(KFOCK),IRREP,IMAT,IERR, 2769 * WORK(KEND1),LWRK1) 2770 IF ((IERR.GT.0) .OR. (IERR.EQ.0 .AND. IRREP.NE.ISYHOP)) THEN 2771 WRITE(LUPRI,*) 'ISYHOP:',ISYHOP 2772 WRITE(LUPRI,*) 'IRREP :',IRREP 2773 WRITE(LUPRI,*) 'IERR :',IERR 2774 WRITE(LUPRI,*) 'LABEL :',LABELH 2775 CALL QUIT('CC_XIETA_DENCON: error reading operator ') 2776 ELSE IF (IERR.LT.0) THEN 2777 CALL DZERO(WORK(KFOCK),N2BST(ISYHOP)) 2778 END IF 2779 2780 ! transform property integrals to Lambda-MO basis 2781 CALL CC_FCKMO(WORK(KFOCK),XLAMDP,XLAMDH,WORK(KEND1),LWRK1, 2782 * ISYHOP,1,1) 2783 2784 CALL CC_FOCK_RESORT(WORK(KFOCKIJ),.TRUE.,WORK(KFOCKIA),.TRUE., 2785 & DUMMY,.FALSE., WORK(KFOCKAB),.TRUE., 2786 & WORK(KFOCK),ISYHOP) 2787 2788 CC_XIETA_DENCON = 2789 & DDOT( NT1AM(ISYHOP),WORK(KDIA),1,WORK(KFOCKIA),1) + 2790 & DDOT(NMATIJ(ISYHOP),WORK(KDIJ),1,WORK(KFOCKIJ),1) + 2791 & DDOT(NMATAB(ISYHOP),WORK(KDAB),1,WORK(KFOCKAB),1) 2792 2793 IF (LOCDBG) THEN 2794 WRITE(LUPRI,*) 'CC_XIETA_DENCON>', ISYHOP,LABELH 2795 WRITE(LUPRI,*) 'NORM^2(DIJ):', 2796 & DDOT(NMATIJ(ISYHOP),WORK(KDIJ),1,WORK(KDIJ),1) 2797 WRITE(LUPRI,*) 'NORM^2(DAB):', 2798 & DDOT(NMATAB(ISYHOP),WORK(KDAB),1,WORK(KDAB),1) 2799 WRITE(LUPRI,*) 'NORM^2(DIA):', 2800 & DDOT(NT1AM(ISYHOP),WORK(KDIA),1,WORK(KDIA),1) 2801 WRITE(LUPRI,*) 'NORM^2(FOCKIJ):', 2802 & DDOT(NMATIJ(ISYHOP),WORK(KFOCKIJ),1,WORK(KFOCKIJ),1) 2803 WRITE(LUPRI,*) 'NORM^2(KFOCKAB):', 2804 & DDOT(NMATAB(ISYHOP),WORK(KFOCKAB),1,WORK(KFOCKAB),1) 2805 WRITE(LUPRI,*) 'NORM^2(FOCKIA):', 2806 & DDOT(NT1AM(ISYHOP),WORK(KFOCKIA),1,WORK(KFOCKIA),1) 2807 WRITE(LUPRI,*) 'DIA x FOCKIA:', 2808 & DDOT(NT1AM(ISYHOP),WORK(KDIA),1,WORK(KFOCKIA),1) 2809 WRITE(LUPRI,*) 'DIJ x FOCKIJ:', 2810 & DDOT(NMATIJ(ISYHOP),WORK(KDIJ),1,WORK(KFOCKIJ),1) 2811 WRITE(LUPRI,*) 'DAB x FOCKAB:', 2812 & DDOT(NMATAB(ISYHOP),WORK(KDAB),1,WORK(KFOCKAB),1) 2813 END IF 2814 2815 CALL QEXIT('CXIETDC') 2816 RETURN 2817 END 2818*=====================================================================* 2819