1! 2! Dalton, a molecular electronic structure program 3! Copyright (C) by the authors of Dalton. 4! 5! This program is free software; you can redistribute it and/or 6! modify it under the terms of the GNU Lesser General Public 7! License version 2.1 as published by the Free Software Foundation. 8! 9! This program is distributed in the hope that it will be useful, 10! but WITHOUT ANY WARRANTY; without even the implied warranty of 11! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 12! Lesser General Public License for more details. 13! 14! If a copy of the GNU LGPL v2.1 was not distributed with this 15! code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html. 16! 17! 18C 19*---------------------------------------------------------------------* 20c/* Deck CC_HMAT */ 21*=====================================================================* 22 SUBROUTINE CC_HMAT( LISTL, ! inp: Zeta vector list 23 & LISTA, ! inp: A amplitude list 24 & LISTB, ! inp: B amplitude list 25 & LISTC, ! inp: C amplitude list 26 & LISTD, ! inp: D amplitude list 27 & NHTRAN, ! inp: nb. of H matrix transf. 28 & MXVEC, ! inp: max. nb. of dot products 29 & IHTRAN, ! inp: index array of H mat. transf. 30 & IHDOTS, ! inp: index array of dot products 31 & HCON, ! out: array for dot products 32 & WORK, ! inp: work space 33 & LWORK, ! scr: length of work space 34 & IOPTRES)! inp: output option 35*---------------------------------------------------------------------* 36* 37* Purpose: calculation of linear transformations of 38* 3 amplitude vectors, T^A, T^B and T^C with the 39* H matrix (third partial derivative of the lagrangian 40* with respect to T) for index combinations for Lambda, 41* T^A, T^B, T^C passed in IHTRAN. 42* 43* epsilon_mu = < Lambda | [[[[H,T^A],T^B],T^C],tau_mu] | HF > 44* 45* return of the result vectors: 46* 47* IOPTRES = 0 : note used (but compare with CC_BMAT to see 48* for which purpose it is reserved) 49* 50* IOPTRES = 1 : the vectors are kept and returned in WORK 51* if possible, start addresses returned in 52* IHTRAN(5,*). N.B.: if WORK is not large 53* enough, CC_HMAT will stop!! 54* 55* IOPTRES = 3 : each result vector is written to its own 56* file by a call to CC_WRRSP, LISTD is used 57* as list type and IHTRAN(5,*) as list index 58* NOTE that IHTRAN(5,*) is in this case input! 59* 60* IOPTRES = 4 : each result vector is added to a vector on 61* file by a call to CC_WARSP, LISTD is used 62* as list type and IHTRAN(5,*) as list index 63* NOTE that IHTRAN(5,*) is in this case input! 64* 65* IOPTRES = 5 : the result vectors are dotted on a array 66* of vectors, the type of the arrays given 67* by LISTD and the indeces from IHDOTS 68* the result of the dot products is returned 69* in the HCON array 70* 71* 72* symmetries/variables: 73* 74* EPSI1, ISYRES : H matrix transformation 75* CTR2, ISYCTR : response lagrangian multipliers 76* T1AMPA, ISYMTA : A response amplitudes 77* T1AMPB, ISYMTB : B response amplitudes 78* T1AMPC, ISYMTC : C response amplitudes 79* T1AMPD, ISYMTD : D response amplitudes 80* 81* uses approximately V^2 O^2 + O^4 + O^3 work space 82* CC2: (ia|jb), CTR2 83* CCSD: (ia|jb), CTR2 84* 85* Written by Christof Haettig, Februar 1998. 86* CC3 noddy version, April 2002, Christof Haettig 87* 88*=====================================================================* 89#if defined (IMPLICIT_NONE) 90 IMPLICIT NONE 91#else 92# include "implicit.h" 93#endif 94#include "priunit.h" 95#include "ccsdsym.h" 96#include "ccsdinp.h" 97#include "ccorb.h" 98#include "ccnoddy.h" 99 100* local parameters: 101 LOGICAL LOCDBG 102 PARAMETER (LOCDBG = .FALSE.) 103 CHARACTER MSGDBG*(17) 104 PARAMETER (MSGDBG='[debug] CC_HMAT> ') 105 106#if defined (SYS_CRAY) 107 REAL ZERO, ONE, TWO, THREE 108#else 109 DOUBLE PRECISION ZERO, ONE, TWO, THREE 110#endif 111 PARAMETER (ZERO = 0.0d0, ONE = 1.0d0, TWO = 2.0d0, THREE = 3.0d0) 112 113 INTEGER KDUM 114 PARAMETER (KDUM = +99 999 999) ! dummy address 115 INTEGER ISYOVOV 116 PARAMETER (ISYOVOV = 1) 117 118 119 CHARACTER*(*) LISTL, LISTA, LISTB, LISTC, LISTD 120 INTEGER LWORK, IOPTRES, MXVEC, NHTRAN 121 INTEGER IHTRAN(5,NHTRAN) 122 INTEGER IHDOTS(MXVEC,NHTRAN) 123 124#if defined (SYS_CRAY) 125 REAL WORK(LWORK) 126 REAL HCON(MXVEC,NHTRAN) 127 REAL DDOT 128#else 129 DOUBLE PRECISION WORK(LWORK) 130 DOUBLE PRECISION HCON(MXVEC,NHTRAN) 131 DOUBLE PRECISION DDOT 132#endif 133 134 CHARACTER MODEL*(10), MODELW*(10) 135 LOGICAL LLSAME 136 INTEGER IZETAV, ITAMPA, ITAMPB, ITAMPC, ITAMPD, ITRAN, IFILE 137 INTEGER KT1AMPA, KT1AMPB, KT1AMPC, KT1AMPD, KEPSI1, KSTART 138 INTEGER ISYRES, ISYCTR, ISYMTA, ISYMTB, ISYMTC, ISYABCD 139 INTEGER KXIAJB, KCTR2, KEND1, LEND1, IOPT, IOPTW, KEND0, LEND0, 140 & KEPSI2, ILEN, IDBL 141 142 INTEGER ILSTSYM 143 144 145*---------------------------------------------------------------------* 146* begin: 147*---------------------------------------------------------------------* 148 149* check coupled cluster model: 150 IF (CCS) THEN 151 MODELW = 'CCS ' 152 IOPTW = 1 153 ELSE IF (CC2) THEN 154 MODELW = 'CC2 ' 155 IOPTW = 1 156 ELSE IF (CCSD) THEN 157 MODELW = 'CCSD ' 158 IOPTW = 1 159 ELSE IF (CC3) THEN 160 MODELW = 'CC3 ' 161 IOPTW = 3 162 ELSE 163 CALL QUIT('Unknown CC model in H matrix transformation') 164 END IF 165 166 IF ( .not. (CCS .or. CC2 .or. CCSD .or. CC3) ) THEN 167 WRITE(LUPRI,'(/a)') ' CC_HMAT called for a Coupled Cluster ' 168 & //'method not implemented in CC_HMAT...' 169 CALL QUIT('Unknown CC method in CC_HMAT.') 170 END IF 171 172* check list combination: 173 IF (LISTA(1:1).NE.'R' 174 & .OR. LISTB(1:1).NE.'R' 175 & .OR. LISTC(1:1).NE.'R' 176 & .OR. .NOT.(LISTD(1:1).EQ.'R' .OR. LISTD(1:1).EQ.'X') 177 & .OR. LISTL(1:1).NE.'L' ) THEN 178 WRITE(LUPRI,'(2A,/A)') 179 & ' In CC_HMAT LISTA, LISTB, LISTC, LISTD', 180 & ' should refer to t-amplitude vectors and LISTL to a ', 181 & ' multiplier vector.' 182 WRITE(LUPRI,'(2A)') 183 & ' LISTL: ',LISTL(1:3), 184 & ' LISTA: ',LISTA(1:3), 185 & ' LISTB: ',LISTB(1:3), 186 & ' LISTC: ',LISTC(1:3), 187 & ' LISTD: ',LISTD(1:3) 188 CALL QUIT('Strange LIST combination in CC_HMAT.') 189 END IF 190 191* check IOPTRES and initialize output array HCON, if used: 192 IF (IOPTRES.EQ.1 .OR. IOPTRES.EQ.3 .OR. IOPTRES.EQ.4 ) THEN 193 CONTINUE 194 ELSE IF (IOPTRES.EQ.5) THEN 195 IF (MXVEC*NHTRAN .GT. 0) CALL DZERO(HCON,MXVEC*NHTRAN) 196 ELSE 197 CALL QUIT('Illegal value of IOPTRES in CC_HMAT.') 198 END IF 199 200* in debug mode print length of work space: 201 IF (LOCDBG) THEN 202 WRITE(LUPRI,'(/1x,a,i15)') 'work space in CC_HMAT:',LWORK 203 Call FLSHFO(LUPRI) 204 END IF 205 206* flush print unit 207 Call FLSHFO(LUPRI) 208 209* initialize start address on the work space: 210 KSTART = 1 211 212*=====================================================================* 213* start loop over all requeste H matrix transformations: 214*=====================================================================* 215 DO ITRAN = 1, NHTRAN 216 217 IZETAV = IHTRAN(1,ITRAN) 218 ITAMPA = IHTRAN(2,ITRAN) 219 ITAMPB = IHTRAN(3,ITRAN) 220 ITAMPC = IHTRAN(4,ITRAN) 221 ITAMPD = IHTRAN(5,ITRAN) 222 IFILE = ITAMPD 223 224* set & check symmetries: 225 ISYCTR = ILSTSYM(LISTL,IZETAV) 226 ISYMTA = ILSTSYM(LISTA,ITAMPA) 227 ISYMTB = ILSTSYM(LISTB,ITAMPB) 228 ISYMTC = ILSTSYM(LISTC,ITAMPC) 229 230 ISYRES = MULD2H(MULD2H(ISYMTA,ISYMTB),MULD2H(ISYCTR,ISYMTC)) 231 232* allocated & initialize single excitation part of result vector EPSI1: 233 KEPSI1 = KSTART 234 KEND0 = KEPSI1 + NT1AM(ISYRES) 235 LEND0 = LWORK - KEND0 236 237 IF (LEND0 .LT. 0) THEN 238 CALL QUIT('Insufficient work space in CC_HMAT. (a)') 239 END IF 240 241 Call DZERO (WORK(KEPSI1), NT1AM(ISYRES)) 242 243*---------------------------------------------------------------------* 244* calculate H matrix transformation (for CCS the result is zero!): 245*---------------------------------------------------------------------* 246 IF (.NOT. CCS) THEN 247 248*---------------------------------------------------------------------* 249* initialize pointer for work space & read single parts of the vectors 250*---------------------------------------------------------------------* 251 KT1AMPA = KEND0 252 KT1AMPB = KT1AMPA + NT1AM(ISYMTA) 253 KT1AMPC = KT1AMPB + NT1AM(ISYMTB) 254 KCTR2 = KT1AMPC + NT1AM(ISYMTC) 255 KXIAJB = KCTR2 + NT2AM(ISYCTR) 256 KEND1 = KXIAJB + NT2AM(ISYOVOV) 257 LEND1 = LWORK - KEND1 258 259 IF (LEND1 .LT. 0) THEN 260 CALL QUIT('Insufficient work space in CC_HMAT. (0)') 261 END IF 262 263* * read packed (ia|jb) integrals: 264 Call CCG_RDIAJB(WORK(KXIAJB),NT2AM(ISYOVOV)) 265 266* read packed multipliers: 267 IOPT = 2 268 Call CC_RDRSP(LISTL,IZETAV,ISYCTR,IOPT,MODEL, 269 & WORK(KDUM),WORK(KCTR2)) 270 271* read A response amplitudes: 272 IOPT = 1 273 Call CC_RDRSP(LISTA,ITAMPA,ISYMTA,IOPT,MODEL, 274 & WORK(KT1AMPA),WORK(KDUM)) 275 276* read B response amplitudes: 277 IOPT = 1 278 Call CC_RDRSP(LISTB,ITAMPB,ISYMTB,IOPT,MODEL, 279 & WORK(KT1AMPB),WORK(KDUM)) 280 281* read C response amplitudes: 282 IOPT = 1 283 Call CC_RDRSP(LISTC,ITAMPC,ISYMTC,IOPT,MODEL, 284 & WORK(KT1AMPC),WORK(KDUM)) 285 286 IF (LOCDBG) THEN 287 Call AROUND('debug_CC_HMAT> response T1 amplitudes B:') 288 WRITE (LUPRI,*) 'List, Index, Symmetry:',LISTA, ITAMPA, ISYMTA 289 Call CC_PRP(WORK(KT1AMPA),WORK(KDUM),ISYMTA,1,0) 290 291 Call AROUND('debug_CC_HMAT> response T1 amplitudes C:') 292 WRITE (LUPRI,*) 'List, Index, Symmetry:',LISTB, ITAMPB, ISYMTB 293 Call CC_PRP(WORK(KT1AMPB),WORK(KDUM),ISYMTB,1,0) 294 295 Call AROUND('debug_CC_HMAT> response T1 amplitudes D:') 296 WRITE (LUPRI,*) 'List, Index, Symmetry:',LISTC, ITAMPC, ISYMTC 297 Call CC_PRP(WORK(KT1AMPC),WORK(KDUM),ISYMTC,1,0) 298 299 Call FLSHFO(LUPRI) 300 END IF 301 302*---------------------------------------------------------------------* 303* calculate C and D terms for the cyclic permutations of B, C, D: 304* (permutation of the first two vectors is treated inside CC_HCD) 305*---------------------------------------------------------------------* 306 CALL CC_HCD( WORK(KXIAJB), ISYOVOV,!inp: (ia|jb) integrals 307 & WORK(KCTR2), ISYCTR, !inp: double of Zeta vector 308 & WORK(KT1AMPA), ISYMTA, !inp: T1 amplitudes for B 309 & WORK(KT1AMPB), ISYMTB, !inp: T1 amplitudes for C 310 & WORK(KT1AMPC), ISYMTC, !inp: T1 amplitudes for D 311 & WORK(KEPSI1), ISYRES, !out: Epsilon result vector 312 & WORK(KEND1), LEND1 )!wrk: work space 313 314 CALL CC_HCD( WORK(KXIAJB), ISYOVOV,!inp: (ia|jb) integrals 315 & WORK(KCTR2), ISYCTR, !inp: double of Zeta vector 316 & WORK(KT1AMPB), ISYMTB, !inp: T1 amplitudes for C 317 & WORK(KT1AMPC), ISYMTC, !inp: T1 amplitudes for D 318 & WORK(KT1AMPA), ISYMTA, !inp: T1 amplitudes for B 319 & WORK(KEPSI1), ISYRES, !out: Epsilon result vector 320 & WORK(KEND1), LEND1 )!wrk: work space 321 322 CALL CC_HCD( WORK(KXIAJB), ISYOVOV,!inp: (ia|jb) integrals 323 & WORK(KCTR2), ISYCTR, !inp: double of Zeta vector 324 & WORK(KT1AMPC), ISYMTC, !inp: T1 amplitudes for D 325 & WORK(KT1AMPA), ISYMTA, !inp: T1 amplitudes for B 326 & WORK(KT1AMPB), ISYMTB, !inp: T1 amplitudes for C 327 & WORK(KEPSI1), ISYRES, !out: Epsilon result vector 328 & WORK(KEND1), LEND1 )!wrk: work space 329 330 END IF ! (.NOT. CCS) 331 332*---------------------------------------------------------------------* 333* add CC3 contributions: 334*---------------------------------------------------------------------* 335 IF (CCSDT) THEN 336 ! allocate and initialize doubles result vector: 337 KEPSI1 = KSTART 338 KEPSI2 = KEPSI1 + NT1AM(ISYRES) 339 KEND0 = KEPSI2 + NT2AM(ISYRES) 340 LEND0 = LWORK - KEND0 341 342 IF (LEND0 .LT. 0) THEN 343 CALL QUIT('Insufficient work space in CC_HMAT. (ccsdt)') 344 END IF 345 346 Call DZERO (WORK(KEPSI2), NT2AM(ISYRES)) 347 348 IF (NODDY_HMAT) THEN 349 350 CALL CCSDT_HMAT_NODDY(LISTL,IZETAV,LISTA,ITAMPA, 351 & LISTB,ITAMPB,LISTC,ITAMPC, 352 & WORK(KEPSI1),WORK(KEPSI2), 353 & WORK(KEND0),LEND0) 354 355 ELSE 356 357 CALL CC3_HMAT(LISTL,IZETAV,LISTA,ITAMPA, 358 & LISTB,ITAMPB,LISTC,ITAMPC, 359 & WORK(KEPSI1),WORK(KEPSI2),ISYRES, 360 & WORK(KEND0),LEND0) 361 362 END IF 363 364 ELSE 365 KEPSI2 = KSTART 366C ! to not get undefined address in calls below /hjaaj Aug 07 367 END IF 368*=====================================================================* 369* store result vectors: 370*=====================================================================* 371 IF ( IOPTRES .EQ. 1 ) THEN 372 373* output returned in work space: just increase start address 374 375 IHTRAN(5,ITRAN) = KSTART 376 KSTART = KSTART + NT1AM(ISYRES) 377 378 ELSE IF (IOPTRES .EQ. 3) THEN 379 380* write to file using CC_WRRSP, pass LISTD as list and 381* IFILE as index 382 383 CALL CC_WRRSP(LISTD,IFILE,ISYRES,IOPTW,MODELW,WORK(KDUM), 384 & WORK(KEPSI1),WORK(KEPSI2),WORK(KEND0),LEND0) 385 386 ELSE IF (IOPTRES .EQ. 4) THEN 387 388* add to vector on file using CC_WRRSP, pass LISTD as list 389* and IFILE as index 390 391 CALL CC_WARSP(LISTD,IFILE,ISYRES,IOPTW,MODELW,WORK(KDUM), 392 & WORK(KEPSI1),WORK(KEPSI2),WORK(KEND0),LEND0) 393 394 ELSE IF (IOPTRES .EQ. 5) THEN 395 396* calculate list of dot products 397 398 IF (IOPTW.GE.2) CALL CCLR_DIASCL(WORK(KEPSI2),TWO,ISYRES) 399 CALL CCDOTRSP(IHDOTS,HCON,IOPTW,LISTD,ITRAN,NHTRAN,MXVEC, 400 & WORK(KEPSI1),WORK(KEPSI2),ISYRES, 401 & WORK(KEND0),LEND0) 402 403 ELSE 404 CALL QUIT('Illegal value for IOPTRES in CC_HMAT.') 405 END IF 406 407*---------------------------------------------------------------------* 408* print debug information: 409*---------------------------------------------------------------------* 410 IF (LOCDBG) THEN 411 WRITE (LUPRI,*) MSGDBG, ' Lambda vector : ',LISTL,IZETAV 412 WRITE (LUPRI,*) MSGDBG, ' A amplitude vector: ',LISTA,ITAMPA 413 WRITE (LUPRI,*) MSGDBG, ' B amplitude vector: ',LISTB,ITAMPB 414 WRITE (LUPRI,*) MSGDBG, ' C amplitude vector: ',LISTC,ITAMPC 415 WRITE (LUPRI,*) MSGDBG, ' epsilon vector:' 416 ILEN = NT1AM(ISYRES) 417 IDBL = 0 418 IF (IOPTW.GE.2) THEN 419 ILEN = ILEN + NT2AM(ISYRES) 420 IDBL = 1 421 END IF 422 Call CC_PRP(WORK(KEPSI1),WORK(KEPSI2),ISYRES,1,IDBL) 423 WRITE (LUPRI,*) MSGDBG, ' Norm^2 of EPSI1:', 424 & DSQRT(DDOT(ILEN,WORK(KEPSI1),1,WORK(KEPSI1),1)) 425 Call FLSHFO(LUPRI) 426 END IF 427 428*---------------------------------------------------------------------* 429 END DO ! ITRAN 430 431 RETURN 432 END 433*=====================================================================* 434* END OF SUBROUTINE CC_HMAT 435*=====================================================================* 436 437*---------------------------------------------------------------------* 438c/* Deck CC_HCD */ 439*=====================================================================* 440 SUBROUTINE CC_HCD( XIAJB, ISYOVOV, ! inp: (ia|jb) integrals 441 & CTR2, ISYCTR, ! inp: double of Zeta vector 442 & T1AMPB, ISYMTB, ! inp: T1 amplitudes for B 443 & T1AMPC, ISYMTC, ! inp: T1 amplitudes for C 444 & T1AMPD, ISYMTD, ! inp: T1 amplitudes for D 445 & EPSI1, ISYRES, ! out: Epsilon result vector 446 & WORK, LWORK )! wrk: work space 447*---------------------------------------------------------------------* 448* 449* Purpose: calculate C and D term contributions for H matrix 450* for a particular permutation of the amplitude vectors 451* 452* Written by Christof Haettig, Februar 1998. 453*---------------------------------------------------------------------* 454#if defined (IMPLICIT_NONE) 455 IMPLICIT NONE 456#else 457# include "implicit.h" 458#endif 459#include "priunit.h" 460#include "ccsdsym.h" 461#include "ccorb.h" 462 463* local parameters: 464 CHARACTER MSGDBG*(16) 465 PARAMETER (MSGDBG='[debug] CC_HCD> ') 466 LOGICAL LOCDBG 467 PARAMETER (LOCDBG = .FALSE.) 468 469 INTEGER ISYMTB, ISYMTC, ISYMTD, ISYCTR, ISYRES, ISYOVOV, LWORK 470 471#if defined (SYS_CRAY) 472 REAL T1AMPB(*) ! dimension (NT1AM(ISYMTB)) 473 REAL T1AMPC(*) ! dimension (NT1AM(ISYMTC)) 474 REAL T1AMPD(*) ! dimension (NT1AM(ISYMTD)) 475 REAL EPSI1(*) ! dimension (NT1AM(ISYRES)) 476 REAL XIAJB(*) ! dimension (NT2AM(ISYOVOV)) 477 REAL CTR2(*) ! dimension (NT2AM(ISYCTR)) 478 REAL WORK(LWORK) 479 REAL DDOT 480#else 481 DOUBLE PRECISION T1AMPB(*) ! dimension (NT1AM(ISYMTB)) 482 DOUBLE PRECISION T1AMPC(*) ! dimension (NT1AM(ISYMTC)) 483 DOUBLE PRECISION T1AMPD(*) ! dimension (NT1AM(ISYMTD)) 484 DOUBLE PRECISION EPSI1(*) ! dimension (NT1AM(ISYRES)) 485 DOUBLE PRECISION XIAJB(*) ! dimension (NT2AM(ISYOVOV)) 486 DOUBLE PRECISION CTR2(*) ! dimension (NT2AM(ISYCTR)) 487 DOUBLE PRECISION WORK(LWORK) 488 DOUBLE PRECISION DDOT 489#endif 490 491 INTEGER ISYMTBC, ISYMA, ISYMLJK, ISYC4O, ISYX4O, ISYMI 492 INTEGER KC4O, KX4O, KEND1, LEND1, NAI, KOFF, KXLJKA, KCLJKA, IOPT 493 494 495*---------------------------------------------------------------------* 496* begin: 497*---------------------------------------------------------------------* 498 ISYMTBC = MULD2H(ISYMTB,ISYMTC) 499 500*---------------------------------------------------------------------* 501* D term: 502*---------------------------------------------------------------------* 503 504* calculate double one-index transformed Zeta vector: 505 ISYC4O = MULD2H(ISYCTR,ISYMTBC) 506 507 KC4O = 1 508 KEND1 = KC4O + N3ORHF(ISYC4O) 509 LEND1 = LWORK - KEND1 510 511 IF ( LEND1 .LT. 0 ) THEN 512 CALL QUIT('Insufficient work space in CC_HCD. (1)') 513 END IF 514 515 IOPT = 3 516 Call CCG_4O(WORK(KC4O),ISYC4O, CTR2,ISYCTR, 517 & T1AMPB,ISYMTB,T1AMPC,ISYMTC, 518 & WORK(KEND1),LEND1,IOPT) 519 520 521* start loop over virtual index a: 522 DO ISYMA = 1, NSYM 523 524 ISYMI = MULD2H(ISYRES,ISYMA) 525 ISYMLJK = MULD2H(ISYOVOV,MULD2H(ISYMTD,ISYMA)) 526 527 KXLJKA = KEND1 528 If ( LEND1 .LT. NMAIJK(ISYMLJK) ) THEN 529 CALL QUIT('Insufficient work space in CC_HCD. (2)') 530 END IF 531 532 DO A = 1, NVIR(ISYMA) 533 534* calculate batch of (l j^D | k a) integrals with fixed a: 535 IOPT = 1 536 Call CCG_TRANS2(WORK(KXLJKA),ISYMLJK,XIAJB,ISYOVOV, 537 & T1AMPD,ISYMTD,A,ISYMA,IOPT) 538 539* contract (l j^D|k a) with [ Zeta(l^C j|k^B i) + Zeta(l^B j|k^C i) ]: 540 DO I = 1, NRHF(ISYMI) 541 NAI = IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I-1) + A 542 KOFF = I3ORHF(ISYMLJK,ISYMI) + NMAIJK(ISYMLJK)*(I-1) 543 EPSI1(NAI) = EPSI1(NAI) + 544 & DDOT(NMAIJK(ISYMLJK),WORK(KC4O+KOFF),1,WORK(KXLJKA),1) 545 END DO 546 547 END DO ! A 548 END DO ! ISYMA 549 550 551*---------------------------------------------------------------------* 552* D term: 553*---------------------------------------------------------------------* 554 ISYX4O = MULD2H(ISYOVOV,ISYMTBC) 555 556 KX4O = 1 557 KEND1 = KX4O + N3ORHF(ISYX4O) 558 LEND1 = LWORK - KEND1 559 560 IF ( LEND1 .LT. 0 ) THEN 561 CALL QUIT('Insufficient work space in CC_HCD. (3)') 562 END IF 563 564* calculate double one-index transformed integrals: 565 IOPT = 3 566 Call CCG_4O(WORK(KX4O),ISYX4O, XIAJB,ISYOVOV, 567 & T1AMPB,ISYMTB,T1AMPC,ISYMTC, 568 & WORK(KEND1),LEND1,IOPT) 569 570C WRITE (LUPRI,*) 'KX4O:' 571C WRITE (LUPRI,'(3X,I5,3X,D25.14)'), 572C & (J+1,WORK(KX4O+J),J=0,N3ORHF(ISYX4O)-1) 573 574* start loop over virtual index a: 575 DO ISYMA = 1, NSYM 576 577 ISYMI = MULD2H(ISYRES,ISYMA) 578 ISYMLJK = MULD2H(ISYCTR,MULD2H(ISYMTD,ISYMA)) 579 580 KCLJKA = KEND1 581 If ( LEND1 .LT. NMAIJK(ISYMLJK) ) THEN 582 CALL QUIT('Insufficient work space in CC_HCD. (4)') 583 END IF 584 585 DO A = 1, NVIR(ISYMA) 586 587 CALL DCOPY(NMAIJK(ISYMLJK),999.99d0,0,WORK(KCLJKA),1) 588* calculate batch of one-index transf. Zeta(l j^D | k a) with fixed a: 589 IOPT = 1 590 Call CCG_TRANS2(WORK(KCLJKA),ISYMLJK,CTR2,ISYCTR, 591 & T1AMPD,ISYMTD,A,ISYMA,IOPT) 592 593* contract Zeta(l j^D|k a) with [ (l^C j|k^B i) + (l^B j|k^C i) ]: 594 DO I = 1, NRHF(ISYMI) 595 NAI = IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I-1) + A 596 KOFF = I3ORHF(ISYMLJK,ISYMI) + NMAIJK(ISYMLJK)*(I-1) 597 EPSI1(NAI) = EPSI1(NAI) + 598 & DDOT(NMAIJK(ISYMLJK),WORK(KX4O+KOFF),1,WORK(KCLJKA),1) 599 600C WRITE (LUPRI,*) 'KX4O LJKA:' 601C WRITE (LUPRI,'(3X,I5,3X,2D25.14)'), 602C & (J+1,WORK(KX4O+KOFF+J),WORK(KCLJKA+J),J=0,NMAIJK(ISYMLJK)-1) 603 END DO 604 605 END DO ! A 606 END DO ! ISYMA 607 608 RETURN 609 END 610*=====================================================================* 611* END OF SUBROUTINE CC_HCD * 612*=====================================================================* 613 614*=====================================================================* 615 SUBROUTINE CC_HTST(WORK,LWORK) 616*=====================================================================* 617C 618C perform finite difference test for H matrix: 619C 620C---------------------------------------------------------------------- 621#if defined (IMPLICIT_NONE) 622 IMPLICIT NONE 623#else 624# include "implicit.h" 625#endif 626#include "priunit.h" 627#include "ccsdinp.h" 628#include "ccorb.h" 629#include "ccsdsym.h" 630#include "ccqrinf.h" 631 632 INTEGER MXHTRAN, MXVEC 633 PARAMETER ( MXHTRAN = 5, MXVEC = 1) 634 635 CHARACTER*(3) LISTB, LISTC, LISTD, LISTL, LISTA 636 CHARACTER*(10) MODEL 637 INTEGER IDLSTB, IDLSTC, IDLSTD, IDLSTL 638 INTEGER KT1AMPB, KT2AMPB, KEND, LEND, LWORK 639 INTEGER ISYMB, ISYMC, ISYMD, ISYML, ISYRES 640 INTEGER KRESLT1, KRESLT2, IOPT, KEPSI1, KEPSI1A, IOPTRES 641 INTEGER IHTRAN(5,MXHTRAN), NHTRAN 642 INTEGER IHDOTS(MXVEC,MXHTRAN) 643 644#if defined (SYS_CRAY) 645 REAL WORK(LWORK) 646 REAL HCON(MXVEC,MXHTRAN) 647 REAL DDOT 648#else 649 DOUBLE PRECISION WORK(LWORK) 650 DOUBLE PRECISION HCON(MXVEC,MXHTRAN) 651 DOUBLE PRECISION DDOT 652#endif 653 654 INTEGER ILSTSYM 655 INTEGER IR1TAMP 656 657*---------------------------------------------------------------------* 658* call G matrix transformation 659*---------------------------------------------------------------------* 660 LISTL = 'L0' 661 LISTA = 'R1' 662 LISTB = 'R1' 663 LISTC = 'R1' 664 LISTD = 'R1' 665 666 IDLSTL = 0 667 IDLSTB = IR1TAMP('ZDIPLEN ',.FALSE.,0.0D0,1) 668 IDLSTC = IR1TAMP('ZDIPLEN ',.FALSE.,0.0D0,1) 669 IDLSTD = IR1TAMP('XDIPLEN ',.FALSE.,0.0D0,1) 670 671 IHTRAN(1,1) = IDLSTL 672 IHTRAN(2,1) = IDLSTB 673 IHTRAN(3,1) = IDLSTC 674 IHTRAN(4,1) = IDLSTD 675 NHTRAN = 1 676 677 ISYML = ILSTSYM(LISTL,IDLSTL) 678 ISYMB = ILSTSYM(LISTB,IDLSTB) 679 ISYMC = ILSTSYM(LISTB,IDLSTC) 680 ISYMD = ILSTSYM(LISTD,IDLSTD) 681 ISYRES = MULD2H(MULD2H(ISYML,ISYMD),MULD2H(ISYMB,ISYMC)) 682 683 KEPSI1 = 1 684 KEPSI1A = KEPSI1 + NT1AM(ISYRES) 685 KEND = KEPSI1A + NT1AM(ISYRES) 686 LEND = LWORK - KEND 687 688C 689C old H matrix routine: 690C 691C WRITE (LUPRI,*) 'CC_HTST: CALL NOW CCCR_K' 692C CALL CCCR_K(LISTL,IDLSTL,LISTB,IDLSTB,LISTC,IDLSTC,LISTD,IDLSTD, 693C & WORK(KEPSI1),ISYRES,WORK(KEND),LEND) 694C WRITE (LUPRI,*) 'CC_HTST: RETURNED FROM CCCR_K' 695 696 WRITE (LUPRI,*) 'CC_HTST: CALL NOW CC_HMAT' 697 IOPTRES = 1 698 CALL CC_HMAT(LISTL,LISTB,LISTC,LISTD,LISTA,NHTRAN,MXVEC, 699 & IHTRAN,IHDOTS,HCON, 700 & WORK(KEPSI1A),LWORK-KEPSI1A,IOPTRES) 701 WRITE (LUPRI,*) 'CC_HTST: RETURNED FROM CC_HMAT' 702 703 KT1AMPB = KEND 704 KT2AMPB = KT1AMPB + NT1AM(ISYMB) 705 KRESLT1 = KT2AMPB + NT2AM(ISYMB) 706 KRESLT2 = KRESLT1 + NT1AM(ISYRES) 707 KEND = KRESLT2 + NT2AM(ISYRES) 708 LEND = LWORK - KEND 709 710 IF (LEND .LT. 0) THEN 711 CALL QUIT('Insufficient work space in CC_HTST.') 712 END IF 713 714 IOPT = 3 715 Call CC_RDRSP(LISTB,IDLSTB,ISYMB,IOPT,MODEL, 716 & WORK(KT1AMPB),WORK(KT2AMPB)) 717 718* zero doubles of B and/or C vector: 719C CALL DZERO(WORK(KT2AMPB),NT2AM(ISYMB)) 720 721 CALL CC_FDH(NT1AM(ISYMB),NT2AM(ISYMB), 722 > LISTC,IDLSTC,LISTD,IDLSTD, 723 > WORK(KT1AMPB), WORK(KRESLT1), 724 > WORK(KEND), LEND) 725 726 IF (CCS) CALL DZERO(WORK(KRESLT2),NT2AM(ISYRES)) 727 728 IF (.TRUE.) THEN 729 WRITE (LUPRI,*) 'finite difference Epsilon vector:' 730 Call CC_PRP(WORK(KRESLT1),WORK(KRESLT2),ISYRES,1,1) 731 WRITE (LUPRI,*) 'analytical Epsilon vector (CCCR_K):' 732 Call CC_PRP(WORK(KEPSI1),WORK,ISYRES,1,0) 733 WRITE (LUPRI,*) 'analytical Epsilon vector (CC_HMAT):' 734 Call CC_PRP(WORK(KEPSI1A),WORK,ISYRES,1,0) 735 END IF 736 737 Call DAXPY(NT1AM(ISYRES),-1.0d0,WORK(KRESLT1),1,WORK(KEPSI1),1) 738 739C 740C compare with result of old H matrix routine: 741C 742C WRITE (LUPRI,*) 'Norm of difference between analytical Epsilon ' 743C > // 'vector (CCCR_K) and the numerical result:' 744C WRITE (LUPRI,*) 'singles excitation part:', 745C > DSQRT(DDOT(NT1AM(ISYRES),WORK(KEPSI1),1,WORK(KEPSI1),1)) 746C WRITE (LUPRI,*) 'double excitation part: ', 747C > DSQRT(DDOT(NT2AM(ISYRES),WORK(KRESLT2),1,WORK(KRESLT2),1)) 748 749 750C 751C compare with result of new H matrix routine: 752C 753 Call DAXPY(NT1AM(ISYRES),-1.0d0,WORK(KRESLT1),1,WORK(KEPSI1A),1) 754 755 WRITE (LUPRI,*) 'Norm of difference between analytical Epsilon ' 756 > // 'vector (CC_HMAT) and the numerical result:' 757 WRITE (LUPRI,*) 'singles excitation part:', 758 > DSQRT(DDOT(NT1AM(ISYRES),WORK(KEPSI1A),1,WORK(KEPSI1A),1)) 759 WRITE (LUPRI,*) 'double excitation part: ', 760 > DSQRT(DDOT(NT2AM(ISYRES),WORK(KRESLT2),1,WORK(KRESLT2),1)) 761 CALL FLSHFO(LUPRI) 762 763 RETURN 764 END 765C---------------------------------------------------------------------- 766 SUBROUTINE CC_FDH(NC1VEC,NC2VEC,LISTB,ITAMPB,LISTC,ITAMPC, 767 & TZAM,RESULT,WORK,LWORK) 768C 769C---------------------------------------------------------------------- 770C Test routine for calculating the CC K*t*t*t vector by 771C finite difference on the H-matrix transformation. 772C C.Haettig and Ove Christiansen oktober 1996, februar 1997 773C 774C---------------------------------------------------------------------- 775C 776#include "implicit.h" 777#include "priunit.h" 778#include "dummy.h" 779#include "maxorb.h" 780#include "iratdef.h" 781#include "ccorb.h" 782#include "aovec.h" 783#include "ccsdinp.h" 784#include "cclr.h" 785#include "ccsdsym.h" 786#include "ccsdio.h" 787#include "leinf.h" 788C 789 DIMENSION WORK(LWORK),ITADR(2) 790 PARAMETER (XHALF = 0.5D00,XMTWO = -2.0D00, DELTA = 1.0D-07) 791 PARAMETER (ONE = 1.0D00, ZERO = 0.0D00 ) 792 CHARACTER*10 MODEL 793C 794 INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 795C 796 IF (CCR12) CALL QUIT('Finite-difference on H-matrix for CCR12 '// 797 & 'not adapted') 798C 799 IF (IPRINT.GT.5) THEN 800 CALL AROUND( 'IN CC_FDH : MAKING FINITE DIFF. CC H-Matrix') 801 ENDIF 802C 803C---------------------------- 804C Work space allocations. 805C---------------------------- 806C 807 ISYMTR = 1 808 ISYMOP = 1 809C 810 NTAMP = NT1AM(ISYMTR) + NT2AM(ISYMTR) 811 NTAMP2 = NTAMP*(NC1VEC + NC2VEC ) 812 KF = 1 813 KRHO1 = KF + NTAMP2 814 KRHO2 = KRHO1 + NT1AMX 815 KC1AM = KRHO2 + MAX(NT2AMX,NT2AM(ISYMTR)) 816 KC2AM = KC1AM + NT1AM(ISYMTR) 817 KEND1 = KC2AM 818 * + MAX(NT2AMX,NT2AM(ISYMTR),NT2AO(ISYMTR), 819 * 2*NT2ORT(ISYMTR)) 820 LWRK1 = LWORK - KEND1 821C 822 KRHO1D = KEND1 823 KRHO2D = KRHO1D + NT1AMX 824 KEND2 = KRHO2D 825 * + MAX(NT2AMX,NT2AM(ISYMTR),NT2AO(ISYMTR), 826 * 2*NT2ORT(ISYMTR)) 827 LWRK2 = LWORK - KEND1 828C 829 IF (IPRINT .GT. 100 ) THEN 830 WRITE(LUPRI,*) ' IN CC_FDH: KF = ',KF 831 WRITE(LUPRI,*) ' IN CC_FDH: KRHO1 = ',KRHO1 832 WRITE(LUPRI,*) ' IN CC_FDH: KRHO2 = ',KRHO2 833 WRITE(LUPRI,*) ' IN CC_FDH: KC1AM = ',KC1AM 834 WRITE(LUPRI,*) ' IN CC_FDH: KC2AM = ',KC2AM 835 WRITE(LUPRI,*) ' IN CC_FDH: KRHO1D = ',KRHO1D 836 WRITE(LUPRI,*) ' IN CC_FDH: KRHO2D = ',KRHO2D 837 WRITE(LUPRI,*) ' IN CC_FDH: KEND2 = ',KEND2 838 WRITE(LUPRI,*) ' IN CC_FDH: LWRK2 = ',LWRK2 839 ENDIF 840 IF (LWRK2.LT.0 ) THEN 841 WRITE(LUPRI,*) 'Too little work space in CC_FDH ' 842 WRITE(LUPRI,*) 'AVAILABLE: LWORK = ',LWORK 843 WRITE(LUPRI,*) 'NEEDED (AT LEAST) = ',KEND2 844 CALL QUIT('TOO LITTLE WORKSPACE IN CC_FDH ') 845 ENDIF 846 KF2 = KF + NC1VEC*NTAMP 847C 848C--------------------- 849C Initializations. 850C--------------------- 851C 852 CALL DZERO(WORK(KC1AM),NT1AMX) 853 CALL DZERO(WORK(KC2AM),NT2AMX) 854 CALL DZERO(WORK(KF),NTAMP2) 855 IF (ABS(DELTA) .GT. 1.0D-15 ) THEN 856 DELTAI = 1.0D00/DELTA 857 ELSE 858 DELTAI = 1 859 ENDIF 860 X11 = 0.0D00 861 X12 = 0.0D00 862 X21 = 0.0D00 863 X22 = 0.0D00 864 XNJ = 0.0D00 865C 866C------------------------------------------------ 867C Read the CC reference amplitudes From disk. 868C------------------------------------------------ 869C 870 IOPT = 3 871 CALL CC_RDRSP('R0',0,1,IOPT,MODEL,WORK(KC1AM),WORK(KC2AM)) 872C 873C---------------------------------------------- 874C Save the CC reference amplitudes on disk. 875C---------------------------------------------- 876C 877 LUTAM = -1 878 CALL GPOPEN(LUTAM,'TAM_SAV','UNKNOWN',' ','UNFORMATTED',IDUMMY, 879 * .FALSE.) 880 REWIND(LUTAM) 881 WRITE(LUTAM) (WORK(KC1AM + I -1 ), I = 1, NT1AMX) 882 WRITE(LUTAM) (WORK(KC2AM + I -1 ), I = 1, NT2AMX) 883 CALL GPCLOSE(LUTAM,'KEEP') 884C 885 IF (IPRINT.GT.125) THEN 886 RHO1N = DDOT(NT1AMX,WORK(KC1AM),1,WORK(KC1AM),1) 887 RHO2N = DDOT(NT2AMX,WORK(KC2AM),1,WORK(KC2AM),1) 888 WRITE(LUPRI,*) 'Norm of T1AM: ',RHO1N 889 WRITE(LUPRI,*) 'Norm of T2AM: ',RHO2N 890 CALL CC_PRP(WORK(KC1AM),WORK(KC2AM),1,1,1) 891 ENDIF 892 RSPIM = .TRUE. 893C 894C------------------------------------ 895C Calculate reference G*T*T vector. 896C------------------------------------ 897C 898 CALL QUIT('CC_FDH has to be fixed because of new G matrix.') 899C CALL CCQR_G('L0',0,LISTB,ITAMPB,LISTC,ITAMPC,WORK(KRHO1D), 900C & ISYMTR,WORK(KRHO2D),LWORK-KRHO2D) 901 902C 903C------------------------- 904C Zero out components. 905C------------------------- 906C 907 IF (LCOR .OR. LSEC) THEN 908C 909 CALL CC_CORE(WORK(KRHO1D),WORK(KRHO2D),ISYMTR) 910C 911 ENDIF 912C 913 IF (IPRINT.GT.2) THEN 914 RHO1N = DDOT(NT1AMX,WORK(KRHO1D),1,WORK(KRHO1D),1) 915 RHO2N = DDOT(NT2AMX,WORK(KRHO2D),1,WORK(KRHO2D),1) 916 WRITE(LUPRI,*) 'Norm of RHO1: ',RHO1N,'ref' 917 WRITE(LUPRI,*) 'Norm of RHO2: ',RHO2N,'ref' 918 ENDIF 919 IF (IPRINT.GT.125) THEN 920 CALL CC_PRP(WORK(KRHO1D),WORK(KRHO2D),1,1,1) 921 ENDIF 922C 923 CALL DCOPY(NT1AMX,WORK(KRHO1D),1,WORK(KRHO1),1) 924 CALL DCOPY(NT2AMX,WORK(KRHO2D),1,WORK(KRHO2),1) 925C 926C============================================= 927C Calculate H-matrix by finite difference. 928C============================================= 929C 930 DO 100 I = 1, NC1VEC 931 WRITE (LUPRI,*) 'singles index:',I 932C 933C---------------------------------------- 934C Add finite displadement to t and 935C calculate new intermediates. 936C---------------------------------------- 937C 938 LUTAM = -1 939 CALL GPOPEN(LUTAM,'TAM_SAV','UNKNOWN',' ','UNFORMATTED', 940 * IDUMMY,.FALSE.) 941 READ(LUTAM) (WORK(KC1AM + J -1 ) , J = 1, NT1AMX) 942 READ(LUTAM) (WORK(KC2AM + J -1 ) , J = 1, NT2AMX) 943 CALL GPCLOSE(LUTAM,'KEEP') 944C 945 TI = SECOND() 946 WORK(KC1AM +I -1) = WORK(KC1AM +I -1 ) + DELTA 947 IF (LCOR .OR. LSEC) THEN 948 CALL CC_CORE(WORK(KC1AM),WORK(KC2AM),ISYMTR) 949 ENDIF 950C 951 IOPT = 3 952 CALL CC_WRRSP('R0',0,1,IOPT,MODEL,DUMMY,WORK(KC1AM), 953 * WORK(KC2AM),WORK(KEND2),LWRK2) 954C 955 RSPIM = .TRUE. 956 CALL CCRHSN(WORK(KRHO1D),WORK(KRHO2D),WORK(KC1AM), 957 * WORK(KC2AM),WORK(KEND2),LWRK2,'XXX') 958C 959C----------------------------------------- 960C Get the CC response vector again. 961C----------------------------------------- 962C 963C CALL DCOPY(NTAMP,TXAM,1,WORK(KC1AM),1) 964C 965C--------------------------------------- 966C For Test zero part of L vector. 967C--------------------------------------- 968C 969C IF ( L1TST ) THEN 970C CALL DZERO(WORK(KC2AM),NT2AMX) 971C ENDIF 972C IF ( L2TST ) THEN 973C CALL DZERO(WORK(KC1AM),NT1AMX) 974C ENDIF 975C 976C------------------ 977C Transform. 978C------------------ 979C 980 CALL QUIT('CC_FDH has to be fixed because of new G matrix.') 981C CALL CCQR_G('L0',0,LISTB,ITAMPB,LISTC,ITAMPC,WORK(KRHO1D), 982C & ISYMTR,WORK(KRHO2D),LWORK-KRHO2D) 983C 984 IF (LCOR .OR. LSEC) THEN 985 CALL CC_CORE(WORK(KRHO1D),WORK(KRHO2D),ISYMTR) 986 ENDIF 987C 988 IF (IPRINT.GT.2) THEN 989 RHO1N = DDOT(NT1AMX,WORK(KRHO1D),1,WORK(KRHO1D),1) 990 RHO2N = DDOT(NT2AMX,WORK(KRHO2D),1,WORK(KRHO2D),1) 991 WRITE(LUPRI,*) 'Norm of RHO1: ',RHO1N,'ai=',I 992 WRITE(LUPRI,*) 'Norm of RHO2: ',RHO2N,'ai=',I 993 ENDIF 994 IF (IPRINT.GT.125) THEN 995 CALL CC_PRP(WORK(KRHO1D),WORK(KRHO2D),1,1,1) 996 ENDIF 997 CALL DAXPY(NT1AMX,-1.0D00,WORK(KRHO1),1,WORK(KRHO1D),1) 998 CALL DAXPY(NT2AMX,-1.0D00,WORK(KRHO2),1,WORK(KRHO2D),1) 999 CALL DSCAL(NT1AMX,DELTAI,WORK(KRHO1D),1) 1000 CALL DSCAL(NT2AMX,DELTAI,WORK(KRHO2D),1) 1001 CALL DCOPY(NT1AMX,WORK(KRHO1D),1, 1002 * WORK(KF+NTAMP*(I-1)),1) 1003 CALL DCOPY(NT2AMX,WORK(KRHO2D),1, 1004 * WORK(KF+NTAMP*(I-1)+NT1AMX),1) 1005 X11 = X11 + DDOT(NT1AMX,WORK(KRHO1D),1,WORK(KRHO1D),1) 1006 X21 = X21 + DDOT(NT2AMX,WORK(KRHO2D),1,WORK(KRHO2D),1) 1007C 1008 TI = SECOND() - TI 1009 IF (IPRINT.GT.5 ) THEN 1010 WRITE(LUPRI,*) ' ' 1011 WRITE(LUPRI,*) 'FDH ROW NR. ',I,' DONE IN ',TI,' SEC.' 1012 ENDIF 1013C 1014 100 CONTINUE 1015C 1016C---------------------------------------------------------------- 1017C Loop over T2 amplitudes. Take care of diagonal t2 elements 1018C is in a different convention in the energy code. 1019C Factor 1/2 from right , and factor 2 from left. 1020C---------------------------------------------------------------- 1021C 1022 DO 200 NAI = 1, NT1AMX 1023 DO 300 NBJ = 1, NAI 1024 I = INDEX(NAI,NBJ) 1025C 1026 IF (I.LE.NC2VEC) THEN 1027 1028 WRITE (LUPRI,*) 'doubles index:',I 1029C 1030C-------------------------------------------- 1031C Add finite displacement to t and 1032C calculate new intermediates. 1033C------------------------------------------- 1034C 1035 LUTAM = -1 1036 CALL GPOPEN(LUTAM,'TAM_SAV','UNKNOWN',' ','UNFORMATTED', 1037 * IDUMMY,.FALSE.) 1038 READ(LUTAM) (WORK(KC1AM + J -1 ) , J = 1, NT1AMX) 1039 READ(LUTAM) (WORK(KC2AM + J -1 ) , J = 1, NT2AMX) 1040 CALL GPCLOSE(LUTAM,'KEEP') 1041C 1042 TI = SECOND() 1043 DELT = DELTA 1044 IF (NAI.EQ.NBJ) DELT = 2*DELTA 1045 WORK(KC2AM + I -1) = WORK(KC2AM+I -1) + DELT 1046 IF (LCOR .OR. LSEC) THEN 1047 CALL CC_CORE(WORK(KC1AM),WORK(KC2AM),ISYMTR) 1048 ENDIF 1049C 1050 IOPT = 3 1051 CALL CC_WRRSP('R0',0,1,IOPT,MODEL,DUMMY,WORK(KC1AM), 1052 * WORK(KC2AM),WORK(KEND2),LWRK2) 1053C 1054 RSPIM = .TRUE. 1055 CALL CCRHSN(WORK(KRHO1D),WORK(KRHO2D),WORK(KC1AM), 1056 * WORK(KC2AM),WORK(KEND2),LWRK2,'XXX') 1057C 1058C------------------------------------------- 1059C Get the CC response vector again. 1060C------------------------------------------- 1061C 1062C CALL DCOPY(NTAMP,TXAM,1,WORK(KC1AM),1) 1063C 1064C----------------------------------------- 1065C For Test zero part of L vector. 1066C----------------------------------------- 1067C 1068C IF ( L1TST ) THEN 1069C CALL DZERO(WORK(KC2AM),NT2AMX) 1070C ENDIF 1071C IF ( L2TST ) THEN 1072C CALL DZERO(WORK(KC1AM),NT1AMX) 1073C ENDIF 1074C 1075C RHO1N = DDOT(NT1AMX,WORK(KC1AM),1,WORK(KC1AM),1) 1076C RHO2N = DDOT(NT2AMX,WORK(KC2AM),1,WORK(KC2AM),1) 1077C IF ( DEBUG ) THEN 1078C WRITE(LUPRI,*) 'Norm of L1AM-inp: ',RHO1N 1079C WRITE(LUPRI,*) 'Norm of L2AM-inp: ',RHO2N 1080C ENDIF 1081C 1082C-------------------- 1083C Transform. 1084C-------------------- 1085C 1086 CALL QUIT('CC_FDH has to be fixed because of new G matrix.') 1087C CALL CCQR_G('L0',0,LISTB,ITAMPB,LISTC,ITAMPC,WORK(KRHO1D), 1088C & ISYMTR,WORK(KRHO2D),LWORK-KRHO2D) 1089C 1090 IF (LCOR .OR. LSEC) THEN 1091 CALL CC_CORE(WORK(KRHO1D),WORK(KRHO2D),ISYMTR) 1092 ENDIF 1093C 1094 IF (IPRINT.GT.2) THEN 1095 RHO1N = DDOT(NT1AMX,WORK(KRHO1D),1,WORK(KRHO1D),1) 1096 RHO2N = DDOT(NT2AMX,WORK(KRHO2D),1,WORK(KRHO2D),1) 1097 WRITE(LUPRI,*) 'Norm of RHO1: ',RHO1N,'aibj=',I 1098 WRITE(LUPRI,*) 'Norm of RHO2: ',RHO2N,'aibj=',I 1099 ENDIF 1100 IF (IPRINT.GT.125) THEN 1101 CALL CC_PRP(WORK(KRHO1D),WORK(KRHO2D),1,1,1) 1102 ENDIF 1103C 1104 CALL DAXPY(NT1AMX,-1.0D00,WORK(KRHO1),1,WORK(KRHO1D),1) 1105 CALL DAXPY(NT2AMX,-1.0D00,WORK(KRHO2),1,WORK(KRHO2D),1) 1106 CALL DSCAL(NT1AMX,DELTAI,WORK(KRHO1D),1) 1107 CALL DSCAL(NT2AMX,DELTAI,WORK(KRHO2D),1) 1108 CALL DCOPY(NT1AMX,WORK(KRHO1D),1, 1109 * WORK(KF2+NTAMP*(I-1)),1) 1110 CALL DCOPY(NT2AMX,WORK(KRHO2D),1, 1111 * WORK(KF2+NTAMP*(I-1)+NT1AMX),1) 1112C 1113 X12 = X12 + DDOT(NT1AMX,WORK(KRHO1D),1,WORK(KRHO1D),1) 1114 X22 = X22 + DDOT(NT2AMX,WORK(KRHO2D),1,WORK(KRHO2D),1) 1115 TI = SECOND() - TI 1116 IF (IPRINT.GT.5 ) THEN 1117 WRITE(LUPRI,*) ' ' 1118 WRITE(LUPRI,*) 'FDG ROW NR. ',I+NT1AMX, 1119 * ' DONE IN ',TI,' SEC.' 1120 ENDIF 1121C 1122 ENDIF 1123C 1124 300 CONTINUE 1125 200 CONTINUE 1126C 1127 WRITE(LUPRI,*) ' ' 1128 WRITE(LUPRI,*) '** FINITE DIFF WITH DELTA ',DELTA, '**' 1129 WRITE(LUPRI,*) ' ' 1130 IF ((IPRINT .GT. 4)) THEN 1131 CALL AROUND( 'FINITE DIFF. CC K*Tx*Ty-Matrix - 11 & 21 PART ') 1132 CALL OUTPUT(WORK(KF),1,NTAMP,1,NC1VEC,NTAMP,NC1VEC,1,LUPRI) 1133 CALL AROUND( 'FINITE DIFF. CC K*Tx*Ty-Matrix - 12 & 22 PART ') 1134 CALL OUTPUT(WORK(KF+NTAMP*NC1VEC),1,NTAMP,1,NC2VEC, 1135 * NTAMP,NC2VEC,1,LUPRI) 1136 ENDIF 1137 IF (.TRUE.) THEN 1138 XNJ = X11 + X12 + X21 + X22 1139 WRITE(LUPRI,*) ' ' 1140 WRITE(LUPRI,*) ' NORM OF FIN. DIFF. K*tx*ty-Matrix.', SQRT(XNJ) 1141 WRITE(LUPRI,*) ' ' 1142 WRITE (LUPRI,*) ' NORM OF 11 PART OF FD. K*tx*ty-mat.: ', 1143 & SQRT(X11) 1144 WRITE (LUPRI,*) ' NORM OF 21 PART OF FD. K*tx*ty-mat.: ', 1145 & SQRT(X21) 1146 WRITE (LUPRI,*) ' NORM OF 12 PART OF FD. K*tx*ty-mat.: ', 1147 & SQRT(X12) 1148 WRITE (LUPRI,*) ' NORM OF 22 PART OF FD. K*tx*ty-mat.: ', 1149 & SQRT(X22) 1150 ENDIF 1151C 1152C-------------------------------------- 1153C Calculate Matrix times Tz vector. 1154C-------------------------------------- 1155C 1156 CALL DGEMV('N',NTAMP,NTAMP,ONE,WORK(KF),NTAMP,TZAM,1, 1157 * ZERO,RESULT,1) 1158C 1159C------------------------------------------------- 1160C Restore the CC reference amplitudes on disk. 1161C------------------------------------------------- 1162C 1163 LUTAM = -1 1164 CALL GPOPEN(LUTAM,'TAM_SAV','UNKNOWN',' ','UNFORMATTED',IDUMMY, 1165 * .FALSE.) 1166 REWIND(LUTAM) 1167 READ(LUTAM) (WORK(KC1AM + I -1 ) , I = 1, NT1AMX) 1168 READ(LUTAM) (WORK(KC2AM + I -1 ) , I = 1, NT2AMX) 1169 CALL GPCLOSE(LUTAM,'DELETE') 1170C 1171 IOPT = 3 1172 CALL CC_WRRSP('R0',0,1,IOPT,MODEL,DUMMY,WORK(KC1AM), 1173 * WORK(KC2AM),WORK(KEND2),LWRK2) 1174C 1175 RSPIM = .TRUE. 1176 CALL CCRHSN(WORK(KRHO1D),WORK(KRHO2D),WORK(KC1AM), 1177 & WORK(KC2AM), 1178 * WORK(KEND2),LWRK2,'XXX') 1179C 1180 IF (IPRINT .GT. 10) THEN 1181 CALL AROUND(' END OF CC_FDH ') 1182 ENDIF 1183C 1184 RETURN 1185 END 1186*=====================================================================* 1187*=====================================================================* 1188 SUBROUTINE CCSDT_HMAT_NODDY(LISTL,IDLSTL,LISTB,IDLSTB, 1189 & LISTC,IDLSTC,LISTD,IDLSTD, 1190 & OMEGA1,OMEGA2,WORK,LWORK) 1191*---------------------------------------------------------------------* 1192* 1193* Purpose: compute triples contribution to G matrix transformation 1194* 1195* (H T^B T^C T^D)^eff_1,2 = (H T^B T^C T^D)_1,2(CCSD) 1196* + (H T^B T^C)_1,2(L_3) 1197* 1198* Written by Christof Haettig, April 2002 1199* based on CCSDT_GMAT_NODDY 1200* 1201*=====================================================================* 1202#if defined (IMPLICIT_NONE) 1203 IMPLICIT NONE 1204#else 1205# include "implicit.h" 1206#endif 1207#include "priunit.h" 1208#include "ccsdinp.h" 1209#include "maxorb.h" 1210#include "ccsdsym.h" 1211#include "ccfield.h" 1212#include "ccorb.h" 1213#include "dummy.h" 1214 1215 LOGICAL LOCDBG 1216 PARAMETER (LOCDBG=.false.) 1217 1218 INTEGER ISYM0 1219 PARAMETER (ISYM0 = 1) 1220 1221 CHARACTER*3 LISTL, LISTB, LISTC, LISTD 1222 INTEGER LWORK, IDLSTL, IDLSTB, IDLSTC, IDLSTD 1223 1224#if defined (SYS_CRAY) 1225 REAL WORK(LWORK), TWO, FREQL, DDOT 1226 REAL OMEGA1(NT1AMX), OMEGA2(NT2AMX) 1227#else 1228 DOUBLE PRECISION WORK(LWORK), TWO, FREQL, DDOT 1229 DOUBLE PRECISION OMEGA1(NT1AMX), OMEGA2(NT2AMX) 1230#endif 1231 PARAMETER ( TWO = 2.0D0 ) 1232 1233 CHARACTER*10 MODEL 1234 INTEGER KXINT, KXIAJB, KYIAJB, KT1AMP0, KLAMP0, KLAMH0, KEND1, 1235 & KFOCK0, LWRK1, KLAMPB, KLAMHB, KLAMPC, KLAMHC, KT1AMPB, 1236 & KLAMPD, KLAMHD, KT1AMPD, KT1AMPC, 1237 & KINT1T0, KINT2T0, 1238 & KINT1SBCD, KINT2SBCD, 1239 & KBDIOVVO, KBDIOOVV, KBDIOOOO, KBDIVVVV, 1240 & KCDIOVVO, KCDIOOVV, KCDIOOOO, KCDIVVVV, 1241 & KBCIOVVO, KBCIOOVV, KBCIOOOO, KBCIVVVV, 1242 & KOME1, KOME2, KL1AM, KL2AM, KL3AM, KT2AM, KFOCKD, 1243 & KSCR1, KDUM 1244 INTEGER ISYMD, ILLL, IDEL, ISYDIS, IJ, NIJ, LUSIFC, INDEX, 1245 & ISYMC, ISYMB, LUFOCK, KEND2, LWRK2, IOPT, ILSTSYM, ISYML 1246 1247 INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 1248 1249 KDUM = 1 1250 CALL QENTER('CCSDT_HMAT_NODDY') 1251 1252 IF (DIRECT) CALL QUIT('CCSDT_HMAT_NODDY: DIRECT NOT IMPLEMENTED') 1253 1254*---------------------------------------------------------------------* 1255* Memory allocation: 1256*---------------------------------------------------------------------* 1257 KSCR1 = 1 1258 KFOCKD = KSCR1 + NT1AMX 1259 KEND1 = KFOCKD + NORBT 1260 1261 KFOCK0 = KEND1 1262 KEND1 = KFOCK0 + NORBT*NORBT 1263 1264 KT1AMP0 = KEND1 1265 KLAMP0 = KT1AMP0 + NT1AMX 1266 KLAMH0 = KLAMP0 + NLAMDT 1267 KEND1 = KLAMH0 + NLAMDT 1268 1269 KT1AMPB = KEND1 1270 KLAMPB = KT1AMPB + NT1AMX 1271 KLAMHB = KLAMPB + NLAMDT 1272 KEND1 = KLAMHB + NLAMDT 1273 1274 KT1AMPC = KEND1 1275 KLAMPC = KT1AMPC + NT1AMX 1276 KLAMHC = KLAMPC + NLAMDT 1277 KEND1 = KLAMHC + NLAMDT 1278 1279 KT1AMPD = KEND1 1280 KLAMPD = KT1AMPD + NT1AMX 1281 KLAMHD = KLAMPD + NLAMDT 1282 KEND1 = KLAMHD + NLAMDT 1283 1284 KXIAJB = KEND1 1285 KYIAJB = KXIAJB + NT1AMX*NT1AMX 1286 KEND1 = KYIAJB + NT1AMX*NT1AMX 1287 1288 KINT1T0 = KEND1 1289 KINT2T0 = KINT1T0 + NT1AMX*NVIRT*NVIRT 1290 KEND1 = KINT2T0 + NRHFT*NRHFT*NT1AMX 1291 1292 KBDIOVVO = KEND1 1293 KBDIOOVV = KBDIOVVO + NRHFT*NVIRT*NVIRT*NRHFT 1294 KBDIOOOO = KBDIOOVV + NRHFT*NVIRT*NVIRT*NRHFT 1295 KBDIVVVV = KBDIOOOO + NRHFT*NRHFT*NRHFT*NRHFT 1296 KEND1 = KBDIVVVV + NVIRT*NVIRT*NVIRT*NVIRT 1297 1298 KCDIOVVO = KEND1 1299 KCDIOOVV = KCDIOVVO + NRHFT*NVIRT*NVIRT*NRHFT 1300 KCDIOOOO = KCDIOOVV + NRHFT*NVIRT*NVIRT*NRHFT 1301 KCDIVVVV = KCDIOOOO + NRHFT*NRHFT*NRHFT*NRHFT 1302 KEND1 = KCDIVVVV + NVIRT*NVIRT*NVIRT*NVIRT 1303 1304 KBCIOVVO = KEND1 1305 KBCIOOVV = KBCIOVVO + NRHFT*NVIRT*NVIRT*NRHFT 1306 KBCIOOOO = KBCIOOVV + NRHFT*NVIRT*NVIRT*NRHFT 1307 KBCIVVVV = KBCIOOOO + NRHFT*NRHFT*NRHFT*NRHFT 1308 KEND1 = KBCIVVVV + NVIRT*NVIRT*NVIRT*NVIRT 1309 1310 KINT1SBCD = KEND1 1311 KINT2SBCD = KINT1SBCD + NT1AMX*NVIRT*NVIRT 1312 KEND1 = KINT2SBCD + NRHFT*NRHFT*NT1AMX 1313 1314 KOME1 = KEND1 1315 KOME2 = KOME1 + NT1AMX 1316 KEND1 = KOME2 + NT1AMX*NT1AMX 1317 1318 KL1AM = KEND1 1319 KL2AM = KL1AM + NT1AMX 1320 KL3AM = KL2AM + NT1AMX*NT1AMX 1321 KEND1 = KL3AM + NT1AMX*NT1AMX*NT1AMX 1322 1323 LWRK1 = LWORK - KEND1 1324 IF (LWRK1 .LT. 0) THEN 1325 CALL QUIT('Insufficient space in CCSDT_HMAT_NODDY') 1326 ENDIF 1327 1328*---------------------------------------------------------------------* 1329* Read SCF orbital energies from file: 1330*---------------------------------------------------------------------* 1331 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY, 1332 & .FALSE.) 1333 REWIND LUSIFC 1334 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 1335 READ (LUSIFC) 1336 READ (LUSIFC) (WORK(KFOCKD+I-1), I=1,NORBT) 1337 CALL GPCLOSE(LUSIFC,'KEEP') 1338 1339*---------------------------------------------------------------------* 1340* Get zeroth-order Lambda matrices: 1341*---------------------------------------------------------------------* 1342 IOPT = 1 1343 Call CC_RDRSP('R0',0,ISYM0,IOPT,MODEL,WORK(KT1AMP0),WORK(KDUM)) 1344 1345 Call LAMMAT(WORK(KLAMP0),WORK(KLAMH0),WORK(KT1AMP0), 1346 & WORK(KEND1),LWRK1) 1347 1348*---------------------------------------------------------------------* 1349* Get response Lambda matrices: 1350*---------------------------------------------------------------------* 1351 ISYMB = ILSTSYM(LISTB,IDLSTB) 1352 IOPT = 1 1353 Call CC_RDRSP(LISTB,IDLSTB,ISYMB,IOPT,MODEL, 1354 & WORK(KT1AMPB),WORK(KDUM)) 1355 1356 CALL CCLR_LAMTRA(WORK(KLAMP0),WORK(KLAMPB), 1357 & WORK(KLAMH0),WORK(KLAMHB),WORK(KT1AMPB),ISYMB) 1358 1359 ISYMC = ILSTSYM(LISTC,IDLSTC) 1360 IOPT = 1 1361 Call CC_RDRSP(LISTC,IDLSTC,ISYMC,IOPT,MODEL, 1362 & WORK(KT1AMPC),WORK(KDUM)) 1363 1364 CALL CCLR_LAMTRA(WORK(KLAMP0),WORK(KLAMPC), 1365 & WORK(KLAMH0),WORK(KLAMHC),WORK(KT1AMPC),ISYMC) 1366 1367 ISYMD = ILSTSYM(LISTD,IDLSTD) 1368 IOPT = 1 1369 Call CC_RDRSP(LISTD,IDLSTD,ISYMD,IOPT,MODEL, 1370 & WORK(KT1AMPD),WORK(KDUM)) 1371 1372 CALL CCLR_LAMTRA(WORK(KLAMP0),WORK(KLAMPD), 1373 & WORK(KLAMH0),WORK(KLAMHD),WORK(KT1AMPD),ISYMD) 1374 1375*---------------------------------------------------------------------* 1376* read zeroth-order AO Fock matrix from file: 1377*---------------------------------------------------------------------* 1378 LUFOCK = -1 1379 CALL GPOPEN(LUFOCK,'CC_FCKH','OLD',' ','UNFORMATTED',IDUMMY, 1380 & .FALSE.) 1381 REWIND(LUFOCK) 1382 READ(LUFOCK) (WORK(KFOCK0-1+I),I=1,N2BST(ISYM0)) 1383 CALL GPCLOSE(LUFOCK,'KEEP') 1384 1385 CALL CC_FCKMO(WORK(KFOCK0),WORK(KLAMP0),WORK(KLAMH0), 1386 & WORK(KEND1),LWRK1,ISYM0,ISYM0,ISYM0) 1387 1388*---------------------------------------------------------------------* 1389* Compute integrals needed for the following contributions: 1390*---------------------------------------------------------------------* 1391 1392 CALL DZERO(WORK(KBDIOVVO),NRHFT*NVIRT*NVIRT*NRHFT) 1393 CALL DZERO(WORK(KBDIOOVV),NRHFT*NVIRT*NVIRT*NRHFT) 1394 CALL DZERO(WORK(KBDIOOOO),NRHFT*NRHFT*NRHFT*NRHFT) 1395 CALL DZERO(WORK(KBDIVVVV),NVIRT*NVIRT*NVIRT*NVIRT) 1396 1397 CALL DZERO(WORK(KCDIOVVO),NRHFT*NVIRT*NVIRT*NRHFT) 1398 CALL DZERO(WORK(KCDIOOVV),NRHFT*NVIRT*NVIRT*NRHFT) 1399 CALL DZERO(WORK(KCDIOOOO),NRHFT*NRHFT*NRHFT*NRHFT) 1400 CALL DZERO(WORK(KCDIVVVV),NVIRT*NVIRT*NVIRT*NVIRT) 1401 1402 CALL DZERO(WORK(KBCIOVVO),NRHFT*NVIRT*NVIRT*NRHFT) 1403 CALL DZERO(WORK(KBCIOOVV),NRHFT*NVIRT*NVIRT*NRHFT) 1404 CALL DZERO(WORK(KBCIOOOO),NRHFT*NRHFT*NRHFT*NRHFT) 1405 CALL DZERO(WORK(KBCIVVVV),NVIRT*NVIRT*NVIRT*NVIRT) 1406 1407 CALL DZERO(WORK(KXIAJB),NT1AMX*NT1AMX) 1408 CALL DZERO(WORK(KYIAJB),NT1AMX*NT1AMX) 1409 1410 CALL DZERO(WORK(KINT1T0),NT1AMX*NVIRT*NVIRT) 1411 CALL DZERO(WORK(KINT2T0),NRHFT*NRHFT*NT1AMX) 1412 1413 CALL DZERO(WORK(KINT1SBCD),NT1AMX*NVIRT*NVIRT) 1414 CALL DZERO(WORK(KINT2SBCD),NRHFT*NRHFT*NT1AMX) 1415 1416 DO ISYMD = 1, NSYM 1417 DO ILLL = 1,NBAS(ISYMD) 1418 IDEL = IBAS(ISYMD) + ILLL 1419 ISYDIS = MULD2H(ISYMD,ISYMOP) 1420 1421C ---------------------------- 1422C Work space allocation no. 2. 1423C ---------------------------- 1424 KXINT = KEND1 1425 KEND2 = KXINT + NDISAO(ISYDIS) 1426 LWRK2 = LWORK - KEND2 1427 IF (LWRK2 .LT. 0) THEN 1428 WRITE(LUPRI,*) 'Need : ',KEND2,'Available : ',LWORK 1429 CALL QUIT('Insufficient space in CCSDT_HMAT_NODDY') 1430 ENDIF 1431 1432C --------------------------- 1433C Read in batch of integrals. 1434C --------------------------- 1435 CALL CCRDAO(WORK(KXINT),IDEL,1,WORK(KEND2),LWRK2, 1436 * WORK(KDUM),DIRECT) 1437 1438C ---------------------------------- 1439C Calculate integrals needed in CC3: 1440C ---------------------------------- 1441 CALL CCSDT_TRAN1(WORK(KINT1T0),WORK(KINT2T0), 1442 & WORK(KLAMP0),WORK(KLAMH0), 1443 & WORK(KXINT),IDEL) 1444 1445 CALL CC3_TRAN2(WORK(KXIAJB),WORK(KYIAJB), 1446 & WORK(KLAMP0),WORK(KLAMH0), 1447 & WORK(KXINT),IDEL) 1448 1449 ! XINT1S = XINT1S + (C-barB K-barC|B-barD D) 1450 ! XINT2S = XINT2S + (C-barB K-barC|L J-barD) 1451 CALL CCSDT_TRAN3_R(WORK(KINT1SBCD),WORK(KINT2SBCD), 1452 & WORK(KLAMP0),WORK(KLAMH0), 1453 & WORK(KLAMPB),WORK(KLAMHC), 1454 & WORK(KLAMPD),WORK(KLAMHD), 1455 & WORK(KXINT),IDEL) 1456 1457 ! XINT1S = XINT1S + (C-barB K-barD|B-barC D) 1458 ! XINT2S = XINT2S + (C-barB K-barD|L J-barC) 1459 CALL CCSDT_TRAN3_R(WORK(KINT1SBCD),WORK(KINT2SBCD), 1460 & WORK(KLAMP0),WORK(KLAMH0), 1461 & WORK(KLAMPB),WORK(KLAMHD), 1462 & WORK(KLAMPC),WORK(KLAMHC), 1463 & WORK(KXINT),IDEL) 1464 1465 ! XINT1S = XINT1S + (C-barC K-barB|B-barD D) 1466 ! XINT2S = XINT2S + (C-barC K-barB|L J-barD) 1467 CALL CCSDT_TRAN3_R(WORK(KINT1SBCD),WORK(KINT2SBCD), 1468 & WORK(KLAMP0),WORK(KLAMH0), 1469 & WORK(KLAMPC),WORK(KLAMHB), 1470 & WORK(KLAMPD),WORK(KLAMHD), 1471 & WORK(KXINT),IDEL) 1472 1473 ! XINT1S = XINT1S + (C-barD K-barB|B-barC D) 1474 ! XINT2S = XINT2S + (C-barD K-barB|L J-barC) 1475 CALL CCSDT_TRAN3_R(WORK(KINT1SBCD),WORK(KINT2SBCD), 1476 & WORK(KLAMP0),WORK(KLAMH0), 1477 & WORK(KLAMPD),WORK(KLAMHB), 1478 & WORK(KLAMPC),WORK(KLAMHC), 1479 & WORK(KXINT),IDEL) 1480 1481 ! XINT1S = XINT1S + (C-barC K-barD|B-barB D) 1482 ! XINT2S = XINT2S + (C-barC K-barD|L J-barB) 1483 CALL CCSDT_TRAN3_R(WORK(KINT1SBCD),WORK(KINT2SBCD), 1484 & WORK(KLAMP0),WORK(KLAMH0), 1485 & WORK(KLAMPC),WORK(KLAMHD), 1486 & WORK(KLAMPB),WORK(KLAMHB), 1487 & WORK(KXINT),IDEL) 1488 1489 ! XINT1S = XINT1S + (C-barD K-barC|B-barB D) 1490 ! XINT2S = XINT2S + (C-barD K-barC|L J-barB) 1491 CALL CCSDT_TRAN3_R(WORK(KINT1SBCD),WORK(KINT2SBCD), 1492 & WORK(KLAMP0),WORK(KLAMH0), 1493 & WORK(KLAMPD),WORK(KLAMHC), 1494 & WORK(KLAMPB),WORK(KLAMHB), 1495 & WORK(KXINT),IDEL) 1496 1497C ---------------------------------------------- 1498C (OV|VO)-BD (OO|VV)-BD (OO|OO)-BD (VV|VV)-BD 1499C ---------------------------------------------- 1500 CALL CCFOP_TRAN1_R(WORK(KBDIOVVO),WORK(KBDIOOVV), 1501 & WORK(KBDIOOOO),WORK(KBDIVVVV), 1502 & WORK(KLAMP0),WORK(KLAMH0), 1503 & WORK(KLAMPB),WORK(KLAMHB), 1504 & WORK(KLAMPD),WORK(KLAMHD), 1505 & WORK(KXINT),IDEL) 1506 1507 CALL CCFOP_TRAN1_R(WORK(KBDIOVVO),WORK(KBDIOOVV), 1508 & WORK(KBDIOOOO),WORK(KBDIVVVV), 1509 & WORK(KLAMP0),WORK(KLAMH0), 1510 & WORK(KLAMPD),WORK(KLAMHD), 1511 & WORK(KLAMPB),WORK(KLAMHB), 1512 & WORK(KXINT),IDEL) 1513 1514C ---------------------------------------------- 1515C (OV|VO)-CD (OO|VV)-CD (OO|OO)-CD (VV|VV)-CD 1516C ---------------------------------------------- 1517 CALL CCFOP_TRAN1_R(WORK(KCDIOVVO),WORK(KCDIOOVV), 1518 & WORK(KCDIOOOO),WORK(KCDIVVVV), 1519 & WORK(KLAMP0),WORK(KLAMH0), 1520 & WORK(KLAMPC),WORK(KLAMHC), 1521 & WORK(KLAMPD),WORK(KLAMHD), 1522 & WORK(KXINT),IDEL) 1523 1524 CALL CCFOP_TRAN1_R(WORK(KCDIOVVO),WORK(KCDIOOVV), 1525 & WORK(KCDIOOOO),WORK(KCDIVVVV), 1526 & WORK(KLAMP0),WORK(KLAMH0), 1527 & WORK(KLAMPD),WORK(KLAMHD), 1528 & WORK(KLAMPC),WORK(KLAMHC), 1529 & WORK(KXINT),IDEL) 1530 1531C ---------------------------------------------- 1532C (OV|VO)-BC (OO|VV)-BC (OO|OO)-BC (VV|VV)-BC 1533C ---------------------------------------------- 1534 CALL CCFOP_TRAN1_R(WORK(KBCIOVVO),WORK(KBCIOOVV), 1535 & WORK(KBCIOOOO),WORK(KBCIVVVV), 1536 & WORK(KLAMP0),WORK(KLAMH0), 1537 & WORK(KLAMPB),WORK(KLAMHB), 1538 & WORK(KLAMPC),WORK(KLAMHC), 1539 & WORK(KXINT),IDEL) 1540 1541 CALL CCFOP_TRAN1_R(WORK(KBCIOVVO),WORK(KBCIOOVV), 1542 & WORK(KBCIOOOO),WORK(KBCIVVVV), 1543 & WORK(KLAMP0),WORK(KLAMH0), 1544 & WORK(KLAMPC),WORK(KLAMHC), 1545 & WORK(KLAMPB),WORK(KLAMHB), 1546 & WORK(KXINT),IDEL) 1547 1548 END DO 1549 END DO 1550 1551*---------------------------------------------------------------------* 1552* Compute L^0_3 multipliers: 1553*---------------------------------------------------------------------* 1554 ISYML = ILSTSYM(LISTL,IDLSTL) 1555 1556 CALL DZERO(WORK(KL3AM),NT1AMX*NT1AMX*NT1AMX) 1557 1558 IF (LISTL(1:3).EQ.'L0 ') THEN 1559 1560 IF (LWRK1 .LT. NT2AMX) 1561 * CALL QUIT('Insufficient space in CCSDT_HMAT_NODDY') 1562 1563 ! read left amplitudes from file and square up the doubles part 1564 IOPT = 3 1565 Call CC_RDRSP(LISTL,IDLSTL,ISYML,IOPT,MODEL, 1566 * WORK(KL1AM),WORK(KEND1)) 1567 CALL CC_T2SQ(WORK(KEND1),WORK(KL2AM),ISYML) 1568 1569 CALL CCSDT_L03AM(WORK(KL3AM),WORK(KINT1T0),WORK(KINT2T0), 1570 * WORK(KXIAJB),WORK(KFOCK0),WORK(KL1AM), 1571 * WORK(KL2AM),WORK(KSCR1),WORK(KFOCKD), 1572 * DUMMY,DUMMY) 1573 1574 ELSE IF (LISTL(1:3).EQ.'L1 ' .OR. LISTL(1:3).EQ.'LE ' .OR. 1575 & LISTL(1:3).EQ.'M1 ' .OR. LISTL(1:3).EQ.'N2 ' .OR. 1576 & LISTL(1:3).EQ.'E0 ' ) THEN 1577 1578 CALL CCSDT_TBAR31_NODDY(WORK(KL3AM),FREQL,LISTL,IDLSTL, 1579 & WORK(KLAMP0),WORK(KLAMH0), 1580 & WORK(KFOCK0),WORK(KFOCKD),WORK(KSCR1), 1581 & WORK(KXIAJB),WORK(KINT1T0),WORK(KINT2T0), 1582 & WORK(KEND1),LWRK1) 1583 1584 ELSE 1585 1586 CALL QUIT('CCSDT_HMAT_NODDY> LISTL NOT AVAILABLE:'//LISTL) 1587 1588 END IF 1589 1590 CALL DSCAL(NT1AMX*NT1AMX*NT1AMX,-1.0D0,WORK(KL3AM),1) 1591 1592 IF (LOCDBG) THEN 1593 WRITE(LUPRI,*) 'NORM^2(L3AM):', 1594 * DDOT(NT1AMX*NT1AMX*NT1AMX,WORK(KL3AM),1,WORK(KL3AM),1) 1595 END IF 1596 1597*---------------------------------------------------------------------* 1598* Compute contribution from <L_3|[[H^BC, T^D_2],\tau_nu_1]|HF> 1599* and <L_3|[[H^BD, T^C_2],\tau_nu_1]|HF> 1600* and <L_3|[[H^CD, T^B_2],\tau_nu_1]|HF> 1601*---------------------------------------------------------------------* 1602 KT2AM = KEND1 1603 KEND1 = KT2AM + NT1AMX*NT1AMX 1604 LWRK1 = LWORK - KEND1 1605 IF (LWRK1 .LT. NT2AMX) THEN 1606 CALL QUIT('Insufficient space in CCSDT_HMAT_NODDY') 1607 ENDIF 1608 1609 CALL DZERO(WORK(KOME1),NT1AMX) 1610 1611 ! read T^D doubles amplitudes from file and square up 1612 ISYMD = ILSTSYM(LISTD,IDLSTD) 1613 IOPT = 2 1614 Call CC_RDRSP(LISTD,IDLSTD,ISYMD,IOPT,MODEL, 1615 & WORK(KDUM),WORK(KEND1)) 1616 Call CCLR_DIASCL(WORK(KEND1),TWO,ISYMD) 1617 CALL CC_T2SQ(WORK(KEND1),WORK(KT2AM),ISYMD) 1618 1619 CALL CC3_L3_OMEGA1_NODDY(WORK(KOME1),WORK(KL3AM), 1620 & WORK(KBCIOOOO),WORK(KBCIOVVO), 1621 & WORK(KBCIOOVV),WORK(KBCIVVVV), 1622 & WORK(KT2AM)) 1623 1624 1625 ! read T^C doubles amplitudes from file and square up 1626 ISYMC = ILSTSYM(LISTC,IDLSTC) 1627 IOPT = 2 1628 Call CC_RDRSP(LISTC,IDLSTC,ISYMC,IOPT,MODEL, 1629 & WORK(KDUM),WORK(KEND1)) 1630 Call CCLR_DIASCL(WORK(KEND1),TWO,ISYMC) 1631 CALL CC_T2SQ(WORK(KEND1),WORK(KT2AM),ISYMC) 1632 1633 CALL CC3_L3_OMEGA1_NODDY(WORK(KOME1),WORK(KL3AM), 1634 & WORK(KBDIOOOO),WORK(KBDIOVVO), 1635 & WORK(KBDIOOVV),WORK(KBDIVVVV), 1636 & WORK(KT2AM)) 1637 1638 1639 ! read T^B doubles amplitudes from file and square up 1640 ISYMB = ILSTSYM(LISTB,IDLSTB) 1641 IOPT = 2 1642 Call CC_RDRSP(LISTB,IDLSTB,ISYMB,IOPT,MODEL, 1643 & WORK(KDUM),WORK(KEND1)) 1644 Call CCLR_DIASCL(WORK(KEND1),TWO,ISYMB) 1645 CALL CC_T2SQ(WORK(KEND1),WORK(KT2AM),ISYMB) 1646 1647 CALL CC3_L3_OMEGA1_NODDY(WORK(KOME1),WORK(KL3AM), 1648 & WORK(KCDIOOOO),WORK(KCDIOVVO), 1649 & WORK(KCDIOOVV),WORK(KCDIVVVV), 1650 & WORK(KT2AM)) 1651 1652 DO I = 1,NT1AMX 1653 OMEGA1(I) = OMEGA1(I) + WORK(KOME1+I-1) 1654 END DO 1655 1656*---------------------------------------------------------------------* 1657* Compute contribution from <L_3|[H^BCD,\tau_nu_2]|HF> 1658*---------------------------------------------------------------------* 1659 CALL DZERO(WORK(KOME2),NT1AMX*NT1AMX) 1660 1661 CALL CC3_L3_OMEGA2_NODDY(WORK(KOME2),WORK(KL3AM), 1662 * WORK(KINT1SBCD),WORK(KINT2SBCD)) 1663 1664 DO I = 1,NT1AMX 1665 DO J = 1,I 1666 IJ = NT1AMX*(I-1) + J 1667 NIJ = INDEX(I,J) 1668 OMEGA2(NIJ) = OMEGA2(NIJ) + WORK(KOME2+IJ-1) 1669 END DO 1670 END DO 1671 1672 CALL QEXIT('CCSDT_HMAT_NODDY') 1673 1674 RETURN 1675 END 1676*---------------------------------------------------------------------* 1677* END OF SUBROUTINE CCSDT_HMAT_NODDY * 1678*---------------------------------------------------------------------* 1679