1! 2! Dalton, a molecular electronic structure program 3! Copyright (C) by the authors of Dalton. 4! 5! This program is free software; you can redistribute it and/or 6! modify it under the terms of the GNU Lesser General Public 7! License version 2.1 as published by the Free Software Foundation. 8! 9! This program is distributed in the hope that it will be useful, 10! but WITHOUT ANY WARRANTY; without even the implied warranty of 11! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 12! Lesser General Public License for more details. 13! 14! If a copy of the GNU LGPL v2.1 was not distributed with this 15! code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html. 16! 17! 18C 19c /* deck cc_famat */ 20*=====================================================================* 21 SUBROUTINE CC_FAMAT( LABELA, ISYMTA, ! inp: label/symmetry A 22 & LISTB, ITAMPB, ! inp: B resp. amplit. 23 & LISTC, IZETVC, ! inp: C resp. zeta vec. 24 & IOPTRES,IFATRAN, ! output option 25 & IFADOTS,FACON, ! indeces/dotproducts 26 & FILFA, ITRAN, ! list,index for dotprods 27 & NFATRAN,MXVEC, ! dimensions for dotprods 28 & NODDY_CCSDT, ! flag for noddy triples 29 & WORK, LWORK )! work space 30*---------------------------------------------------------------------* 31* 32* Purpose: transformation of a response vector with a F matrix 33* where the hamiltonian has been substituted by a 34* perturbation operator 35* 36* F^C{A} * t^B = <lambda^C|[[A,t^B],tau_nu]|CC> 37* 38* 39* symmetries/variables: 40* 41* ISYRES : result vector GAMMA1, GAMMA2 42* ISYCTR : lagrangian multipliers (zeta vector) CTR1, CTR2 43* ISYMTA : A perturbation 44* ISYMTB : B response amplitudes 45* 46* Note: the single and double excitation parts of the result GAMMA2 47* are returned at the beginning of the work space in 48* WORK(1)... WORK(NT1AM(ISYRES)) 49* WORK(NT1AM(ISYRES)+1)... WORK(NT1AM(ISYRES)+NT2AM(ISYRES)) 50* (double excitation part will be stored in packed form) 51* 52* Written by Christof Haettig, October 1996. 53* CC3 noddy version, April 2002, Christof Haettig 54* Added CCR12, June 2005, Christian Neiss 55* 56*=====================================================================* 57#if defined (IMPLICIT_NONE) 58 IMPLICIT NONE 59#else 60# include "implicit.h" 61#endif 62#include "priunit.h" 63#include "iratdef.h" 64#include "cbieri.h" 65#include "mxcent.h" 66#include "eribuf.h" 67#include "maxorb.h" 68#include "distcl.h" 69#include "ccorb.h" 70#include "ccisao.h" 71#include "ccsdsym.h" 72#include "ccsdinp.h" 73#include "cclists.h" 74#include "r12int.h" 75#include "ccr12int.h" 76 77* local parameters: 78 CHARACTER MSGDBG*(18) 79 PARAMETER (MSGDBG='[debug] CC_FAMAT> ') 80#if defined (SYS_CRAY) 81 REAL ONE, TWO, THREE 82#else 83 DOUBLE PRECISION ONE, TWO, THREE 84#endif 85 PARAMETER (ONE = 1.0d0, TWO = 2.0d0, THREE = 3.0d0) 86 LOGICAL LOCDBG 87 PARAMETER (LOCDBG = .FALSE.) 88 INTEGER KDUM 89 PARAMETER (KDUM = +99 999 999) ! dummy address 90 91 CHARACTER*8 LABELA 92 CHARACTER*10 MODEL, MODELW 93 CHARACTER LISTB*(*), LISTC*(*), FILFA*(*) 94 LOGICAL NODDY_CCSDT 95 INTEGER ISYRES, ISYCTR, ISYMTA, ISYMTB, LWORK 96 INTEGER ITAMPB, IZETVC 97 98 INTEGER ISYTATB, ISYMJ, ISYMB, ISYMXY, ISYMI, ISYMA, IRREP, ISYM 99 INTEGER KBTAOO, KBTAVV, KPERTA, KT1AMPB, KT2AMPB, KCTR1, KCTR2 100 INTEGER KCTMO, KT1AMP0, KLAMDP0, KLAMDH0, KOFF1, KOFF2, KSCR 101 INTEGER KEND1, KEND2, KEND3, KEND1A, KXBMAT, KYBMAT, IERR 102 INTEGER LEND1, LEND2, LEND3, LEND1A, KGAMMA1, KGAMMA2, KEND0 103 INTEGER IOPT, MAXJ, NIJ, NJI, NAB, NBA, KEMAT1, KEMAT2, IVEC, 104 & LEND0, KGAMMA1EFF, KGAMMA2EFF, IOPTW, IOPTWE, IFILE 105 106 INTEGER NFATRAN, MXVEC, IOPTRES, ITRAN 107 INTEGER IFATRAN(MXDIM_FATRAN,NFATRAN), IFADOTS(MXVEC,NFATRAN) 108 INTEGER KGAMMAR12, KCTRR12, KT12AMB, KXINTTRI, KXINTSQ, IDUM 109 INTEGER LUNIT, IAN, KLAMDPB, KLAMDHB, LENMOD, IOPTWR12 110 111#if defined (SYS_CRAY) 112 REAL FREQB, FREQA 113 REAL SWAP, DUMMY 114 REAL WORK(LWORK), FACON(MXVEC,NFATRAN) 115#else 116 DOUBLE PRECISION FREQB, FREQA 117 DOUBLE PRECISION SWAP, DUMMY 118 DOUBLE PRECISION WORK(LWORK), FACON(MXVEC,NFATRAN) 119#endif 120 121 INTEGER ILSTSYM 122 CHARACTER*(3) APROXR12 123 124 125 CALL QENTER('CC_FAMAT') 126*---------------------------------------------------------------------* 127* begin: 128*---------------------------------------------------------------------* 129 IF ( .not. (CCS .or. CC2 .or. CCSD .or. CC3) ) THEN 130 WRITE (LUPRI,'(/1x,a)') 'CC_FAMAT called for a Coupled Cluster ' 131 & //'method not implemented in CC_FAMAT...' 132 CALL QUIT('Unknown CC method in CC_FAMAT.') 133 END IF 134 135*---------------------------------------------------------------------* 136* set & check symmetries: 137*---------------------------------------------------------------------* 138 ISYMTB = ILSTSYM(LISTB,ITAMPB) ! B 139 ISYCTR = ILSTSYM(LISTC,IZETVC) ! C 140 ISYTATB = MULD2H(ISYMTA,ISYMTB) ! A x B 141 ISYRES = MULD2H(ISYCTR,ISYTATB) ! A x B x C 142 143 IF (LOCDBG) THEN 144 WRITE (LUPRI,*) 'LISTB,ITAMPB,ISYMTB:',LISTB,ITAMPB,ISYMTB 145 WRITE (LUPRI,*) 'LISTC,IZETVC,ISYCTR:',LISTC,IZETVC,ISYCTR 146 END IF 147 148 149 IF (ISYMOP .NE. 1) THEN 150 WRITE (LUPRI,'(/1x,a)') 'non-total-symmetric MO integrals?!... ' 151 & //'CCLR_G has never been debugged for that!...' 152 END IF 153 154 IF (MULD2H(ISYCTR,ISYTATB) .NE. ISYRES) THEN 155 CALL QUIT('Symmetry mismatch in CC_FAMAT.') 156 END IF 157 158*---------------------------------------------------------------------* 159* flush print unit 160*---------------------------------------------------------------------* 161 Call FLSHFO(LUPRI) 162 163 IF (LOCDBG) THEN 164 WRITE (LUPRI,'(/1x,a,i15)') 'work space in CC_FAMAT:',LWORK 165 END IF 166*---------------------------------------------------------------------* 167* initialize pointer for work space and allocate memory for 168* 1) single excitation part of the result vector 169* 2) one-index transformed perturbation integrals A^B (occ/occ block) 170* 3) one-index transformed perturbation integrals A^B (vir/vir block) 171* 4) perturbation integrals A 172* 5) singles part of response amplitudes T1^B 173* 6) singles part of zeroth order lagrangian multipliers 174*---------------------------------------------------------------------* 175 KGAMMA1 = 1 176 KEND0 = KGAMMA1 + NT1AM(ISYRES) 177 LEND0 = LWORK - KEND0 178 179 KBTAOO = KEND0 180 KBTAVV = KBTAOO + NMATIJ(ISYTATB) 181 KPERTA = KBTAVV + NMATAB(ISYTATB) 182 KT1AMPB = KPERTA + NT1AM(ISYMTA) 183 KCTR1 = KT1AMPB + NT1AM(ISYMTB) 184 KEND1 = KCTR1 + NT1AM(ISYCTR) 185 LEND1 = LWORK - KEND1 186 187 IF (LEND1 .LT. 0) THEN 188 CALL QUIT('Insufficient work space in CC_FAMAT.') 189 END IF 190 191*---------------------------------------------------------------------* 192* initialize single excitation part of result vector GAMMA1: 193*---------------------------------------------------------------------* 194 Call DZERO (WORK(KGAMMA1), NT1AM(ISYRES)) 195 196*---------------------------------------------------------------------* 197* for CCS and zeroth-order zeta vector all contributions vanish: 198*---------------------------------------------------------------------* 199 IF (CCS .AND. LISTC(1:2).EQ.'L0') GOTO 9999 200 201*---------------------------------------------------------------------* 202* read singles parts for B response amplitudes and zeta vector: 203*---------------------------------------------------------------------* 204 IOPT = 1 205 CALL CC_RDRSP(LISTB,ITAMPB,ISYMTB,IOPT,MODEL, 206 & WORK(KT1AMPB),WORK(KDUM) ) 207 208 209 IOPT = 1 210 Call CC_RDRSP(LISTC,IZETVC,ISYCTR,IOPT,MODEL, 211 & WORK(KCTR1),WORK(KDUM)) 212 213 IF (LOCDBG) THEN 214 CAll AROUND('response T amplitudes B:') 215 WRITE (LUPRI,*) 'LIST/INDEX:',LISTB,ITAMPB 216 WRITE (LUPRI,*) 'Symmetry: ',ISYMTB 217 CAll CC_PRP(WORK(KT1AMPB),WORK(KDUM),ISYMTB,1,0) 218 CALL AROUND('CC lagrange multipliers') 219 CALL CC_PRP(WORK(KCTR1), WORK(KDUM), ISYCTR, 1, 0) 220 END IF 221 222*---------------------------------------------------------------------* 223* read & resort one-electron integrals for operator A: 224*---------------------------------------------------------------------* 225 KCTMO = KEND1 226 KT1AMP0 = KCTMO + N2BST(ISYMTA) 227 KLAMDP0 = KT1AMP0 + NT1AM(ISYMOP) 228 KLAMDH0 = KLAMDP0 + NLAMDT 229 KEND1A = KLAMDH0 + NLAMDT 230 LEND1A = LWORK - KEND1A 231 232 IF (LEND1A .LT. 0) THEN 233 CALL QUIT('Insufficient work space in CC_FAMAT.') 234 END IF 235 236* read the AO integrals: 237 CALL CCPRPAO(LABELA,.TRUE.,WORK(KCTMO),IRREP,ISYM,IERR, 238 & WORK(KEND1A),LEND1A) 239 IF (IERR.NE.0 .OR. IRREP.NE.ISYMTA) THEN 240 CALL QUIT('CC_FAMAT: error while reading operator '//LABELA) 241 END IF 242 243 244* get MO coefficients: 245 CALL DZERO(WORK(KT1AMP0),NT1AMX) 246 CALL LAMMAT(WORK(KLAMDP0),WORK(KLAMDH0),WORK(KT1AMP0), 247 & WORK(KEND1A),LEND1A) 248 249* transform one-electron integrals in place: 250 CALL CC_FCKMO(WORK(KCTMO),WORK(KLAMDP0),WORK(KLAMDH0), 251 & WORK(KEND1A),LEND1A,ISYMTA,1,1) 252 253* resort occupied/virtual block to T1 like storage: 254 CALL DZERO(WORK(KPERTA),NT1AM(ISYMTA)) 255 DO ISYMJ = 1, NSYM 256 ISYMB = MULD2H(ISYMJ,ISYMTA) 257 258 DO J = 1, NRHF(ISYMJ) 259 DO B = 1, NVIR(ISYMB) 260 KOFF1 = IT1AM(ISYMB,ISYMJ) + NVIR(ISYMB)*(J-1) + B 261 KOFF2 = IFCVIR(ISYMJ,ISYMB) + NORB(ISYMJ)*(B-1) + J 262 263 WORK(KPERTA-1+KOFF1) = WORK(KCTMO-1+KOFF2) 264 END DO 265 END DO 266 END DO 267 268 IF (LOCDBG) THEN 269 WRITE (LUPRI,*) MSGDBG,' A integrals in MO basis:' 270 WRITE (LUPRI,*) MSGDBG,' label, symmetry:',LABELA,ISYMTA 271 Call CC_PRP(WORK(KPERTA),WORK(KDUM),ISYMTA,1,0) 272 END IF 273 274*---------------------------------------------------------------------* 275* calculate A perturbation integrals one-index transformed with 276* the B response amplitudes T1^B: 277*---------------------------------------------------------------------* 278* occ/occ block: 279 Call CCG_1ITROO(WORK(KBTAOO), ISYTATB, 280 & WORK(KPERTA), ISYMTA, 281 & WORK(KT1AMPB),ISYMTB ) 282 283* vir/vir block: 284 Call CCG_1ITRVV(WORK(KBTAVV), ISYTATB, 285 & WORK(KPERTA), ISYMTA, 286 & WORK(KT1AMPB),ISYMTB ) 287 288*=====================================================================* 289* CCS part: < Zeta_1 | [tA^B, tau_1] | HF> 290*=====================================================================* 291* do one-index transformation with Zeta vector: 292 IOPT = 2 293 Call CCG_1ITRVO(WORK(KGAMMA1),ISYRES,WORK(KBTAOO),WORK(KBTAVV), 294 & ISYTATB,WORK(KCTR1),ISYCTR,IOPT ) 295 296 IF (LOCDBG) THEN 297 WRITE (LUPRI,*) MSGDBG,'one-index trans. A (occ/occ block):' 298 WRITE (LUPRI,'(5f12.6)') (WORK(KBTAOO-1+I),I=1,NMATIJ(ISYTATB)) 299 WRITE (LUPRI,*) MSGDBG,'one-index trans. A (vir/vir block):' 300 WRITE (LUPRI,'(2f12.6)') (WORK(KBTAVV-1+I),I=1,NMATAB(ISYTATB)) 301 WRITE (LUPRI,*) MSGDBG, 302 * 'contrib. of one-index trans. A to GAMMA:' 303 Call CC_PRP(WORK(KGAMMA1),WORK(KDUM),ISYRES,1,0) 304 END IF 305 306*---------------------------------------------------------------------* 307* end of CCS part 308*---------------------------------------------------------------------* 309 310 If (CCS) GOTO 9999 311 312*=====================================================================* 313* CC2/CCSD part for the singles: <Zeta_2| [[A,T2B], tau_1] |HF> 314*=====================================================================* 315 316*---------------------------------------------------------------------* 317* memory allocation: 318* 1) double excitation part of response amplitudes T2B (packed) 319* 2) double excitation part of zeta vector (squared) 320* 3) double excitation part of zeta vector (packed) 321* N.B. we account here for the fact, that the packed double excitation 322* part of the result vector will be returned at the beginning of the 323* work space, so we make sure, that there is enough space before 324* the zeta vector to store there later on GAMMA2 325*---------------------------------------------------------------------* 326 KT2AMPB = KEND1 327 KCTR2 = KT2AMPB + MAX( NT2AM(ISYMTB), NT2AM(ISYRES) ) 328 KEND2 = KCTR2 + NT2SQ(ISYCTR) 329 LEND2 = LWORK - KEND2 330 331 IF (LEND2 .LT. NT2AM(ISYCTR) ) THEN 332 CALL QUIT('Insufficient work space in CC_FAMAT.') 333 END IF 334 335*---------------------------------------------------------------------* 336* read response amplitudes T2B and scale the diagonal: 337*---------------------------------------------------------------------* 338 IOPT = 2 339 CALL CC_RDRSP(LISTB,ITAMPB,ISYMTB,IOPT,MODEL, 340 & WORK(KDUM),WORK(KT2AMPB) ) 341 342 CAll CCLR_DIASCL(WORK(KT2AMPB),TWO,ISYMTB) 343 344 IF (LOCDBG) THEN 345 WRITE (LUPRI,*) MSGDBG, 'B response amplitudes:' 346 Call CC_PRP(WORK(KT1AMPB),WORK(KT2AMPB),ISYMTB,1,1) 347 END IF 348 349*---------------------------------------------------------------------* 350* read packed lagrangian multipliers and square them up: 351*---------------------------------------------------------------------* 352 IOPT = 2 353 Call CC_RDRSP(LISTC,IZETVC,ISYCTR,IOPT,MODEL, 354 & WORK(KDUM),WORK(KEND2)) 355 356 CALL CC_T2SQ (WORK(KEND2), WORK(KCTR2), ISYCTR) 357 358*---------------------------------------------------------------------* 359* calculate X^B and Y^B intermediates: 360*---------------------------------------------------------------------* 361 ISYMXY = MULD2H(ISYCTR,ISYMTB) 362 363 KXBMAT = KEND2 364 KYBMAT = KXBMAT + NMATIJ(ISYMXY) 365 KSCR = KYBMAT + NMATAB(ISYMXY) 366 KEND3 = KSCR + NT1AM(ISYRES) 367 LEND3 = LWORK - KEND3 368 369 If (LEND3 .LT. 0) THEN 370 CALL QUIT('Insufficient work space in CC_FAMAT.') 371 END IF 372 373* calculate X^C & Y^C intermediate: 374 Call CC_XI(WORK(KXBMAT),WORK(KCTR2), ISYCTR, 375 & WORK(KT2AMPB),ISYMTB,WORK(KEND3),LEND3) 376 377 Call CC_YI(WORK(KYBMAT),WORK(KCTR2), ISYCTR, 378 & WORK(KT2AMPB),ISYMTB,WORK(KEND3),LEND3) 379 380 IF (LOCDBG) THEN 381 WRITE (LUPRI,*) MSGDBG,'response X intermediate:' 382 WRITE (LUPRI,'(5f12.6)') (WORK(KXBMAT-1+I),I=1,NMATIJ(ISYMXY)) 383 WRITE (LUPRI,*) MSGDBG,'response Y intermediate:' 384 WRITE (LUPRI,'(2f12.6)') (WORK(KYBMAT-1+I),I=1,NMATAB(ISYMXY)) 385 END IF 386 387* calculate XY^B: XY^B_ij = X^B_ji, XY^B_bd = -Y^B_bd 388* 1.) transpose X^B intermediate 389 DO ISYMI = 1, NSYM 390 ISYMJ = MULD2H(ISYMI,ISYMXY) 391 IF (ISYMJ .LE. ISYMI) THEN 392 DO I = 1, NRHF(ISYMI) 393 MAXJ = NRHF(ISYMJ) 394 IF (ISYMJ .EQ. ISYMI) MAXJ = I-1 395 DO J = 1, MAXJ 396 NIJ = IMATIJ(ISYMI,ISYMJ) + NRHF(ISYMI)*(J-1) + I 397 NJI = IMATIJ(ISYMJ,ISYMI) + NRHF(ISYMJ)*(I-1) + J 398 SWAP = WORK(KXBMAT-1+NIJ) 399 WORK(KXBMAT-1+NIJ) = WORK(KXBMAT-1+NJI) 400 WORK(KXBMAT-1+NJI) = SWAP 401 END DO 402 END DO 403 END IF 404 END DO 405 406* 2.) multiply Y^B intermediate with -1: 407 Call DSCAL(NMATAB(ISYMXY), -ONE, WORK(KYBMAT), 1) 408 409 410* do one-index transformation of XY^B with A integrals: 411 IOPT = 2 412 Call CCG_1ITRVO(WORK(KSCR),ISYRES, 413 & WORK(KXBMAT),WORK(KYBMAT),ISYMXY, 414 & WORK(KPERTA),ISYMTA, IOPT ) 415 416* add contribution to GAMMA1: 417 Call DAXPY (NT1AM(ISYRES), ONE, WORK(KSCR),1, WORK(KGAMMA1), 1) 418 419 IF (LOCDBG) THEN 420 WRITE (LUPRI,*) MSGDBG,'A integrals:' 421 WRITE (LUPRI,'(5f12.6)') (WORK(KPERTA-1+I),I=1,NT1AM(ISYMTA)) 422 WRITE (LUPRI,*) MSGDBG, 'CC2/CCSD contribution to singles part:' 423 Call CC_PRP(WORK(KSCR),WORK(KDUM),ISYRES,1,0) 424 WRITE (LUPRI,*) MSGDBG, 'GAMMA1 now:' 425 Call CC_PRP(WORK(KGAMMA1),WORK(KDUM),ISYRES,1,0) 426 END IF 427 428*---------------------------------------------------------------------* 429 430*=====================================================================* 431* CC2/CCSD part for the doubles: <Zeta_2| [[A,T1B], tau_2] |HF> 432*=====================================================================* 433 434*---------------------------------------------------------------------* 435* reorganize work space, so that the result vector GAMMA2 can be 436* stored at the early beginning of the work space 437*---------------------------------------------------------------------* 438 KGAMMA1 = 1 439 KGAMMA2 = KGAMMA1 + NT1AM(ISYRES) 440 KEND0 = KGAMMA2 + NT2AM(ISYRES) 441 LEND0 = LWORK - KEND0 442 443 IF (KEND0 .GT. KCTR2) THEN 444 CALL QUIT('memory organization mixed up in CC_FAMAT.') 445 END IF 446 447 KEMAT1 = KCTR2 + NT2SQ(ISYCTR) 448 KEMAT2 = KEMAT1 + NMATAB(ISYTATB) 449 KEND3 = KEMAT2 + NMATIJ(ISYTATB) 450 LEND3 = LWORK - KEND3 451 452 IF ( LEND3 .LT. 0 ) THEN 453 CALL QUIT('Insufficient work space in CC_FAMAT.') 454 END IF 455 456*---------------------------------------------------------------------* 457* transpose tA^B(a b) --> EMAT1(b a) 458*---------------------------------------------------------------------* 459 DO ISYMA = 1, NSYM 460 ISYMB = MULD2H(ISYMA,ISYTATB) 461 DO A = 1, NVIR(ISYMA) 462 DO B = 1, NVIR(ISYMB) 463 NAB = IMATAB(ISYMA,ISYMB) + NVIR(ISYMA)*(B-1) + A 464 NBA = IMATAB(ISYMB,ISYMA) + NVIR(ISYMB)*(A-1) + B 465 466 WORK(KEMAT1 - 1 + NBA) = WORK(KBTAVV - 1 + NAB) 467 END DO 468 END DO 469 END DO 470 471 472*---------------------------------------------------------------------* 473* transpose tA^B(i j) --> EMAT2(j i) 474*---------------------------------------------------------------------* 475 DO ISYMI = 1, NSYM 476 ISYMJ = MULD2H(ISYMI,ISYTATB) 477 DO I = 1, NRHF(ISYMI) 478 DO J = 1, NRHF(ISYMJ) 479 NIJ = IMATIJ(ISYMI,ISYMJ) + NRHF(ISYMI)*(J-1) + I 480 NJI = IMATIJ(ISYMJ,ISYMI) + NRHF(ISYMJ)*(I-1) + J 481 482 WORK(KEMAT2 - 1 + NJI) = WORK(KBTAOO - 1 + NIJ) 483 END DO 484 END DO 485 END DO 486 487 IF (LOCDBG) THEN 488 WRITE (LUPRI,*) MSGDBG,'one-index trans. A (occ/occ block):' 489 WRITE (LUPRI,'(5f12.6)') (WORK(KBTAOO-1+I),I=1,NMATIJ(ISYTATB)) 490 WRITE (LUPRI,*) MSGDBG,'one-index trans. A (vir/vir block):' 491 WRITE (LUPRI,'(2f12.6)') (WORK(KBTAVV-1+I),I=1,NMATAB(ISYTATB)) 492 WRITE (LUPRI,*) MSGDBG,'EMAT2:' 493 WRITE (LUPRI,'(5f12.6)') (WORK(KEMAT2-1+I),I=1,NMATIJ(ISYTATB)) 494 WRITE (LUPRI,*) MSGDBG,'EMAT1:' 495 WRITE (LUPRI,'(2f12.6)') (WORK(KEMAT1-1+I),I=1,NMATAB(ISYTATB)) 496 END IF 497 498*---------------------------------------------------------------------* 499* combine EMAT1/EMAT2 with lagrangian multipliers: 500* (note: this overwrites the intermedites stored at the beginning 501* of the work space...) 502*---------------------------------------------------------------------* 503* initialize GAMMA2: 504 CALL DZERO(WORK(KGAMMA2),NT2AM(ISYRES)) 505 506* do the caculation: 507 Call CCRHS_E(WORK(KGAMMA2),WORK(KCTR2),WORK(KEMAT1), 508 & WORK(KEMAT2), WORK(KEND3), LEND3, ISYCTR, ISYTATB) 509 510*---------------------------------------------------------------------* 511 IF (LOCDBG) THEN 512 WRITE (LUPRI,*) MSGDBG, 'GAMMA:' 513 Call CC_PRP(WORK(KGAMMA1),WORK(KGAMMA2),ISYRES,1,1) 514 END IF 515 516*=====================================================================* 517* CC-R12 part: 518*=====================================================================* 519 IF (CCR12) THEN 520 IF (.NOT.IANR12.EQ.1) THEN 521 CALL QUIT('Sorry, only CC-R12 Ansatz 1 implemented yet!') 522 ENDIF 523* allocate memory: 524 KGAMMAR12 = KEND0 525 KEND0 = KGAMMAR12 + NTR12AM(ISYRES) 526 LEND0 = LWORK - KEND0 527C 528 KCTRR12 = KEND0 529 KT12AMB = KCTRR12 + NTR12SQ(ISYCTR) 530 KXINTTRI= KT12AMB + MAX(NTR12SQ(ISYMTB),NTR12SQ(ISYRES)) 531 KXINTSQ = KXINTTRI+ NR12R12P(1) 532 KSCR = KXINTSQ + NR12R12SQ(1) 533 KCTMO = KSCR + NTR12SQ(ISYRES) 534 KT1AMP0 = KCTMO + N2BST(ISYMTA) 535 KLAMDP0 = KT1AMP0 + NT1AMX 536 KLAMDH0 = KLAMDP0 + NLAMDT 537 KT1AMPB = KLAMDH0 + NLAMDT 538 KLAMDPB = KT1AMPB + NT1AM(ISYMTB) 539 KLAMDHB = KLAMDPB + NGLMDT(ISYMTB) 540 KEND1 = KLAMDHB + NGLMDT(ISYMTB) 541 LEND1 = LWORK - KEND1 542 IF (LEND1.LT.0) THEN 543 CALL QUIT('Insufficient work space for R12 in CC_FAMAT') 544 END IF 545 546*---------------------------------------------------------------------* 547* read in everything needed: 548*---------------------------------------------------------------------* 549* read the AO integrals: 550 CALL CCPRPAO(LABELA,.TRUE.,WORK(KCTMO),IRREP,ISYM,IERR, 551 & WORK(KEND1),LEND1) 552 IF (IERR.NE.0 .OR. IRREP.NE.ISYMTA) THEN 553 CALL QUIT('CC_FAMAT: error while reading operator '//LABELA) 554 END IF 555 556* get MO coefficients: 557 CALL DZERO(WORK(KT1AMP0),NT1AMX) 558 CALL LAMMAT(WORK(KLAMDP0),WORK(KLAMDH0),WORK(KT1AMP0), 559 & WORK(KEND1),LEND1) 560* transform MO matrices with C1: 561 IOPT = 1 562 Call CC_RDRSP(LISTB,ITAMPB,ISYMTB,IOPT,MODEL,WORK(KT1AMPB), 563 & DUMMY) 564 CALL CCLR_LAMTRA(WORK(KLAMDP0),WORK(KLAMDPB),WORK(KLAMDH0), 565 & WORK(KLAMDHB),WORK(KT1AMPB),ISYMTB) 566 567* read R12 Lagrangian multipliers: 568 CALL CC_R12GETCT(WORK(KCTRR12),ISYCTR,2,2.0D0*BRASCL,.FALSE., 569 & 'N',IDUM,IDUM,IDUM,LISTC,IZETVC, 570 & WORK(KEND1),LEND1) 571 572* read R12 response amplitudes: 573 CALL CC_R12GETCT(WORK(KT12AMB),ISYMTB,2,KETSCL,.FALSE.,'N', 574 & IDUM,IDUM,IDUM,LISTB,ITAMPB,WORK(KEND1),LEND1) 575 576* read R12 overlap matrix: 577 LUNIT = -1 578 CALL GPOPEN(LUNIT,FCCR12X,'OLD',' ','UNFORMATTED', 579 & IDUM,.FALSE.) 580 REWIND(LUNIT) 581 8888 READ(LUNIT) IAN 582 READ(LUNIT) (WORK(KXINTTRI-1+I), I=1, NR12R12P(1)) 583 IF (IAN.NE.IANR12) GOTO 8888 584 CALL GPCLOSE(LUNIT,'KEEP') 585 IOPT = 2 586 CALL CCR12UNPCK2(WORK(KXINTTRI),1,WORK(KXINTSQ),'N',IOPT) 587 588*---------------------------------------------------------------------* 589* Calculate singles contribution: 590*---------------------------------------------------------------------* 591 CALL CC_R12ETAA(WORK(KGAMMA1),ISYRES,WORK(KCTRR12),ISYCTR, 592 & WORK(KT12AMB),ISYMTB,WORK(KXINTSQ),WORK(KCTMO), 593 & ISYMTA,WORK(KLAMDP0),WORK(KLAMDH0),.FALSE., 594 & WORK(KEND1),LEND1) 595 596*---------------------------------------------------------------------* 597* Calculate R12 doubles contribution: 598*---------------------------------------------------------------------* 599 CALL DZERO(WORK(KT12AMB),NTR12SQ(ISYRES)) 600 CALL CC_R12XI2A(WORK(KSCR),ISYRES,WORK(KCTRR12),ISYCTR, 601 & WORK(KCTMO),ISYMTA,WORK(KLAMDP0),WORK(KLAMDHB), 602 & ISYMTB,'T',WORK(KEND1),LEND1) 603 CALL CC_R12XI2B(WORK(KT12AMB),'N',WORK(KXINTSQ),1,'N', 604 & WORK(KSCR),ISYRES,'N',-1.0D0) 605 IOPT = 1 606 CALL CCR12PCK2(WORK(KGAMMAR12),ISYRES,.FALSE.,WORK(KT12AMB), 607 & 'N',IOPT) 608 CALL CCLR_DIASCLR12(WORK(KGAMMAR12),0.5D0*KETSCL,ISYRES) 609 610 END IF 611 612*=====================================================================* 613* CC3 part: 614*=====================================================================* 615 IF (CCSDT) THEN 616 IF (IOPTRES.LT.5) THEN 617 KGAMMA1EFF = KEND0 618 KGAMMA2EFF = KGAMMA1EFF + NT1AM(ISYRES) 619 KEND0 = KGAMMA2EFF + NT2AM(ISYRES) 620 END IF 621 LEND0 = LWORK - KEND0 622 623 IF ( LEND0 .LT. 0 ) THEN 624 CALL QUIT('Insufficient work space in CC_FAMAT.') 625 END IF 626 627 IF (IOPTRES.GE.1 .AND. IOPTRES.LE.4) THEN 628 CALL QUIT('CC_FAMAT needs to be fixed for this case') 629 ! in this case we need to find the frequency associated 630 ! with the A perturbation, such that we can construct 631 ! the correct effective rhs vector... 632 ! FREQA = 633 ELSE IF (IOPTRES.EQ.5) THEN 634 FREQA = 0.0D0 635 ELSE 636 CALL QUIT('Illegal value for IOPTRES in CC_FAMAT.') 637 END IF 638 639 IF (IOPTRES.NE.5 .OR. NODDY_CCSDT) THEN 640 641 CALL CCSDT_FAMAT_NODDY(LISTC,IZETVC,LISTB,ITAMPB,IOPTRES, 642 & LABELA,FREQA, 643 & WORK(KGAMMA1),WORK(KGAMMA2), 644 & WORK(KGAMMA1EFF),WORK(KGAMMA2EFF), 645 & IFADOTS,FACON,FILFA,ITRAN, 646 & NFATRAN,MXVEC,WORK(KEND0),LEND0 ) 647 648 END IF 649 650 END IF 651 652*=====================================================================* 653* final section: depending on IOPTRES store vector in memory or on 654* file or contract it with some other vectors: 655* 656* the memory has to be organized as follows: 657* kgamma1 singles result vector 658* kgamma2 doubles result vector 659* kgammar12 r12 doubles result vector 660* kgamma1eff effective singles result vector for CC3 661* kgamma2eff effective doubles result vector for CC3 662* kend0 start of unused work space 663*=====================================================================* 6649999 CONTINUE 665 666 LEND0 = LWORK - KEND0 667 668 IF (CCS) THEN 669 MODELW = 'CCS ' 670 IOPTW = 1 671 ELSE IF (CC2) THEN 672 MODELW = 'CC2 ' 673 IOPTW = 3 674 ELSE IF (CCSD) THEN 675 MODELW = 'CCSD ' 676 IOPTW = 3 677 ELSE IF (CC3) THEN 678 MODELW = 'CC3 ' 679 IOPTW = 3 680 IOPTWE = 24 681 ELSE 682 CALL QUIT('Unknown coupled cluster model in CC_FAMAT.') 683 END IF 684 IF (CCR12) THEN 685 APROXR12 = ' ' 686 CALL CCSD_MODEL(MODELW,LENMOD,10,MODELW,10,APROXR12) 687 IOPTWR12 = 32 688 END IF 689 690 691 IF (IOPTRES.EQ.0 .OR. IOPTRES.EQ.1) THEN 692 CALL QUIT('IOPTRES=0,1 not implemented in CC_FAMAT.') 693 694 ELSE IF (IOPTRES.EQ.3) THEN 695 IFILE = IFATRAN(4,ITRAN) 696 IF (ILSTSYM(FILFA,IFILE).NE.ISYRES) THEN 697 CALL QUIT('Symmetry mismatch for result vector in CC_FAMAT.') 698 END IF 699 CALL CC_WRRSP(FILFA,IFILE,ISYRES,IOPTW,MODELW,DUMMY, 700 & WORK(KGAMMA1),WORK(KGAMMA2),WORK(KEND0),LEND0) 701 IF (CCR12) THEN 702 CALL CC_WRRSP(FILFA,IFILE,ISYRES,IOPTWR12, 703 & MODELW,DUMMY,DUMMY,WORK(KGAMMAR12), 704 & WORK(KEND0),LEND0) 705 END IF 706 IF (CCSDT) THEN 707 CALL CC_WRRSP(FILFA,IFILE,ISYRES,IOPTWE,MODELW,DUMMY, 708 & WORK(KGAMMA1EFF),WORK(KGAMMA2EFF),WORK(KEND0),LEND0) 709 END IF 710 ELSE IF (IOPTRES.EQ.4) THEN 711 IFILE = IFATRAN(4,ITRAN) 712 IF (ILSTSYM(FILFA,IFILE).NE.ISYRES) THEN 713 CALL QUIT('Symmetry mismatch for result vector in CC_FAMAT.') 714 END IF 715 CALL CC_WARSP(FILFA,IFILE,ISYRES,IOPTW,MODELW,DUMMY, 716 & WORK(KGAMMA1),WORK(KGAMMA2),WORK(KEND0),LEND0) 717 IF (CCR12) THEN 718 CALL CC_WARSP(FILFA,IFILE,ISYRES,IOPTWR12, 719 & MODELW,DUMMY,DUMMY,WORK(KGAMMAR12), 720 & WORK(KEND0),LEND0) 721 END IF 722 IF (CCSDT) THEN 723 CALL CC_WARSP(FILFA,IFILE,ISYRES,IOPTWE,MODELW,DUMMY, 724 & WORK(KGAMMA1EFF),WORK(KGAMMA2EFF),WORK(KEND0),LEND0) 725 END IF 726 ELSE IF (IOPTRES.EQ.5) THEN 727 IF (LOCDBG) THEN 728 IVEC = 1 729 WRITE(LUPRI,*) 'FACON TRIPLES CONTRIBUTION:' 730 DO WHILE (IFADOTS(IVEC,ITRAN).NE.0 .AND. IVEC.LE.MXVEC) 731 WRITE (LUPRI,*) 732 & 'FACON:',IVEC,ITRAN,FACON(IVEC,ITRAN),IOPTW 733 IVEC = IVEC + 1 734 END DO 735 END IF 736 737 IF (.NOT.CCS) CALL CCLR_DIASCL(WORK(KGAMMA2),TWO,ISYRES) 738 CALL CCDOTRSP(IFADOTS,FACON,IOPTW,FILFA,ITRAN,NFATRAN,MXVEC, 739 & WORK(KGAMMA1),WORK(KGAMMA2),ISYRES, 740 & WORK(KEND0),LEND0) 741 742 IF (CCR12) THEN 743 CALL CCDOTRSP(IFADOTS,FACON,IOPTWR12,FILFA,ITRAN,NFATRAN,MXVEC, 744 & DUMMY,WORK(KGAMMAR12),ISYRES, 745 & WORK(KEND0),LEND0) 746 END IF 747 748 IF (LOCDBG) THEN 749 IVEC = 1 750 DO WHILE (IFADOTS(IVEC,ITRAN).NE.0 .AND. IVEC.LE.MXVEC) 751 WRITE (LUPRI,*) 752 & 'FACON:',IVEC,ITRAN,FACON(IVEC,ITRAN),IOPTW 753 IVEC = IVEC + 1 754 END DO 755 END IF 756 ELSE 757 CALL QUIT('Illegal value for IOPTRES in CC_FAMAT.') 758 END IF 759 760 CALL QEXIT('CC_FAMAT') 761 RETURN 762 END 763*=====================================================================* 764* END OF SUBROUTINE CC_FAMAT * 765*=====================================================================* 766*=====================================================================* 767 SUBROUTINE CCSDT_FAMAT_NODDY(LISTL,IDLSTL,LISTB,IDLSTB,IOPTRES, 768 & LABELA,FREQA, 769 & OMEGA1,OMEGA2, 770 & OMEGA1EFF,OMEGA2EFF, 771 & IDOTS,DOTPROD,LISTDP,ITRAN, 772 & NFATRAN,MXVEC,WORK,LWORK ) 773*---------------------------------------------------------------------* 774* 775* Purpose: compute triples contribution to F{A} transformed vector 776* 777* (F{A} T^B)^eff_1,2 = (F{A} T^B)_1,2(CCSD) + (F{A} T^B)_1,2(T^B_3) 778* - (F{A} T^B)_3 A_3;1,2 (w_3 - w)^1 779* 780* 781* Written by Christof Haettig, April 2002 782* based on CCSDT_ETA_NODDY 783* Extensions for cubic response, CCH, May 2003 784* 785*=====================================================================* 786 IMPLICIT NONE 787#include "priunit.h" 788#include "ccsdinp.h" 789#include "maxorb.h" 790#include "ccsdsym.h" 791#include "ccfield.h" 792#include "ccorb.h" 793#include "dummy.h" 794 795 LOGICAL LOCDBG 796 PARAMETER (LOCDBG=.FALSE.) 797 798 INTEGER ISYM0 799 PARAMETER (ISYM0 = 1) 800 801 CHARACTER LABELA*8 802 CHARACTER*3 LISTL, LISTDP, LISTB 803 INTEGER LWORK, IDLSTL, IDLSTB, IOPTRES, ITRAN, MXVEC, NFATRAN 804 INTEGER IDOTS(MXVEC,NFATRAN) 805 806#if defined (SYS_CRAY) 807 REAL DOTPROD(MXVEC,NFATRAN), DDOT 808 REAL WORK(LWORK), FREQL, FREQA, FREQB, FREQF, FREQC 809 REAL OMEGA1(*), OMEGA2(*) 810 REAL OMEGA1EFF(*), OMEGA2EFF(*) 811 REAL SIXTH, ONE, TWO, TCON, DCON, SCON, FF, SIGN 812#else 813 DOUBLE PRECISION DOTPROD(MXVEC,NFATRAN), DDOT 814 DOUBLE PRECISION WORK(LWORK), FREQL, FREQA, FREQB, FREQE, FREQC 815 DOUBLE PRECISION OMEGA1(*), OMEGA2(*) 816 DOUBLE PRECISION OMEGA1EFF(*), OMEGA2EFF(*) 817 DOUBLE PRECISION SIXTH, ONE, TWO, TCON, DCON, SCON, FF, SIGN 818#endif 819 PARAMETER(SIXTH=1.0D0/6.0D0, ONE=1.0D0, TWO=2.0D0) 820 821 CHARACTER*10 MODEL 822 LOGICAL L2INCL 823 INTEGER INDEX, LUSIFC, IOPT, ISYMD, ILLL, IDEL, ISYDIS, NIJ, IJ, 824 & IVEC, IDLSTC, ISYMC, LUFOCK, ILSTSYM, ISYML 825 INTEGER KSCR1, KFOCKD, KEND1, KT1AMP0, KLAMP0, KLAMH0, 826 & KINT1T0, KINT2T0, KINT1S0, KINT2S0, KXIAJB, KYIAJB, 827 & K0IOVVO, K0IOOVV, K0IOOOO, K0IVVVV, KOME1, KOME2, KDUM, 828 & KXINT, KEND2, LWRK2, KL1AM, KL2AM, KL3AM, KT3AM, KT2AM, 829 & KEND3, LWRK3, KINT1SC, KINT2SC, KLAMPC, KLAMHC,KFOCKAB, 830 & KFOCKC, LWRK1, KE3AM, KTC3AM, KTC1AM, KTC2AM, 831 & ISYMA, KINT1SB, KINT2SB, KLAMPB, KLAMHB, KFOCKB, ISYMB, 832 & KFOCKA, KFOCKA_AO, KFOCK0, IRREP, ISYM, IERR, KFOCKAB1, 833 & KFIELD, KFIELDAO, KT1AMB 834 835 INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 836 837 CALL QENTER('CCSDT_FAMAT_NODDY') 838 839 KDUM = 1 840 IF (DIRECT) CALL QUIT('CCSDT_FAMAT_NODDY: DIRECT NOT IMPLEMENTED') 841 842 IF (LOCDBG) THEN 843 WRITE(LUPRI,*) 'CCSDT_FAMAT_NODDY> left vector:',LISTL,IDLSTL 844 WRITE(LUPRI,*) 'CCSDT_FAMAT_NODDY> right vector:',LISTB,IDLSTB 845 WRITE(LUPRI,*) 'CCSDT_FAMAT_NODDY> operator :',LABELA 846 WRITE(LUPRI,*) 'CCSDT_FAMAT_NODDY> frequency :',FREQA 847 WRITE(LUPRI,*) 'CCSDT_FAMAT_NODDY> result vector on entry:' 848 CALL CC_PRP(OMEGA1,OMEGA2,1,1,1) 849 END IF 850 851*---------------------------------------------------------------------* 852* Memory allocation: 853*---------------------------------------------------------------------* 854 KSCR1 = 1 855 KFOCKD = KSCR1 + NT1AMX 856 KFOCK0 = KFOCKD + NORBT 857 KEND1 = KFOCK0 + NORBT*NORBT 858 859 KFOCKA = KEND1 860 KFOCKA_AO = KFOCKA + NORBT*NORBT 861 KEND1 = KFOCKA_AO + NORBT*NORBT 862 863 IF (NONHF) THEN 864 KFIELD = KEND1 865 KFIELDAO = KFIELD + NORBT*NORBT 866 KEND1 = KFIELDAO + NORBT*NORBT 867 END IF 868 869 KT1AMP0 = KEND1 870 KLAMP0 = KT1AMP0 + NT1AMX 871 KLAMH0 = KLAMP0 + NLAMDT 872 KEND1 = KLAMH0 + NLAMDT 873 874 KINT1T0 = KEND1 875 KINT2T0 = KINT1T0 + NT1AMX*NVIRT*NVIRT 876 KEND1 = KINT2T0 + NRHFT*NRHFT*NT1AMX 877 878 KINT1S0 = KEND1 879 KINT2S0 = KINT1S0 + NT1AMX*NVIRT*NVIRT 880 KEND1 = KINT2S0 + NRHFT*NRHFT*NT1AMX 881 882 KXIAJB = KEND1 883 KYIAJB = KXIAJB + NT1AMX*NT1AMX 884 KEND1 = KYIAJB + NT1AMX*NT1AMX 885 886 K0IOVVO = KEND1 887 K0IOOVV = K0IOVVO + NRHFT*NVIRT*NVIRT*NRHFT 888 K0IOOOO = K0IOOVV + NRHFT*NVIRT*NVIRT*NRHFT 889 K0IVVVV = K0IOOOO + NRHFT*NRHFT*NRHFT*NRHFT 890 KEND1 = K0IVVVV + NVIRT*NVIRT*NVIRT*NVIRT 891 892 KOME1 = KEND1 893 KOME2 = KOME1 + NT1AMX 894 KEND1 = KOME2 + NT1AMX*NT1AMX 895 896 KFOCKAB1= KEND1 897 KFOCKAB = KFOCKAB1+ NORBT*NORBT 898 KEND1 = KFOCKAB + NORBT*NORBT 899 900 LWRK1 = LWORK - KEND1 901 IF (LWRK1 .LT. 0) THEN 902 CALL QUIT('Insufficient space in CCSDT_FAMAT_NODDY') 903 ENDIF 904 905*---------------------------------------------------------------------* 906* Read SCF orbital energies from file: 907*---------------------------------------------------------------------* 908 LUSIFC = -1 909 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY, 910 & .FALSE.) 911 REWIND LUSIFC 912 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 913 READ (LUSIFC) 914 READ (LUSIFC) (WORK(KFOCKD+I-1), I=1,NORBT) 915 CALL GPCLOSE(LUSIFC,'KEEP') 916 917*---------------------------------------------------------------------* 918* Get zeroth-order Lambda matrices: 919*---------------------------------------------------------------------* 920 IOPT = 1 921 Call CC_RDRSP('R0',0,ISYM0,IOPT,MODEL,WORK(KT1AMP0),WORK(KDUM)) 922 923 Call LAMMAT(WORK(KLAMP0),WORK(KLAMH0),WORK(KT1AMP0), 924 & WORK(KEND1),LWRK1) 925 926*---------------------------------------------------------------------* 927* read zeroth-order AO Fock matrix from file: 928*---------------------------------------------------------------------* 929 LUFOCK = -1 930 CALL GPOPEN(LUFOCK,'CC_FCKH','OLD',' ','UNFORMATTED',IDUMMY, 931 & .FALSE.) 932 REWIND(LUFOCK) 933 READ(LUFOCK) (WORK(KFOCK0-1+I),I=1,N2BST(ISYM0)) 934 CALL GPCLOSE(LUFOCK,'KEEP') 935 936 CALL CC_FCKMO(WORK(KFOCK0),WORK(KLAMP0),WORK(KLAMH0), 937 & WORK(KEND1),LWRK1,ISYM0,ISYM0,ISYM0) 938 939*---------------------------------------------------------------------* 940* If needed get external field: 941*---------------------------------------------------------------------* 942 IF (NONHF) THEN 943 CALL DZERO(WORK(KFIELDAO),NORBT*NORBT) 944 DO I = 1, NFIELD 945 FF = EFIELD(I) 946 CALL CC_ONEP(WORK(KFIELDAO),WORK(KEND1),LWRK1,FF,1,LFIELD(I)) 947 ENDDO 948 CALL DCOPY(NORBT*NORBT,WORK(KFIELDAO),1,WORK(KFIELD),1) 949 950 ! calculate external field in zero-order lambda basis 951 CALL CC_FCKMO(WORK(KFIELD),WORK(KLAMP0),WORK(KLAMH0), 952 * WORK(KEND1),LWRK1,1,1,1) 953 954 IF (LOCDBG) WRITE(LUPRI,*) 'NORM^2(FIELD):', 955 & DDOT(NORBT*NORBT,WORK(KFIELD),1,WORK(KFIELD),1) 956 ENDIF 957 958*---------------------------------------------------------------------* 959* Get property integrals and transform them to the MO basis: 960*---------------------------------------------------------------------* 961 ISYMA = 1 ! since this code is limited to C1 symmetry... 962 963 CALL CCPRPAO(LABELA,.TRUE.,WORK(KFOCKA_AO),IRREP,ISYM,IERR, 964 & WORK(KEND1),LWRK1) 965 IF (IERR.NE.0 .OR. IRREP.NE.ISYMA) THEN 966 CALL QUIT('CCSDT_FA_NODDY: error reading operator '//LABELA) 967 END IF 968 969 CALL DCOPY(NORBT*NORBT,WORK(KFOCKA_AO),1,WORK(KFOCKA),1) 970 971 CALL CC_FCKMO(WORK(KFOCKA),WORK(KLAMP0),WORK(KLAMH0), 972 & WORK(KEND1),LWRK1,ISYMA,ISYM0,ISYM0) 973 974*---------------------------------------------------------------------* 975* Compute some integral intermediates: 976*---------------------------------------------------------------------* 977 978 CALL DZERO(WORK(KINT1T0),NT1AMX*NVIRT*NVIRT) 979 CALL DZERO(WORK(KINT2T0),NRHFT*NRHFT*NT1AMX) 980 981 CALL DZERO(WORK(KINT1S0),NT1AMX*NVIRT*NVIRT) 982 CALL DZERO(WORK(KINT2S0),NRHFT*NRHFT*NT1AMX) 983 984 CALL DZERO(WORK(KXIAJB), NT1AMX*NT1AMX) 985 CALL DZERO(WORK(KYIAJB), NT1AMX*NT1AMX) 986 987 CALL DZERO(WORK(K0IOVVO),NRHFT*NVIRT*NVIRT*NRHFT) 988 CALL DZERO(WORK(K0IOOVV),NRHFT*NVIRT*NVIRT*NRHFT) 989 CALL DZERO(WORK(K0IOOOO),NRHFT*NRHFT*NRHFT*NRHFT) 990 CALL DZERO(WORK(K0IVVVV),NVIRT*NVIRT*NVIRT*NVIRT ) 991 992 993 DO ISYMD = 1, NSYM 994 DO ILLL = 1,NBAS(ISYMD) 995 IDEL = IBAS(ISYMD) + ILLL 996 ISYDIS = MULD2H(ISYMD,ISYMOP) 997 998C ---------------------------- 999C Work space allocation no. 2. 1000C ---------------------------- 1001 KXINT = KEND1 1002 KEND2 = KXINT + NDISAO(ISYDIS) 1003 LWRK2 = LWORK - KEND2 1004 IF (LWRK2 .LT. 0) THEN 1005 WRITE(LUPRI,*) 'Need : ',KEND2,'Available : ',LWORK 1006 CALL QUIT('Insufficient space in CCSDT_FAMAT_NODDY') 1007 ENDIF 1008 1009C --------------------------- 1010C Read in batch of integrals. 1011C --------------------------- 1012 CALL CCRDAO(WORK(KXINT),IDEL,1,WORK(KEND2),LWRK2, 1013 * WORK(KDUM),DIRECT) 1014 1015C ---------------------------------- 1016C Calculate integrals needed in CC3: 1017C ---------------------------------- 1018 CALL CCSDT_TRAN1(WORK(KINT1T0),WORK(KINT2T0), 1019 & WORK(KLAMP0),WORK(KLAMH0), 1020 & WORK(KXINT),IDEL) 1021 1022 CALL CC3_TRAN2(WORK(KXIAJB),WORK(KYIAJB), 1023 & WORK(KLAMP0),WORK(KLAMH0), 1024 & WORK(KXINT),IDEL) 1025 1026 CALL CCSDT_TRAN3(WORK(KINT1S0),WORK(KINT2S0),WORK(KLAMP0), 1027 & WORK(KLAMH0),WORK(KXINT),IDEL) 1028 1029 CALL CCFOP_TRAN1_R(WORK(K0IOVVO),WORK(K0IOOVV), 1030 & WORK(K0IOOOO),WORK(K0IVVVV), 1031 & WORK(KLAMP0),WORK(KLAMH0), 1032 & WORK(KLAMP0),WORK(KLAMH0), 1033 & WORK(KLAMP0),WORK(KLAMH0), 1034 & WORK(KXINT),IDEL) 1035 1036 END DO 1037 END DO 1038 1039*---------------------------------------------------------------------* 1040* Some more memory allocations: 1041*---------------------------------------------------------------------* 1042 KL1AM = KEND1 1043 KL2AM = KL1AM + NT1AMX 1044 KL3AM = KL2AM + NT1AMX*NT1AMX 1045 KEND2 = KL3AM + NT1AMX*NT1AMX*NT1AMX 1046 LWRK2 = LWORK - KEND2 1047 1048 KT3AM = KEND2 1049 KT2AM = KT3AM + NT1AMX*NT1AMX*NT1AMX 1050 KEND3 = KT2AM + NT1AMX*NT1AMX 1051 LWRK3 = LWORK - KEND3 1052 1053 IF (LWRK3 .LT. NT2AMX) THEN 1054 CALL QUIT('Insufficient space in CCSDT_FAMAT_NODDY') 1055 ENDIF 1056 1057 ! read T^B doubles amplitudes from file and square up 1058 ISYMB = ILSTSYM(LISTB,IDLSTB) 1059 IOPT = 2 1060 Call CC_RDRSP(LISTB,IDLSTB,ISYMB,IOPT,MODEL, 1061 & WORK(KDUM),WORK(KEND3)) 1062 Call CCLR_DIASCL(WORK(KEND3),TWO,ISYMB) 1063 CALL CC_T2SQ(WORK(KEND3),WORK(KT2AM),ISYMB) 1064 1065 ! read L^0 multipliers from file and square up doubles part 1066 ISYML = ILSTSYM(LISTL,IDLSTL) 1067 IOPT = 3 1068 Call CC_RDRSP(LISTL,IDLSTL,ISYML,IOPT,MODEL, 1069 & WORK(KL1AM),WORK(KT3AM)) 1070 CALL CC_T2SQ(WORK(KT3AM),WORK(KL2AM),ISYM0) 1071 1072*---------------------------------------------------------------------* 1073* Compute triples amplitude response: 1074*---------------------------------------------------------------------* 1075 KINT1SB = KEND3 1076 KINT2SB = KINT1SB + NT1AMX*NVIRT*NVIRT 1077 KEND3 = KINT2SB + NRHFT*NRHFT*NT1AMX 1078 1079 KT1AMB = KEND3 1080 KLAMPB = KT1AMB + NT1AMX 1081 KLAMHB = KLAMPB + NLAMDT 1082 KFOCKB = KLAMHB + NLAMDT 1083 KEND3 = KFOCKB + NORBT*NORBT 1084 1085 LWRK3 = LWORK - KEND3 1086 IF (LWRK2 .LT. 0) THEN 1087 CALL QUIT('Insufficient space in CCSDT_FAMAT_NODDY') 1088 ENDIF 1089 1090 IF (LISTB(1:3).EQ.'R1 ' .OR. LISTB(1:3).EQ.'RE ' .OR. 1091 & LISTB(1:3).EQ.'RC ' ) THEN 1092 1093 CALL CCSDT_T31_NODDY(WORK(KT3AM),LISTB,IDLSTB,FREQB,.FALSE., 1094 & .FALSE.,WORK(KINT1S0),WORK(KINT2S0), 1095 & .FALSE.,WORK(KDUM),WORK(KDUM), 1096 & .FALSE.,WORK(KDUM),WORK(KDUM), 1097 & WORK(KINT1SB),WORK(KINT2SB), 1098 & WORK(KLAMPB),WORK(KLAMHB),WORK(KFOCKB), 1099 & WORK(KLAMP0),WORK(KLAMH0),WORK(KFOCK0), 1100 & WORK(KDUM),WORK(KFOCKD), 1101 & WORK(KEND3),LWRK3) 1102 1103 ELSE IF (LISTB(1:3).EQ.'R2 ' .OR. LISTB(1:3).EQ.'ER1') THEN 1104 1105 CALL CCSDT_T32_NODDY(WORK(KT3AM),LISTB,IDLSTB,FREQB, 1106 & WORK(KINT1S0),WORK(KINT2S0), 1107 & WORK(KLAMP0),WORK(KLAMH0),WORK(KFOCK0), 1108 & WORK(KFOCKD),WORK(KFIELDAO),WORK(KFIELD), 1109 & WORK(KSCR1),WORK(KEND3),LWRK3) 1110 1111 IOPT = 1 1112 Call CC_RDRSP(LISTB,IDLSTB,ISYMB,IOPT,MODEL,WORK(KT1AMB),DUMMY) 1113 1114 CALL CCLR_LAMTRA(WORK(KLAMP0),WORK(KLAMPB),WORK(KLAMH0), 1115 & WORK(KLAMHB),WORK(KT1AMB),ISYMB) 1116 ELSE 1117 1118 CALL QUIT('Unknown list '//LISTB//' in CCSDT_FAMAT_NODDY.') 1119 1120 END IF 1121 1122*---------------------------------------------------------------------* 1123* Compute triples multipliers L_3: 1124*---------------------------------------------------------------------* 1125 IF (LISTL(1:3).EQ.'L0 ') THEN 1126 1127 FREQL = 0.0D0 1128 1129 CALL DZERO(WORK(KL3AM),NT1AMX*NT1AMX*NT1AMX) 1130 1131 IF (NONHF .AND. LWRK3.LT.NT1AMX*NT1AMX*NT1AMX) 1132 * CALL QUIT('Out of memory in CCSDT_FAMAT_NODDY.') 1133 1134 ! remember that CCSDT_L03AM returns -L3 !! 1135 CALL CCSDT_L03AM(WORK(KL3AM),WORK(KINT1T0),WORK(KINT2T0), 1136 * WORK(KXIAJB),WORK(KFOCK0),WORK(KL1AM), 1137 * WORK(KL2AM),WORK(KSCR1),WORK(KFOCKD), 1138 * WORK(KFIELD),WORK(KEND3)) 1139 1140 CALL DSCAL(NT1AMX*NT1AMX*NT1AMX,-1.0D0,WORK(KL3AM),1) 1141 1142 ELSE IF (LISTL(1:3).EQ.'L1 ' .OR. LISTL(1:3).EQ.'LE ' .OR. 1143 & LISTL(1:3).EQ.'M1 ' .OR. LISTL(1:3).EQ.'N2 ' .OR. 1144 & LISTL(1:3).EQ.'E0 ' ) THEN 1145 1146 CALL DZERO(WORK(KL3AM),NT1AMX*NT1AMX*NT1AMX) 1147 1148 CALL CCSDT_TBAR31_NODDY(WORK(KL3AM),FREQL,LISTL,IDLSTL, 1149 & WORK(KLAMP0),WORK(KLAMH0), 1150 & WORK(KFOCK0),WORK(KFOCKD),WORK(KSCR1), 1151 & WORK(KXIAJB),WORK(KINT1T0),WORK(KINT2T0), 1152 & WORK(KEND3),LWRK3) 1153 1154 CALL DSCAL(NT1AMX*NT1AMX*NT1AMX,-1.0D0,WORK(KL3AM),1) 1155 1156 ELSE 1157 1158 ! FREQL = ?? 1159 1160 CALL QUIT('CCSDT_FAMAT_NODDY> LISTL NOT AVAILABLE:'//LISTL) 1161 1162 END IF 1163 1164*---------------------------------------------------------------------* 1165* Compute contribution from <L_3|[[A,T^B_3],\tau_nu_1|HF>: 1166*---------------------------------------------------------------------* 1167 CALL DZERO(WORK(KOME1),NT1AMX) 1168 1169 CALL CCSDT_E1AM(WORK(KOME1),WORK(KL3AM),WORK(KT3AM),WORK(KFOCKA)) 1170 1171 DO I = 1,NT1AMX 1172 OMEGA1(I) = OMEGA1(I) + WORK(KOME1+I-1) 1173 END DO 1174 1175ctest 1176C LUFOCK = -1 1177C CALL GPOPEN(LUFOCK,'CCTEST'//LISTL(1:2)//LISTB(1:2), 1178C & 'UNKNOWN',' ','FORMATTED',IDUMMY,.FALSE.) 1179C REWIND(LUFOCK) 1180C WRITE(LUFOCK,'(A)') 'OMEGA1:' 1181C WRITE(LUFOCK,'(F20.16)') (WORK(KOME1-1+I),I=1,NT1AMX) 1182C WRITE(LUFOCK,'(A)') 'FOCKA:' 1183C WRITE(LUFOCK,'(F20.16)') (WORK(KFOCKA-1+I),I=1,NORBT*NORBT) 1184C WRITE(LUFOCK,'(A)') 'L3AM:' 1185C WRITE(LUFOCK,'(F20.16)')(WORK(KL3AM-1+I),I=1,NT1AMX*NT1AMX*NT1AMX) 1186C WRITE(LUFOCK,'(A)') 'T3AM:' 1187C WRITE(LUFOCK,'(F20.16)')(WORK(KT3AM-1+I),I=1,NT1AMX*NT1AMX*NT1AMX) 1188C CALL GPCLOSE(LUFOCK,'KEEP') 1189ctest 1190 IF (LOCDBG) THEN 1191 WRITE(LUPRI,*) 'CCSDT_FAMAT_NODDY> Contribution to F{A} T^B:' 1192 CALL CC_PRP(WORK(KOME1),WORK,1,1,0) 1193 END IF 1194 1195*---------------------------------------------------------------------* 1196* Compute contribution from <L_3|[[A,T^B_2],\tau_nu_2]|HF> 1197*---------------------------------------------------------------------* 1198 CALL DZERO(WORK(KOME2),NT1AMX*NT1AMX) 1199 1200 CALL CCSDT_E2AM(WORK(KOME2),WORK(KL3AM),WORK(KT2AM),WORK(KFOCKA)) 1201 1202 DO I = 1,NT1AMX 1203 DO J = 1,I 1204 IJ = NT1AMX*(I-1) + J 1205 NIJ = INDEX(I,J) 1206 OMEGA2(NIJ) = OMEGA2(NIJ) + WORK(KOME2+IJ-1) 1207 END DO 1208 END DO 1209 1210*---------------------------------------------------------------------* 1211* Compute [A,T^B_1] by AO-to-MO transformation of A with 1212* the response Lambda matrices: 1213*---------------------------------------------------------------------* 1214 CALL DCOPY(NORBT*NORBT,WORK(KFOCKA_AO),1,WORK(KFOCKAB1),1) 1215 CALL DCOPY(NORBT*NORBT,WORK(KFOCKA_AO),1,WORK(KFOCKAB),1) 1216 1217 CALL CC_FCKMO(WORK(KFOCKAB1),WORK(KLAMPB),WORK(KLAMH0), 1218 & WORK(KEND3),LWRK3,ISYMA,ISYMB,ISYM0) 1219 1220 CALL CC_FCKMO(WORK(KFOCKAB),WORK(KLAMP0),WORK(KLAMHB), 1221 & WORK(KEND3),LWRK3,ISYMA,ISYM0,ISYMB) 1222 1223 CALL DAXPY(NORBT*NORBT,ONE,WORK(KFOCKAB1),1,WORK(KFOCKAB),1) 1224 1225*---------------------------------------------------------------------* 1226* Compute triples result vector 1227* <L_3|[[A,T^B_1],\tau_nu_3]|HF> , 1228*---------------------------------------------------------------------* 1229 ! overwrite T3 vector 1230 KE3AM = KT3AM 1231 1232 CALL DZERO(WORK(KE3AM),NT1AMX*NT1AMX*NT1AMX) 1233 1234 L2INCL = .FALSE. 1235 CALL CCSDT_E3AM(WORK(KE3AM),WORK(KDUM),WORK(KL3AM), 1236 & WORK(KFOCKAB),L2INCL) 1237 1238*---------------------------------------------------------------------* 1239* Now we split: 1240* for IOPTRES < 5 we compute the effective result vector 1241* for IOPTRES = 5 we compute the contractions F{A} T^B T^C 1242*---------------------------------------------------------------------* 1243 IF (IOPTRES.GE.1 .AND. IOPTRES.LE.4) THEN 1244 1245 IOPT = 2 1246 Call CC_RDRSP('R0 ',0,ISYM0,IOPT,MODEL,WORK(KDUM),WORK(KEND3)) 1247 CALL CC_T2SQ(WORK(KEND3),WORK(KT2AM),ISYM0) 1248 1249 CALL DCOPY(NT1AMX,OMEGA1,1,OMEGA1EFF,1) 1250 CALL DCOPY(NT2AMX,OMEGA2,1,OMEGA2EFF,1) 1251 1252 FREQE = FREQL + FREQA + FREQB 1253 1254 CALL CC_LHPART_NODDY(OMEGA1EFF,OMEGA2EFF,WORK(KE3AM),-FREQE, 1255 & WORK(KFOCKD),WORK(KFIELD), 1256 & WORK(K0IOOOO),WORK(K0IOVVO), 1257 & WORK(K0IOOVV),WORK(K0IVVVV), 1258 & WORK(KT2AM),WORK(KINT1S0),WORK(KINT2S0), 1259 & WORK(KEND3),LWRK3) 1260 1261 ELSE IF (IOPTRES.EQ.5) THEN 1262 1263 SIGN = +1.0D0 1264 1265 CALL CCDOTRSP_NODDY(WORK(KOME1),WORK(KOME2),WORK(KE3AM),SIGN, 1266 & ITRAN,LISTDP,IDOTS,DOTPROD,MXVEC, 1267 & WORK(KLAMP0),WORK(KLAMH0), 1268 & WORK(KFOCK0),WORK(KFOCKD), 1269 & WORK(KXIAJB), WORK(KYIAJB), 1270 & WORK(KINT1T0),WORK(KINT2T0), 1271 & WORK(KINT1S0),WORK(KINT2S0), 1272 & 'CCSDT_FAMAT_NODDY',LOCDBG,LOCDBG, 1273 & .FALSE.,WORK(KEND3),LWRK3) 1274 1275 ELSE 1276 CALL QUIT('Illegal value for IOPTRES IN CCSDT_FAMAT_NODDY') 1277 END IF 1278 1279 CALL QEXIT('CCSDT_FAMAT_NODDY') 1280 1281 RETURN 1282 END 1283*---------------------------------------------------------------------* 1284* END OF SUBROUTINE CCSDT_FAMAT_NODDY * 1285*---------------------------------------------------------------------* 1286*=====================================================================* 1287 SUBROUTINE CCSDT_FA_SETUP(IFATRAN,IFADOTS,NFATRAN,MXVEC, 1288 & IDXL_FADEN,IDXB_FADEN,IDXC_FADEN, 1289 & NFADEN,MXFADEN, 1290 & IDXL_INTER,IDXR_INTER,LSTR_INTER, 1291 & NINTER,MXINTER, 1292 & LISTL,LISTO,LISTB,LISTC) 1293*---------------------------------------------------------------------* 1294* 1295* Purpose: setup loop structures to compute intermediates and 1296* densities needed to calculate contractions of the 1297* form F{A} t^B t^C 1298* 1299* Christof Haettig, May 2003 1300*---------------------------------------------------------------------* 1301 IMPLICIT NONE 1302#include "priunit.h" 1303#include "ccorb.h" 1304#include "ccsdsym.h" 1305#include "cclists.h" 1306#include "ccroper.h" 1307 1308 CHARACTER*(25) MSGDBG 1309 PARAMETER (MSGDBG = '[debug] CCSDT_FA_SETUP> ') 1310 LOGICAL LOCDBG 1311 PARAMETER (LOCDBG = .FALSE.) 1312 1313* input: 1314 CHARACTER*3 LISTL, LISTB, LISTC, LISTO 1315 INTEGER MXFADEN, NFATRAN, MXVEC, MXINTER 1316 INTEGER IFATRAN(MXDIM_FATRAN,NFATRAN), IFADOTS(MXVEC,NFATRAN) 1317 1318* output: 1319 CHARACTER*3 LSTR_INTER(*) 1320 INTEGER IDXL_FADEN(*), IDXB_FADEN(*), IDXC_FADEN(*), 1321 & IDXL_INTER(*), IDXR_INTER(*) 1322 INTEGER NFADEN, NINTER 1323 1324* local: 1325 CHARACTER*8 LABELA 1326 LOGICAL LPDBSA 1327 INTEGER ITRAN, IDLSTL, IOPERA, IDLSTB, IFILE, IRELAX, ISYRES, 1328 & IVEC, IDLSTC, ISYMTC, ISYCTR, ISYMTA, ISYMTB, 1329 & IFADEN, IDXLB, IDXLC, IDX 1330 1331* external functions: 1332 INTEGER ILSTSYM 1333 1334C ------ 1335C begin: 1336C ------ 1337C IF (LISTL.NE.'L0') THEN 1338C CALL QUIT('CCSDT_FA_DEN can not yet treat '//LISTL// 1339C & ' type vectors.') 1340C END IF 1341 1342C ------------------------------------------------ 1343C set up list of effective F{A} t^B t^C densities: 1344C ------------------------------------------------ 1345 NFADEN = 0 1346 DO ITRAN = 1, NFATRAN 1347 1348 IDLSTL = IFATRAN(1,ITRAN) 1349 IOPERA = IFATRAN(2,ITRAN) 1350 IDLSTB = IFATRAN(3,ITRAN) 1351 IFILE = IFATRAN(4,ITRAN) 1352 IRELAX = IFATRAN(5,ITRAN) 1353 1354 LABELA = LBLOPR(IOPERA) 1355 LPDBSA = LPDBSOP(IOPERA) 1356 1357 ISYCTR = ILSTSYM(LISTL,IDLSTL) 1358 ISYMTA = ILSTSYM(LISTO,IOPERA) 1359 ISYMTB = ILSTSYM(LISTB,IDLSTB) 1360 1361 ISYRES = MULD2H(MULD2H(ISYMTA,ISYMTB),ISYCTR) 1362 1363 IF ( IRELAX.GT.0 .OR. LPDBSA ) THEN 1364 CALL QUIT('Relaxed perturbations not yet implemented in '// 1365 & 'CCSDT_FA_DEN.') 1366 END IF 1367 1368 IVEC = 1 1369 DO WHILE(IFADOTS(IVEC,ITRAN).NE.0 .AND. IVEC.LE.MXVEC) 1370 IDLSTC = IFADOTS(IVEC,ITRAN) 1371 ISYMTC = ILSTSYM(LISTC,IDLSTC) 1372 1373 IF (ISYMTC.NE.ISYRES) THEN 1374 CALL QUIT('Symmetry mismatch in CCSDT_FA_DEN') 1375 END IF 1376 1377 ! check if density is already on to-do list 1378 IFADEN = 0 1379 DO IDX = 1, NFADEN 1380 IF (IDLSTL .EQ.IDXL_FADEN(IDX) .AND. 1381 & IDLSTB .EQ.IDXB_FADEN(IDX) .AND. 1382 & IDLSTC .EQ.IDXC_FADEN(IDX) ) THEN 1383 IFADEN = IDX 1384 END IF 1385 IF (IFADEN.EQ.0 .AND. LISTB.EQ.LISTC .AND. 1386 & IDLSTL .EQ.IDXL_FADEN(IDX) .AND. 1387 & IDLSTC .EQ.IDXB_FADEN(IDX) .AND. 1388 & IDLSTB .EQ.IDXC_FADEN(IDX) ) THEN 1389 IFADEN = IDX 1390 END IF 1391 END DO 1392 1393 ! if not found, put on to-do list for densities 1394 IF (IFADEN.EQ.0) THEN 1395 NFADEN = NFADEN + 1 1396 IF (NFADEN.GT.MXFADEN) CALL QUIT('NFADEN out of range') 1397 IDXL_FADEN(NFADEN) = IDLSTL 1398 IDXB_FADEN(NFADEN) = IDLSTB 1399 IDXC_FADEN(NFADEN) = IDLSTC 1400 END IF 1401 1402 IVEC = IVEC + 1 1403 END DO 1404 1405 END DO 1406 1407C ------------------------------------------------ 1408C set up list of intermediates that depend on 1409C pairs (L,B) and (L,C): 1410C ------------------------------------------------ 1411 NINTER = 0 1412 DO IFADEN = 1, NFADEN 1413 1414 ! check if pair (L,B) is already on the list 1415 IDXLB = 0 1416 DO IDX = 1, NINTER 1417 IF (IDXL_FADEN(IFADEN).EQ.IDXL_INTER(IDX) .AND. 1418 & IDXB_FADEN(IFADEN).EQ.IDXR_INTER(IDX) .AND. 1419 & LISTB .EQ.LSTR_INTER(IDX) ) THEN 1420 IDXLB = IDX 1421 END IF 1422 END DO 1423 1424 ! if not found, put on to-do list for intermediates 1425 IF (IDXLB.EQ.0) THEN 1426 NINTER = NINTER + 1 1427 IF (NINTER.GT.MXINTER) CALL QUIT('NINTER out of range') 1428 IDXL_INTER(NINTER) = IDXL_FADEN(IFADEN) 1429 IDXR_INTER(NINTER) = IDXB_FADEN(IFADEN) 1430 LSTR_INTER(NINTER) = LISTB 1431 END IF 1432 1433 ! check if pair (L,C) is already on the list 1434 IDXLC = 0 1435 DO IDX = 1, NINTER 1436 IF (IDXL_FADEN(IFADEN).EQ.IDXL_INTER(IDX) .AND. 1437 & IDXC_FADEN(IFADEN).EQ.IDXR_INTER(IDX) .AND. 1438 & LISTC .EQ.LSTR_INTER(IDX) ) THEN 1439 IDXLC = IDX 1440 END IF 1441 END DO 1442 1443 ! if not found, put on to-do list for intermediates 1444 IF (IDXLC.EQ.0) THEN 1445 NINTER = NINTER + 1 1446 IF (NINTER.GT.MXINTER) CALL QUIT('NINTER out of range') 1447 IDXL_INTER(NINTER) = IDXL_FADEN(IFADEN) 1448 IDXR_INTER(NINTER) = IDXC_FADEN(IFADEN) 1449 LSTR_INTER(NINTER) = LISTC 1450 END IF 1451 1452 END DO 1453 1454C ------------------------- 1455C if requested print lists: 1456C ------------------------- 1457 IF (LOCDBG) THEN 1458 WRITE(LUPRI,'(//,A)') 'list of F{A} T^B T^C densities:' 1459 WRITE(LUPRI,'(A)') '-------------------------------' 1460 WRITE(LUPRI,'(5X,A)') ' IDX L B C ' 1461 DO IFADEN = 1, NFADEN 1462 WRITE(LUPRI,'(5X,4I5)') IFADEN, IDXL_FADEN(IFADEN), 1463 & IDXB_FADEN(IFADEN), IDXC_FADEN(IFADEN) 1464 END DO 1465 WRITE(LUPRI,'(//,A)') 'list intermediates:' 1466 WRITE(LUPRI,'(A)') '-------------------' 1467 WRITE(LUPRI,'(5X,A)') ' IDX L R ' 1468 DO IFADEN = 1, NINTER 1469 WRITE(LUPRI,'(5X,2I5,3X,A3,I3)') IFADEN, IDXL_INTER(IFADEN), 1470 & LSTR_INTER(IFADEN), IDXR_INTER(IFADEN) 1471 END DO 1472 WRITE(LUPRI,'(//)') 1473 END IF 1474 1475 RETURN 1476 END 1477*---------------------------------------------------------------------* 1478* END OF SUBROUTINE CCSDT_FA_SETUP * 1479*---------------------------------------------------------------------* 1480*=====================================================================* 1481 SUBROUTINE CCSDT_FA_DEN(LISTL,LISTO,LISTB,LISTC,FADEN_NODDY, 1482 & NFATRAN,MXVEC,IFATRAN,IFADOTS,FACON, 1483 & FNDELD,FNCKJD,FNDKBC,FNTOC, 1484 & FN3VI,FNDKBC3,FN3FOPX,FN3FOP2X, 1485 & IDXL_INTER,IDXR_INTER, 1486 & LSTR_INTER,NINTETA, 1487 & IDXL_FADEN,IDXB_FADEN, 1488 & IDXC_FADEN,NFADEN, 1489 & IADR_INTER,IADR_DEN, 1490 & WORK,LWORK) 1491*---------------------------------------------------------------------* 1492* 1493* Purpose: compute triples contributions to to F{A} matrix 1494* contractions via effective densities 1495* 1496* Written by Christof Haettig, Mai 2003, Friedrichstal 1497* 1498*=====================================================================* 1499 IMPLICIT NONE 1500#include "ccsdsym.h" 1501#include "priunit.h" 1502#include "ccorb.h" 1503#include "dummy.h" 1504#include "ccr1rsp.h" 1505#include "ccroper.h" 1506#include "ccexci.h" 1507#include "cclists.h" 1508 1509 LOGICAL LOCDBG 1510 PARAMETER ( LOCDBG = .FALSE. ) 1511 1512 INTEGER ISYM0 1513 PARAMETER (ISYM0 = 1) 1514 1515 CHARACTER FNINTER*12, FNR2TP*12, FNDEFF*12 1516 PARAMETER (FNINTER='CCFADENINTER', FNR2TP = 'CCFADEN_R2TP', 1517 * FNDEFF ='CCFADEN_DEFF') 1518 1519* input/output: 1520 CHARACTER*3 LISTL, LISTO, LISTB, LISTC, LSTR_INTER(*) 1521 CHARACTER*(*) FNDELD, FNCKJD, FNDKBC, FNTOC, FN3VI, FNDKBC3 1522 CHARACTER*(*) FN3FOPX, FN3FOP2X 1523 LOGICAL FADEN_NODDY 1524 INTEGER NINTETA, NFADEN, MXVEC, NFATRAN, LWORK 1525 INTEGER IDXL_INTER(*), IDXR_INTER(*), IDXL_FADEN(*), 1526 * IDXB_FADEN(*), IDXC_FADEN(*), 1527 * IADR_INTER(3,*), IADR_DEN(*), 1528 * IFADOTS(MXVEC,*), IFATRAN(MXDIM_FATRAN,*) 1529 1530#if defined (SYS_CRAY) 1531 REAL FACON(MXVEC,NFATRAN) 1532 REAL WORK(*) 1533#else 1534 DOUBLE PRECISION FACON(MXVEC,NFATRAN) 1535 DOUBLE PRECISION WORK(*) 1536#endif 1537 1538* local: 1539 CHARACTER MODEL*(10), LISTR*3, LABELR*8, LABELA*8 1540 LOGICAL LPDBSA, CALL_ETA_FA_DEN, NEED_L1AM, NEED_L2TP, 1541 & NEED_T2TP, NEED_R2TP 1542 INTEGER ISYM, ILEN, ISYMFN, ISYMIM, NIMFN(8) 1543 INTEGER LUFOCK, LUINTER, LUR2TP, 1544 & LUDELD, LUDKBC, LUTOC, LU3VI, 1545 & LUDKBC3, LU3FOPX, LU3FOP2X, LUDEFF, LUCKJD 1546 INTEGER KEND1, KT1AM, KLAMP0, KLAMH0, KFOCK0, KT2TP, KL1AM, 1547 & KL2TP, IOPT, IDLSTL, ISYML, KEND2, LWRK2, IADRINT, 1548 & IDLSTR, ISYMR, INTETA, KR2TP, KEND3, LWRK3, ISYETA, 1549 & KDAB, KDIJ, KMMAT, ISYDEN, KDIA1, KDIA2, KDIA3, KDIA4, 1550 & KDIA5, KDIA6, IFADEN, IADRIA, ITRAN, IOPERA, IDLSTB, 1551 & IFILE, IRELAX, ISYMA, ISYMB, ISYRES, IVEC, IDLSTC, 1552 & ISYMC, KFOCK, KFOCKIA, IRREP, IERR, IMAT, LWRK1, IDX, 1553 & KEND0 1554 1555#if defined (SYS_CRAY) 1556 REAL FREQR, TRIPCON 1557 REAL TWO 1558#else 1559 DOUBLE PRECISION FREQR, TRIPCON 1560 DOUBLE PRECISION TWO 1561#endif 1562 PARAMETER ( TWO = 2.0D0 ) 1563 1564* external functions: 1565 INTEGER ILSTSYM 1566#if defined (SYS_CRAY) 1567 REAL DDOT 1568#else 1569 DOUBLE PRECISION DDOT 1570#endif 1571 1572 1573*---------------------------------------------------------------------* 1574* some initializations: 1575*---------------------------------------------------------------------* 1576 DO ISYM = 1, NSYM 1577 ILEN = 0 1578 DO ISYMFN = 1, NSYM 1579 ISYMIM = MULD2H(ISYM,ISYMFN) 1580 ILEN = ILEN + NMATIJ(ISYMIM)*NT1AM(ISYMFN) 1581 END DO 1582 NIMFN(ISYM) = ILEN 1583 END DO 1584 1585 IF (LOCDBG) THEN 1586 WRITE(LUPRI,*) 'Entered CCSDT_FA_DEN...' 1587 WRITE(LUPRI,*) 'LISTL,LISTB,LISTC:',LISTL,LISTB,LISTC 1588 END IF 1589*---------------------------------------------------------------------* 1590* allocate some zeroth-order intermediates: 1591*---------------------------------------------------------------------* 1592 KEND0 = 1 1593 1594 KT1AM = KEND0 1595 KLAMP0 = KT1AM + NT1AMX 1596 KLAMH0 = KLAMP0 + NLAMDT 1597 KEND0 = KLAMH0 + NLAMDT 1598 1599 KFOCK0 = KEND0 1600 KEND1 = KFOCK0 + N2BAST 1601 1602 LWRK1 = LWORK - KEND1 1603 1604 IF (LWRK1 .LT. 0) THEN 1605 CALL QUIT('Insufficient work space CCSDT_FA_DEN (1)') 1606 END IF 1607 1608*---------------------------------------------------------------------* 1609* read zeroth-order amplitudes and compute lambda matrices: 1610*---------------------------------------------------------------------* 1611 IOPT = 1 1612 CALL CC_RDRSP('R0 ',0,ISYM0,IOPT,MODEL,WORK(KT1AM),DUMMY) 1613 1614 CALL LAMMAT(WORK(KLAMP0),WORK(KLAMH0),WORK(KT1AM), 1615 & WORK(KEND1),LWRK1) 1616 1617*---------------------------------------------------------------------* 1618* get zeroth-order AO Fock in Lambda-MO basis: 1619*---------------------------------------------------------------------* 1620 LUFOCK = -1 1621 CALL GPOPEN(LUFOCK,'CC_FCKH','OLD',' ','UNFORMATTED',IDUMMY, 1622 & .FALSE.) 1623 REWIND(LUFOCK) 1624 READ(LUFOCK) (WORK(KFOCK0-1+I),I=1,N2BST(ISYM0)) 1625 CALL GPCLOSE(LUFOCK,'KEEP') 1626 1627 CALL CC_FCKMO(WORK(KFOCK0),WORK(KLAMP0),WORK(KLAMH0), 1628 & WORK(KEND1),LWRK1,ISYM0,ISYM0,ISYM0) 1629 1630*---------------------------------------------------------------------* 1631* open some files: 1632*---------------------------------------------------------------------* 1633 LUINTER = -1 1634 CALL WOPEN2(LUINTER,FNINTER,64,0) 1635 1636 LUDELD = -1 1637 LUCKJD = -1 1638 LUDKBC = -1 1639 LUTOC = -1 1640 LU3VI = -1 1641 LUDKBC3 = -1 1642 LU3FOPX = -1 1643 LU3FOP2X= -1 1644 1645 CALL WOPEN2(LUDELD,FNDELD,64,0) 1646 CALL WOPEN2(LUCKJD,FNCKJD,64,0) 1647 CALL WOPEN2(LUDKBC,FNDKBC,64,0) 1648 CALL WOPEN2(LUTOC,FNTOC,64,0) 1649 CALL WOPEN2(LU3VI,FN3VI,64,0) 1650 CALL WOPEN2(LUDKBC3,FNDKBC3,64,0) 1651 CALL WOPEN2(LU3FOPX,FN3FOPX,64,0) 1652 CALL WOPEN2(LU3FOP2X,FN3FOP2X,64,0) 1653 1654*---------------------------------------------------------------------* 1655* loop over the intermediates that need to be computed: 1656*---------------------------------------------------------------------* 1657 IADRINT = 1 1658 1659 DO INTETA = 1, NINTETA 1660 1661 IF (LISTL.EQ.'L0 ') THEN 1662 IDLSTL = 0 1663 ISYML = 1 1664 ELSE 1665 IDLSTL = IDXL_INTER(INTETA) 1666 ISYML = ILSTSYM(LISTL,IDLSTL) 1667 END IF 1668 1669 LISTR = LSTR_INTER(INTETA) 1670 IDLSTR = IDXR_INTER(INTETA) 1671 ISYMR = ILSTSYM(LISTR,IDLSTR) 1672 1673 ! check for modules which needed extra arrays / preparations 1674 CALL_ETA_FA_DEN = ( (LISTL.EQ.'L0 ' .OR. LISTL.EQ.'LE ') .AND. 1675 * (LISTR.EQ.'R1 ' .OR. LISTR.EQ.'RE ' 1676 * .OR. LISTR.EQ.'R2 ') ) 1677 1678 NEED_L1AM = CALL_ETA_FA_DEN 1679 NEED_L2TP = CALL_ETA_FA_DEN 1680 NEED_T2TP = CALL_ETA_FA_DEN 1681 NEED_R2TP = CALL_ETA_FA_DEN 1682 1683 1684c ------------------------------- 1685c allocate space for multipliers: 1686c ------------------------------- 1687 KEND2 = KEND1 1688 1689 IF (NEED_T2TP) THEN 1690 KT2TP = KEND2 1691 KEND2 = KT2TP + NT2SQ(ISYM0) 1692 END IF 1693 1694 IF (NEED_L1AM) THEN 1695 KL1AM = KEND2 1696 KEND2 = KL1AM + NT1AM(ISYML) 1697 END IF 1698 1699 IF (NEED_L2TP) THEN 1700 KL2TP = KEND2 1701 KEND2 = KL2TP + NT2SQ(ISYML) 1702 END IF 1703 1704 LWRK2 = LWORK - KEND2 1705 IF (LWRK2 .LT. NT2SQ(ISYML)) THEN 1706 CALL QUIT('Insufficient work space CCSDT_FA_DEN (3)') 1707 END IF 1708 1709c ------------------------------------------------ 1710c read and resort zeroth-order doubles amplitudes: 1711c ------------------------------------------------ 1712 IF (NEED_T2TP) THEN 1713 IF (LWRK2 .LT. NT2SQ(ISYM0)) 1714 & CALL QUIT('Insufficient work space CCSDT_FA_DEN (2a)') 1715 1716 IOPT = 2 1717 CALL CC_RDRSP('R0',0,ISYM0,IOPT,MODEL,DUMMY,WORK(KT2TP)) 1718 CALL CC_T2SQ(WORK(KT2TP),WORK(KEND2),ISYM0) 1719 CALL CC3_T2TP(WORK(KT2TP),WORK(KEND2),ISYM0) 1720 IF (LOCDBG) WRITE(LUPRI,*) 'Norm of T2TP ', 1721 * DDOT(NT2SQ(ISYM0),WORK(KT2TP),1,WORK(KT2TP),1) 1722 END IF 1723 1724c ------------------------- 1725c read singles multipliers: 1726c ------------------------- 1727 IF (NEED_L1AM) THEN 1728 IOPT = 1 1729 CALL CC_RDRSP(LISTL,IDLSTL,ISYML,IOPT,MODEL,WORK(KL1AM),DUMMY) 1730 END IF 1731 1732c ------------------------------------ 1733c read and resort doubles multipliers: 1734c ------------------------------------ 1735 IF (NEED_L2TP) THEN 1736 IF (LWRK2 .LT. NT2SQ(ISYM0)) 1737 & CALL QUIT('Insufficient work space CCSDT_FA_DEN (2b)') 1738 1739 IOPT = 3 1740 CALL CC_RDRSP(LISTL,IDLSTL,ISYML,IOPT,MODEL, 1741 * WORK(KL1AM),WORK(KL2TP)) 1742 CALL CC_T2SQ(WORK(KL2TP),WORK(KEND2),ISYML) 1743 CALL CC3_T2TP(WORK(KL2TP),WORK(KEND2),ISYML) 1744 IF (LOCDBG) WRITE(LUPRI,*) 'Norm of L2TP ', 1745 * DDOT(NT2SQ(ISYML),WORK(KL2TP),1,WORK(KL2TP),1) 1746 END IF 1747 1748c --------------------------------------------- 1749c read and resort doubles response multipliers: 1750c --------------------------------------------- 1751 IF (NEED_R2TP) THEN 1752 LUR2TP = -1 1753 CALL WOPEN2(LUR2TP,FNR2TP,64,0) 1754 1755 KR2TP = KEND2 1756 KEND3 = KR2TP + NT2SQ(ISYMR) 1757 LWRK3 = LWORK - KEND3 1758 IF (LWRK3.LT.NT2AM(ISYMR)) 1759 & CALL QUIT('Out of memory in CCSDT_FA_DEN. (4a)') 1760 1761 IOPT = 2 1762 CALL CC_RDRSP(LISTR,IDLSTR,ISYMR,IOPT,MODEL, 1763 & DUMMY,WORK(KEND3)) 1764 CALL CCLR_DIASCL(WORK(KEND3),TWO,ISYMR) 1765 CALL CC_T2SQ(WORK(KEND3),WORK(KR2TP),ISYMR) 1766 1767 CALL PUTWA2(LUR2TP,FNR2TP,WORK(KR2TP),1,NT2SQ(ISYMR)) 1768 END IF 1769 1770c ------------------------------------------------------ 1771c allocate memory for intermediates and initialize them: 1772c ------------------------------------------------------ 1773 ISYETA = MULD2H(ISYML,ISYMR) 1774 1775 KDAB = KEND2 1776 KDIJ = KDAB + NMATAB(ISYETA) 1777 KMMAT = KDIJ + NMATIJ(ISYETA) 1778 KEND3 = KMMAT + NIMFN(ISYETA) 1779 1780 LWRK3 = LWORK - KEND3 1781 IF (LWRK3.LT.0) CALL QUIT('Out of memory in CCSDT_FA_DEN. (4)') 1782 1783 CALL DZERO(WORK(KDAB), NMATAB(ISYETA)) 1784 CALL DZERO(WORK(KDIJ), NMATIJ(ISYETA)) 1785 CALL DZERO(WORK(KMMAT),NIMFN(ISYETA)) 1786 1787c ------------------------------------------------------ 1788c compute triples intermediates: D(ab), D(ij), M(imfn) 1789c ------------------------------------------------------ 1790 IF ( ( (LISTL(1:3).EQ.'L0 ') .OR. (LISTL(1:3).EQ.'LE ') ) .AND. 1791 * ( (LISTR(1:3).EQ.'R1 ') .OR. (LISTR(1:3).EQ.'RE ') ) ) THEN 1792 1793 IF (LISTR(1:3).EQ.'R1 ') THEN 1794 FREQR = FRQLRT(IDLSTR) 1795 LABELR = LRTLBL(IDLSTR) 1796 IF (LORXLRT(IDLSTR)) CALL QUIT('NO ORBITAL RELAXED '// 1797 & 'PERTURBATION IMPLEMENTED IN CCSDT_FA_DEN.') 1798 ELSE IF (LISTR(1:3).EQ.'RE ') THEN 1799 FREQR = EIGVAL(IDLSTR) 1800 LABELR = '- none -' 1801 ELSE 1802 CALL QUIT('Illegal right vector in CCSDT_FA_DEN.') 1803 END IF 1804 1805 1806 IF (FADEN_NODDY) THEN 1807 WRITE(LUPRI,*) 'No noddy code for the case in CCSDT_FA_DEN:' 1808 WRITE(LUPRI,*) 'LISTL,LISTR:',LISTL,LISTR 1809 WRITE(LUPRI,*) 'the real code will be used instead...' 1810 END IF 1811 1812 !frequency of LISTL is handled inside 1813 CALL CCSDT_ETA_FA_DEN(LISTL,IDLSTL,ISYML, 1814 * LISTR,IDLSTR,ISYMR,FREQR,LABELR, 1815 * .FALSE.,0, 1816 * WORK(KLAMP0),WORK(KLAMH0),WORK(KFOCK0), 1817 * WORK(KT2TP),WORK(KL2TP),WORK(KL1AM), 1818 * ISYETA,WORK(KDAB),WORK(KDIJ), 1819 * .FALSE.,DUMMY, 1820 * .FALSE.,DUMMY,.TRUE.,WORK(KMMAT), 1821 * WORK(KEND3),LWRK3, 1822 * LUDELD, FNDELD, LUCKJD, FNCKJD, 1823 * LUDKBC, FNDKBC, LUTOC, FNTOC, 1824 * LU3VI, FN3VI, LUDKBC3,FNDKBC3, 1825 * LU3FOPX, FN3FOPX, LU3FOP2X,FN3FOP2X, 1826 * LUR2TP, FNR2TP ) 1827C 1828 1829 ELSE IF ( LISTL.EQ.'L0 '.AND. 1830 * (LISTR.EQ.'R2 '.OR. LISTR.EQ.'ER1') ) THEN 1831 1832C 1833 1834 IF (FADEN_NODDY) THEN 1835 1836 CALL CCSDT_ADEN_NODDY(LISTL,IDLSTL,LISTR,IDLSTR, 1837 * WORK(KLAMP0),WORK(KLAMH0), 1838 * WORK(KFOCK0), 1839 * WORK(KDIJ),WORK(KDAB), 1840 * .FALSE.,DUMMY,.TRUE.,WORK(KMMAT), 1841 * WORK(KEND3),LWRK3, 1842 * LUDELD,FNDELD,LUCKJD,FNCKJD,LUDKBC, 1843 * FNDKBC,LUTOC,FNTOC,LU3VI,FN3VI, 1844 * LUDKBC3,FNDKBC3,LU3FOPX,FN3FOPX, 1845 * LU3FOP2X,FN3FOP2X) 1846C 1847 ELSE 1848 1849 !FIX LISTL WITH IF STATEMENTS INSIDE !!!! 1850 CALL CC3_ADEN_CUB_T0('L1 ',1,LISTR,IDLSTR, 1851 * WORK(KLAMP0),WORK(KLAMH0),WORK(KFOCK0), 1852 * WORK(KDIJ),WORK(KDAB),ISYETA, 1853 * .TRUE.,WORK(KMMAT), 1854 * WORK(KEND3),LWRK3, 1855 * LUDELD,FNDELD,LUCKJD,FNCKJD,LUDKBC, 1856 * FNDKBC,LUTOC,FNTOC,LU3VI,FN3VI, 1857 * LUDKBC3,FNDKBC3,LU3FOPX,FN3FOPX, 1858 * LU3FOP2X,FN3FOP2X) 1859C 1860 END IF 1861 1862 ELSE IF ( (LISTL.EQ.'L1 '.OR. LISTL.EQ.'LE ') .AND. 1863 * (LISTR.EQ.'R1 '.OR. LISTR.EQ.'RE ') ) THEN 1864 1865 IF (FADEN_NODDY) THEN 1866 CALL CCSDT_ADEN_NODDY(LISTL,IDLSTL,LISTR,IDLSTR, 1867 * WORK(KLAMP0),WORK(KLAMH0), 1868 * WORK(KFOCK0), 1869 * WORK(KDIJ),WORK(KDAB), 1870 * .FALSE.,DUMMY,.TRUE.,WORK(KMMAT), 1871 * WORK(KEND3),LWRK3, 1872 * LUDELD,FNDELD,LUCKJD,FNCKJD,LUDKBC, 1873 * FNDKBC,LUTOC,FNTOC,LU3VI,FN3VI, 1874 * LUDKBC3,FNDKBC3,LU3FOPX,FN3FOPX, 1875 * LU3FOP2X,FN3FOP2X) 1876C 1877 1878 ELSE 1879 CALL CC3_ADEN(LISTL,IDLSTL,LISTR,IDLSTR, 1880 * WORK(KLAMP0),WORK(KLAMH0),WORK(KFOCK0), 1881 * WORK(KDIJ),WORK(KDAB), 1882 * .FALSE.,DUMMY,ISYETA, 1883 * .TRUE.,WORK(KMMAT), 1884 * WORK(KEND3),LWRK3, 1885 * LUDELD,FNDELD,LUCKJD,FNCKJD,LUDKBC, 1886 * FNDKBC,LUTOC,FNTOC,LU3VI,FN3VI, 1887 * LUDKBC3,FNDKBC3,LU3FOPX,FN3FOPX, 1888 * LU3FOP2X,FN3FOP2X) 1889C 1890 END IF 1891 1892 ELSE 1893 CALL QUIT('Encountered non-implemented case in CCSDT_FA_DEN.') 1894 END IF 1895 1896c ------------------------------------------------- 1897c put intermediates on file and remember addresses: 1898c ------------------------------------------------- 1899 IADR_INTER(1,INTETA) = IADRINT 1900 CALL PUTWA2(LUINTER,FNINTER,WORK(KDAB),IADRINT,NMATAB(ISYETA)) 1901 IADRINT = IADRINT + NMATAB(ISYETA) 1902 1903 IADR_INTER(2,INTETA) = IADRINT 1904 CALL PUTWA2(LUINTER,FNINTER,WORK(KDIJ),IADRINT,NMATIJ(ISYETA)) 1905 IADRINT = IADRINT + NMATIJ(ISYETA) 1906 1907 IADR_INTER(3,INTETA) = IADRINT 1908 CALL PUTWA2(LUINTER,FNINTER,WORK(KMMAT),IADRINT,NIMFN(ISYETA)) 1909 IADRINT = IADRINT + NIMFN(ISYETA) 1910 1911 1912 IF (NEED_R2TP) CALL WCLOSE2(LUR2TP,FNR2TP,'DELETE') 1913 1914 END DO 1915 1916*---------------------------------------------------------------------* 1917* close/delete some file: 1918*---------------------------------------------------------------------* 1919 CALL WCLOSE2(LUDELD,FNDELD,'KEEP') 1920 CALL WCLOSE2(LUCKJD,FNCKJD,'KEEP') 1921 CALL WCLOSE2(LUTOC,FNTOC,'KEEP') 1922 CALL WCLOSE2(LUDKBC,FNDKBC,'KEEP') 1923 CALL WCLOSE2(LU3VI,FN3VI,'KEEP') 1924 CALL WCLOSE2(LUDKBC3,FNDKBC3,'KEEP') 1925 CALL WCLOSE2(LU3FOPX,FN3FOPX,'KEEP') 1926 CALL WCLOSE2(LU3FOP2X,FN3FOP2X,'KEEP') 1927 1928*---------------------------------------------------------------------* 1929* now compute the F{A} densities from the intermediates: 1930*---------------------------------------------------------------------* 1931 LUDEFF = -1 1932 CALL WOPEN2(LUDEFF,FNDEFF,64,0) 1933 1934 IADRIA = 1 1935 1936 DO IFADEN = 1, NFADEN 1937 1938 IDLSTL = IDXL_FADEN(IFADEN) 1939 IDLSTB = IDXB_FADEN(IFADEN) 1940 IDLSTC = IDXC_FADEN(IFADEN) 1941 1942 ISYML = ILSTSYM(LISTL,IDLSTL) 1943 ISYMB = ILSTSYM(LISTB,IDLSTB) 1944 ISYMC = ILSTSYM(LISTC,IDLSTC) 1945 1946 ISYDEN = MULD2H(MULD2H(ISYML,ISYMB),ISYMC) 1947 1948 KDIA1 = KEND0 1949 KEND1 = KDIA1 + NT1AM(ISYDEN) 1950 1951 IF (LOCDBG) THEN 1952 KDIA2 = KEND1 1953 KDIA3 = KDIA2 + NT1AM(ISYDEN) 1954 KDIA4 = KDIA3 + NT1AM(ISYDEN) 1955 KDIA5 = KDIA4 + NT1AM(ISYDEN) 1956 KDIA6 = KDIA5 + NT1AM(ISYDEN) 1957 KEND1 = KDIA6 + NT1AM(ISYDEN) 1958 ELSE 1959 KDIA2 = KDIA1 1960 KDIA3 = KDIA1 1961 KDIA4 = KDIA1 1962 KDIA5 = KDIA1 1963 KDIA6 = KDIA1 1964 END IF 1965 1966 LWRK1 = LWORK - KEND1 1967 IF (LWRK1 .LT. 0) THEN 1968 CALL QUIT('Insufficient work space CCSDT_FA_DEN (5)') 1969 END IF 1970 1971 CALL DZERO(WORK(KDIA1),NT1AM(ISYDEN)) 1972 IF (LOCDBG) THEN 1973 CALL DZERO(WORK(KDIA2),NT1AM(ISYDEN)) 1974 CALL DZERO(WORK(KDIA3),NT1AM(ISYDEN)) 1975 CALL DZERO(WORK(KDIA4),NT1AM(ISYDEN)) 1976 CALL DZERO(WORK(KDIA5),NT1AM(ISYDEN)) 1977 CALL DZERO(WORK(KDIA6),NT1AM(ISYDEN)) 1978 END IF 1979 1980 CALL CCSDT_FA_DEN1(WORK(KDIA1),WORK(KDIA2),WORK(KDIA3), 1981 & IDLSTL,ISYML, 1982 & LISTB,IDLSTB,ISYMB, 1983 & LISTC,IDLSTC,ISYMC, 1984 & IDXL_INTER,IDXR_INTER,LSTR_INTER, 1985 & IADR_INTER,LUINTER,FNINTER, 1986 & NIMFN,NINTETA,WORK(KEND1),LWRK1) 1987 1988 CALL CCSDT_FA_DEN1(WORK(KDIA4),WORK(KDIA5),WORK(KDIA6), 1989 & IDLSTL,ISYML, 1990 & LISTC,IDLSTC,ISYMC, 1991 & LISTB,IDLSTB,ISYMB, 1992 & IDXL_INTER,IDXR_INTER,LSTR_INTER, 1993 & IADR_INTER,LUINTER,FNINTER, 1994 & NIMFN,NINTETA,WORK(KEND1),LWRK1) 1995 1996 IF (LOCDBG) THEN 1997 WRITE(LUPRI,*) 'triples contribution 1 to F{A} density:', 1998 & DDOT(NT1AM(ISYDEN),WORK(KDIA1),1,WORK(KDIA1),1) 1999 CALL CC_PRP(WORK(KDIA1),DUMMY,1,1,0) 2000 WRITE(LUPRI,*) 'triples contribution 2 to F{A} density:', 2001 & DDOT(NT1AM(ISYDEN),WORK(KDIA2),1,WORK(KDIA2),1) 2002 CALL CC_PRP(WORK(KDIA2),DUMMY,1,1,0) 2003 WRITE(LUPRI,*) 'triples contribution 3 to F{A} density:', 2004 & DDOT(NT1AM(ISYDEN),WORK(KDIA3),1,WORK(KDIA3),1) 2005 CALL CC_PRP(WORK(KDIA3),DUMMY,1,1,0) 2006 WRITE(LUPRI,*) 'triples contribution 4 to F{A} density:', 2007 & DDOT(NT1AM(ISYDEN),WORK(KDIA4),1,WORK(KDIA4),1) 2008 CALL CC_PRP(WORK(KDIA4),DUMMY,1,1,0) 2009 WRITE(LUPRI,*) 'triples contribution 5 to F{A} density:', 2010 & DDOT(NT1AM(ISYDEN),WORK(KDIA5),1,WORK(KDIA5),1) 2011 CALL CC_PRP(WORK(KDIA5),DUMMY,1,1,0) 2012 WRITE(LUPRI,*) 'triples contribution 6 to F{A} density:', 2013 & DDOT(NT1AM(ISYDEN),WORK(KDIA6),1,WORK(KDIA6),1) 2014 CALL CC_PRP(WORK(KDIA6),DUMMY,1,1,0) 2015 2016 ! sum up all 6 contributions 2017 CALL DAXPY(NT1AM(ISYDEN),1.0D0,WORK(KDIA2),1,WORK(KDIA1),1) 2018 CALL DAXPY(NT1AM(ISYDEN),1.0D0,WORK(KDIA3),1,WORK(KDIA1),1) 2019 CALL DAXPY(NT1AM(ISYDEN),1.0D0,WORK(KDIA4),1,WORK(KDIA1),1) 2020 CALL DAXPY(NT1AM(ISYDEN),1.0D0,WORK(KDIA5),1,WORK(KDIA1),1) 2021 CALL DAXPY(NT1AM(ISYDEN),1.0D0,WORK(KDIA6),1,WORK(KDIA1),1) 2022 2023 WRITE(LUPRI,*) 'complete F{A} density:', 2024 & DDOT(NT1AM(ISYDEN),WORK(KDIA1),1,WORK(KDIA1),1) 2025 CALL CC_PRP(WORK(KDIA1),DUMMY,1,1,0) 2026 END IF 2027 2028 ! store density on file 2029 IADR_DEN(IFADEN) = IADRIA 2030 CALL PUTWA2(LUDEFF,FNDEFF,WORK(KDIA1),IADRIA,NT1AM(ISYDEN)) 2031 IADRIA = IADRIA + NT1AM(ISYDEN) 2032 2033 END DO 2034 2035 CALL WCLOSE2(LUINTER,FNINTER,'DELETE') 2036 2037*---------------------------------------------------------------------* 2038* now compute from the densities the contractions F{A} t^B t^C: 2039*---------------------------------------------------------------------* 2040 DO ITRAN = 1, NFATRAN 2041 2042 IDLSTL = IFATRAN(1,ITRAN) 2043 IOPERA = IFATRAN(2,ITRAN) 2044 IDLSTB = IFATRAN(3,ITRAN) 2045 IFILE = IFATRAN(4,ITRAN) 2046 IRELAX = IFATRAN(5,ITRAN) 2047 2048 LABELA = LBLOPR(IOPERA) 2049 LPDBSA = LPDBSOP(IOPERA) 2050 2051 ISYML = ILSTSYM(LISTL,IDLSTL) 2052 ISYMA = ILSTSYM(LISTO,IOPERA) 2053 ISYMB = ILSTSYM(LISTB,IDLSTB) 2054 2055 ISYRES = MULD2H(MULD2H(ISYMA,ISYMB),ISYML) 2056 2057 IF ( IRELAX.GT.0 .OR. LPDBSA ) THEN 2058 CALL QUIT('Relaxed perturbations not yet implemented in '// 2059 & 'CCSDT_FA_DEN.') 2060 END IF 2061 2062 IVEC = 1 2063 DO WHILE(IFADOTS(IVEC,ITRAN).NE.0 .AND. IVEC.LE.MXVEC) 2064 IDLSTC = IFADOTS(IVEC,ITRAN) 2065 ISYMC = ILSTSYM(LISTC,IDLSTC) 2066 2067 IF (ISYMC.NE.ISYRES) THEN 2068 CALL QUIT('Symmetry mismatch in CCSDT_FA_DEN') 2069 END IF 2070 2071 ! find index of the density needed 2072 IFADEN = 0 2073 DO IDX = 1, NFADEN 2074 IF (IDLSTL .EQ.IDXL_FADEN(IDX) .AND. 2075 & IDLSTB .EQ.IDXB_FADEN(IDX) .AND. 2076 & IDLSTC .EQ.IDXC_FADEN(IDX) ) THEN 2077 IFADEN = IDX 2078 END IF 2079 IF (IFADEN.EQ.0 .AND. LISTB.EQ.LISTC .AND. 2080 & IDLSTL .EQ.IDXL_FADEN(IDX) .AND. 2081 & IDLSTC .EQ.IDXB_FADEN(IDX) .AND. 2082 & IDLSTB .EQ.IDXC_FADEN(IDX) ) THEN 2083 IFADEN = IDX 2084 END IF 2085 END DO 2086 IF (IFADEN.LE.0) CALL QUIT('Fatal error in CCSDT_FA_DEN.') 2087 2088 ISYDEN = MULD2H(MULD2H(ISYML,ISYMB),ISYMC) 2089 2090 IF (ISYDEN.EQ.ISYMA) THEN 2091 KDIA1 = KEND0 2092 KFOCK = KDIA1 + NT1AM(ISYDEN) 2093 KFOCKIA = KFOCK + N2BST(ISYDEN) 2094 KEND1 = KFOCKIA + NT1AM(ISYDEN) 2095 LWRK1 = LWORK - KEND1 2096 2097 IADRIA = IADR_DEN(IFADEN) 2098 CALL GETWA2(LUDEFF,FNDEFF,WORK(KDIA1),IADRIA,NT1AM(ISYDEN)) 2099 2100 CALL CCPRPAO(LABELA,.TRUE.,WORK(KFOCK),IRREP,IMAT,IERR, 2101 * WORK(KEND1),LWRK1) 2102 IF ((IERR.GT.0) .OR. (IERR.EQ.0 .AND. IRREP.NE.ISYMA)) THEN 2103 WRITE(LUPRI,*) 'ISYMA :',ISYMA 2104 WRITE(LUPRI,*) 'IRREP :',IRREP 2105 WRITE(LUPRI,*) 'IERR :',IERR 2106 WRITE(LUPRI,*) 'LABEL :',LABELA 2107 CALL QUIT('CCSDT_FA_DEN: error reading operator ') 2108 ELSE IF (IERR.LT.0) THEN 2109 CALL DZERO(WORK(KFOCK),N2BST(ISYMA)) 2110 END IF 2111 2112 ! transform property integrals to Lambda-MO basis and resort 2113 CALL CC_FCKMO(WORK(KFOCK),WORK(KLAMP0),WORK(KLAMH0), 2114 & WORK(KEND1),LWRK1,ISYMA,1,1) 2115 CALL CC_FOCK_RESORT(DUMMY,.FALSE.,WORK(KFOCKIA),.TRUE., 2116 & DUMMY,.FALSE.,DUMMY,.FALSE., 2117 & WORK(KFOCK),ISYMA) 2118 2119 TRIPCON = DDOT(NT1AM(ISYMA),WORK(KDIA1),1,WORK(KFOCKIA),1) 2120 2121 FACON(IVEC,ITRAN) = FACON(IVEC,ITRAN) + TRIPCON 2122 END IF 2123 2124 IVEC = IVEC + 1 2125 END DO 2126 2127 END DO 2128 2129 RETURN 2130 END 2131*---------------------------------------------------------------------* 2132* END OF SUBROUTINE CCSDT_FA_DEN * 2133*---------------------------------------------------------------------* 2134*=====================================================================* 2135 SUBROUTINE CCSDT_FA_DEN1(DIA1,DIA2,DIA3,IDLSTL,ISYML, 2136 & LISTB,IDLSTB,ISYMB, 2137 & LISTC,IDLSTC,ISYMC, 2138 & IDXL_INTER,IDXR_INTER,LSTR_INTER, 2139 & IADR_INTER,LUINTER,FNINTER, 2140 & NIMFN,NINTER,WORK,LWORK) 2141*---------------------------------------------------------------------* 2142* 2143* Purpose: compute contribution to F{A} density from 2144* precomputed intermediates on file 2145* 2146* Written by Christof Haettig, Mai 2003, Friedrichstal 2147* 2148*=====================================================================* 2149 IMPLICIT NONE 2150#include "priunit.h" 2151#include "ccorb.h" 2152#include "ccsdsym.h" 2153#include "dummy.h" 2154 2155 CHARACTER FNINTER*(*) 2156 CHARACTER*3 LSTR_INTER(*), LISTB, LISTC 2157 2158* input/output: 2159 INTEGER LWORK, NINTER, NIMFN(8) 2160 INTEGER IDLSTL, ISYML, IDLSTB, ISYMB, IDLSTC, ISYMC, LUINTER 2161 INTEGER IDXL_INTER(*), IDXR_INTER(*), IADR_INTER(3,*) 2162 2163#if defined (SYS_CRAY) 2164 REAL DIA1(*), DIA2(*), DIA3(*) 2165 REAL WORK(*) 2166 REAL ONE, TWO 2167#else 2168 DOUBLE PRECISION DIA1(*), DIA2(*), DIA3(*) 2169 DOUBLE PRECISION WORK(*) 2170 DOUBLE PRECISION ONE, TWO 2171#endif 2172 PARAMETER ( ONE=1.0D0, TWO=2.0D0 ) 2173 2174* local: 2175 CHARACTER CDUMMY*1, MODEL*10 2176 INTEGER ISYETA, ISYDEN, KDAE, KDIJ, KMMAT, KC1AM, KC2TP, KEND1, 2177 & LWRK1, IDXLB, IDX, IADRINT, IOPT, IOPTT2, ISYMA, ISYMI, 2178 & ISYME, ISYMJ, KOFFDEA, KOFFDIJ, KOFFDIA, KOFFTEI, 2179 & KOFFTAJ, NVIRA, NVIRE, NRHFI 2180 2181 2182 2183c ------------------ 2184c memory allocation: 2185c ------------------ 2186 ISYETA = MULD2H(ISYML, ISYMB) 2187 ISYDEN = MULD2H(ISYETA,ISYMC) 2188 2189 KDAE = 1 2190 KDIJ = KDAE + NMATAB(ISYETA) 2191 KMMAT = KDIJ + NMATIJ(ISYETA) 2192 KC1AM = KMMAT + NIMFN(ISYETA) 2193 KC2TP = KC1AM + NT1AM(ISYMC) 2194 KEND1 = KC2TP + NT2SQ(ISYMC) 2195 2196 LWRK1 = LWORK - KEND1 2197 IF (LWRK1.LT.NT2SQ(ISYMC)) THEN 2198 CALL QUIT('Out of memory in CCSDT_FA_DEN1.') 2199 END IF 2200 2201c ----------------------------------------------------- 2202c restore intermediates depending on L and B from file 2203c ----------------------------------------------------- 2204 ! get index for the intermediates depending on L and B: 2205 IDXLB = 0 2206 DO IDX = 1, NINTER 2207 IF (IDLSTL.EQ.IDXL_INTER(IDX) .AND. 2208 & IDLSTB.EQ.IDXR_INTER(IDX) .AND. 2209 & LISTB .EQ.LSTR_INTER(IDX) ) THEN 2210 IDXLB = IDX 2211 END IF 2212 END DO 2213 IF (IDXLB.LE.0) CALL QUIT('Fatal error in CCSDT_FA_DEN1.') 2214 2215 ! read intermediates 2216 IADRINT = IADR_INTER(1,IDXLB) 2217 CALL GETWA2(LUINTER,FNINTER,WORK(KDAE),IADRINT,NMATAB(ISYETA)) 2218 2219 IADRINT = IADR_INTER(2,IDXLB) 2220 CALL GETWA2(LUINTER,FNINTER,WORK(KDIJ),IADRINT,NMATIJ(ISYETA)) 2221 2222 IADRINT = IADR_INTER(3,IDXLB) 2223 CALL GETWA2(LUINTER,FNINTER,WORK(KMMAT),IADRINT,NIMFN(ISYETA)) 2224 2225c -------------------------------------------- 2226c read C vector and resort doubles amplitudes: 2227c -------------------------------------------- 2228 IOPT = 3 2229 CALL CC_RDRSP(LISTC,IDLSTC,ISYMC,IOPT,MODEL, 2230 & WORK(KC1AM),WORK(KC2TP)) 2231 CALL CCLR_DIASCL(WORK(KC2TP),TWO,ISYMC) 2232 CALL CC_T2SQ(WORK(KC2TP),WORK(KEND1),ISYMC) 2233 CALL CC3_T2TP(WORK(KC2TP),WORK(KEND1),ISYMC) 2234 2235 2236C --------------------------------------------- 2237C D(ia) <-- D(ia) + sum_fnm t(am,fn) y^M(imfn): 2238C --------------------------------------------- 2239C ! turn sign for this contribution to D(ia) 2240C CALL DSCAL(NT2SQ(ISYMC),-1.0D0,WORK(KC2TP),1) 2241 2242 IOPTT2 = 0 2243 CALL CCSDT_ETA_TM2(DIA3,ISYDEN,WORK(KMMAT),ISYETA, 2244 & WORK(KC2TP),IDUMMY,CDUMMY,ISYMC,IOPTT2, 2245 & WORK(KEND1),LWRK1) 2246 2247 2248C ---------------------------------------------------------------- 2249C D(ia) <-- D(ia) - sum_e y^t(ei) D^0(ea) + sum_j y^t(aj) D^0(ij): 2250C ---------------------------------------------------------------- 2251 DO ISYMA = 1, NSYM 2252 ISYMI = MULD2H(ISYDEN,ISYMA) 2253 ISYME = MULD2H(ISYETA,ISYMA) 2254 ISYMJ = MULD2H(ISYETA,ISYMI) 2255 2256 KOFFDEA = KDAE + IMATAB(ISYME,ISYMA) 2257 KOFFDIJ = KDIJ + IMATIJ(ISYMI,ISYMJ) 2258 KOFFDIA = 1 + IT1AM(ISYMA,ISYMI) 2259 KOFFTEI = KC1AM + IT1AM(ISYME,ISYMI) 2260 KOFFTAJ = KC1AM + IT1AM(ISYMA,ISYMJ) 2261 2262 NVIRA = MAX(NVIR(ISYMA),1) 2263 NVIRE = MAX(NVIR(ISYME),1) 2264 NRHFI = MAX(NRHF(ISYMI),1) 2265 2266 CALL DGEMM('T','N',NVIR(ISYMA),NRHF(ISYMI),NVIR(ISYME), 2267 & -ONE,WORK(KOFFDEA),NVIRE,WORK(KOFFTEI),NVIRE, 2268 & ONE,DIA1(KOFFDIA),NVIRA) 2269 2270 CALL DGEMM('N','T',NVIR(ISYMA),NRHF(ISYMI),NRHF(ISYMJ), 2271 & ONE,WORK(KOFFTAJ),NVIRA,WORK(KOFFDIJ),NRHFI, 2272 & ONE,DIA2(KOFFDIA),NVIRA) 2273 2274 END DO 2275 2276 RETURN 2277 END 2278*---------------------------------------------------------------------* 2279* END OF SUBROUTINE CCSDT_FA_DEN1 * 2280*---------------------------------------------------------------------* 2281! -- end of cc_famat.F -- 2282