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_CMAT */ 21*=====================================================================* 22 SUBROUTINE CC_CMAT( ICTRAN, NCTRAN, LISTA, LISTB, LISTC, 23 & IOPTRES, FILCMA, ICDOTS,CCONS, MXVEC, 24 & WORK, LWORK ) 25*---------------------------------------------------------------------* 26* 27* Purpose: AO-direct calculation of a linear transformation of three 28* CC amplitude vectors, T^A, T^B and T^C, with the coupled 29* cluster C matrix (second derivative of the CC jacobian 30* with respect to the amplitudes) 31* 32* The linear transformations are calculated for a list 33* of T^A, T^B and T^C vectors: 34* 35* LISTA -- type of T^A vectors 36* LISTB -- type of T^B vectors 37* LISTC -- type of T^C vectors 38* ICTRAN(1,*) -- indeces of T^A vectors 39* ICTRAN(2,*) -- indeces of T^B vectors 40* ICTRAN(3,*) -- indeces of T^C vectors 41* ICTRAN(4,*) -- indeces or addresses of result vectors 42* NCTRAN -- number of requested transformations 43* FILCMA -- file name / list type of result vectors 44* or list type of vectors to be dotted on 45* ICDOTS -- indeces of vectors to be dotted on 46* CCONS -- contains the dot products on return 47* 48* return of the result vectors: 49* 50* IOPTRES = 0 : all result vectors are written to a direct 51* access file, FILCMA is used as file name 52* the start addresses of the vectors are 53* returned in ICTRAN(4,*) 54* 55* IOPTRES = 1 : the vectors are kept and returned in WORK 56* if possible, start addresses returned in 57* ICTRAN(4,*). N.B.: if WORK is not large 58* enough iopt is automatically reset to 0!!! 59* 60* IOPTRES = 3 : each result vector is written to its own 61* file by a call to CC_WRRSP, FILCMA is used 62* as list type and ICTRAN(4,*) as list index 63* NOTE that ICTRAN(4,*) is in this case input! 64* 65* IOPTRES = 4 : each result vector is added to a vector on 66* file by a call to CC_WRRSP, FILCMA is used 67* as list type and ICTRAN(4,*) as list index 68* NOTE that ICTRAN(4,*) is in this case input! 69* 70* IOPTRES = 5 : the result vectors are dotted on a array 71* of vectors, the type of the arrays given 72* by FILCMA and the indeces from ICDOTS 73* the result of the dot products is returned 74* in the CCONS array 75* 76* 77* Written by Christof Haettig, april-june 1997. 78* 79* included CC-R12: Christian Neiss, nov. 2005 80* 81*=====================================================================* 82#if defined (IMPLICIT_NONE) 83 IMPLICIT NONE 84#else 85# include "implicit.h" 86#endif 87#include "priunit.h" 88#include "ccsdinp.h" 89#include "ccsdsym.h" 90#include "maxorb.h" 91#include "mxcent.h" 92#include "ccorb.h" 93#include "cciccset.h" 94#include "cbieri.h" 95#include "distcl.h" 96#include "iratdef.h" 97#include "eritap.h" 98#include "ccisao.h" 99#include "ccfield.h" 100#include "aovec.h" 101#include "blocks.h" 102#include "r12int.h" 103#include "ccr12int.h" 104 105* local parameters: 106 CHARACTER MSGDBG*(17) 107 PARAMETER (MSGDBG='[debug] CC_CMAT> ') 108 109 LOGICAL LOCDBG 110 PARAMETER (LOCDBG = .FALSE.) 111 112 INTEGER KDUM 113 PARAMETER( KDUM = +99 999 999 ) ! dummy address for work space 114 INTEGER ISYM0 115 PARAMETER( ISYM0 = 1 ) ! symmetry of the reference state 116 INTEGER ISYOVOV 117 PARAMETER( ISYOVOV = 1 ) ! symmetry of (ia|jb) integrals 118 119 INTEGER LUCMAT 120 121 CHARACTER*(*) LISTA, LISTB, LISTC, FILCMA 122 INTEGER IOPTRES 123 INTEGER NCTRAN, MXVEC, LWORK 124 INTEGER ICTRAN(4,NCTRAN) 125 INTEGER ICDOTS(MXVEC,NCTRAN) 126 127#if defined (SYS_CRAY) 128 REAL WORK(LWORK) 129 REAL ZERO, ONE, TWO, HALF 130 REAL CCONS(MXVEC,NCTRAN) 131#else 132 DOUBLE PRECISION WORK(LWORK) 133 DOUBLE PRECISION ZERO, ONE, TWO, HALF 134 DOUBLE PRECISION CCONS(MXVEC,NCTRAN) 135#endif 136 PARAMETER (ZERO = 0.0d0, ONE = 1.0d0, TWO = 2.0d0, HALF = 0.5d0) 137 138 139 CHARACTER*(10) MODEL, MODELW 140 CHARACTER*(3) LISTXA, LISTXB, LISTXC 141 CHARACTER*(8) CBAFIL, DBAFIL 142 INTEGER INDEXA(MXCORB_CC) 143 INTEGER IOPTW, IOPT, ITRAN, IADRTH 144 INTEGER LENALL, LEN, IOFFCD, ICYCLE, ILLL, IERROR 145 INTEGER KEND0, LEND0, KENDE2, LENDE2, KENDE3, LENDE3, KEND, LEND 146 INTEGER KLIAJB, KT2AMPA, KT2AMPC, KT2AMPB, KTHETA1, KTHETA2 147 INTEGER KT1AMPA, KT1AMPB, KT1AMPC, KFOCKA, KFOCKB, KFOCKC 148 INTEGER KFABVV, KFABOO, KFBCVV, KFBCOO, KFCAVV, KFCAOO, KTHETA0 149 INTEGER KX4O, KXIAJB, KKINTC, KCBAR, KDBAR, LUCBAR, LUDBAR 150 INTEGER IDLSTA, IDLSTB, IDLSTC, ISYRES, ISYMAB, ISYCBAR, ISYDBAR 151 INTEGER ISYMTA, ISYMTB, ISYMTC, ISYMFA, ISYMFB, ISYMFC, ISYX4O 152 INTEGER ISYFAB, ISYFBC, ISYFCA, ISYMD1, NTOT, ISYMKC 153 INTEGER NTOSYM, KCCFB1, KINDXB, KFREE, LFREE, KENDSV, LWRKSV 154 INTEGER KODCL1, KODBC1, KRDBC1, KODPP1, KRDPP1, KRECNR 155 INTEGER KODCL2, KODBC2, KRDBC2, KODPP2, KRDPP2, NUMDIS 156 INTEGER IDEL2, ISYDEL, IDEL, KXINT 157 INTEGER KT1AMP0, KLAMDH0, KLAMDP0, KLAMDHA, KLAMDPA 158 INTEGER KLAMDHB, KLAMDPB, KLAMDPC, KLAMDHC 159 INTEGER KENDF1, LENDF1, KENDF2, LENDF2, KEND1, LEND1 160 INTEGER KEND2, LEND2, KEND3, LEND3, KENDA3, LENDA3 161 INTEGER KENDB3, LENDB3, IOPTTCME 162 INTEGER IDUM, IAN, KTR12AM, KVINT, LUNIT 163 164#if defined (SYS_CRAY) 165 REAL XNORM 166#else 167 DOUBLE PRECISION XNORM 168#endif 169 170 171 172 173* external functions: 174 INTEGER ILSTSYM 175 176#if defined (SYS_CRAY) 177 REAL DDOT 178#else 179 DOUBLE PRECISION DDOT 180#endif 181 182 CALL QENTER('CC_CMAT') 183 184*---------------------------------------------------------------------* 185* begin: 186*---------------------------------------------------------------------* 187 IF (LOCDBG) THEN 188 Call AROUND('ENTERED CC_CMAT') 189 IF (DIRECT) WRITE(LUPRI,'(/1X,A)') 'AO direct transformation' 190 WRITE (LUPRI,*) MSGDBG,'LISTA : ',LISTA(1:2) 191 WRITE (LUPRI,*) MSGDBG,'LISTB : ',LISTB(1:2) 192 WRITE (LUPRI,*) MSGDBG,'LISTC : ',LISTC(1:2) 193 WRITE (LUPRI,*) MSGDBG,'FILCMA: ',FILCMA 194 WRITE (LUPRI,*) MSGDBG,'NCTRAN: ',NCTRAN 195 WRITE (LUPRI,*) MSGDBG,'IOPTRES:',IOPTRES 196 WRITE (LUPRI,*) MSGDBG,'MXVEC:',MXVEC 197 WRITE (LUPRI,*) MSGDBG,'LWORK:',LWORK 198 CALL FLSHFO(LUPRI) 199 END IF 200 201 IF (CCSDT) THEN 202 WRITE(LUPRI,'(/1x,a)') 'C matrix transformations not ' 203 & //'implemented for triples yet...' 204 CALL QUIT('Triples not implemented for C '// 205 & 'matrix transformations') 206 END IF 207 208 IF ( .not. (CCS .or. CC2 .or. CCSD) ) THEN 209 WRITE(LUPRI,'(/1x,a)') 'CC_CMAT called for a Coupled Cluster ' 210 & //'method not implemented in CC_CMAT...' 211 CALL QUIT('Unknown CC method in CC_CMAT.') 212 END IF 213 214 IF ( LISTA(1:1).NE.'R' 215 & .OR. LISTB(1:1).NE.'R' 216 & .OR. LISTC(1:1).NE.'R' ) THEN 217 WRITE(LUPRI,*) 'LISTA, LISTB and LISTC must refer to', 218 & ' t-amplituded vectors in CC_CMAT.' 219 CALL QUIT('Illegal LISTA or LISTB or LISTC in CC_CMAT.') 220 END IF 221 222 IF (ISYMOP .NE. 1) THEN 223 WRITE(LUPRI,*) 'ISYMOP = ',ISYMOP 224 WRITE(LUPRI,*) 'CC_CMAT is not implemented for ISYMOP.NE.1' 225 CALL QUIT('CC_CMAT is not implemented for ISYMOP.NE.1') 226 END IF 227 228C IF (NCTRAN .GT. MAXSIM) THEN 229C WRITE(LUPRI,*) 'NCTRAN = ', NCTRAN 230C WRITE(LUPRI,*) 'MAXSIM = ', MAXSIM 231C WRITE(LUPRI,*) 'number of requested transformation is larger' 232C WRITE(LUPRI,*) 'than the maximum number of allowed ', 233C & 'simultaneous transformation.' 234C CALL QUIT('Error in CC_CMAT: NCTRAN is larger than MAXSIM.') 235C END IF 236 237* check return option for the result vectors: 238 IF (IOPTRES .EQ. 0 .OR. IOPTRES .EQ. 1) THEN 239 240 LUCMAT = -1 241 CALL WOPEN2(LUCMAT,FILCMA,64,0) 242 243 ELSE IF (IOPTRES.EQ.3 .OR. IOPTRES.EQ.4) THEN 244 CONTINUE 245 ELSE IF (IOPTRES.EQ.5) THEN 246 IF (MXVEC*NCTRAN.NE.0) CALL DZERO(CCONS,MXVEC*NCTRAN) 247 ELSE 248 CALL QUIT('Illegal value of IOPTRES in CC_CMAT.') 249 END IF 250 251* construct 'MODEL' string for CC_WRRSP routine and set write option: 252 IF (CCS) THEN 253 MODELW = 'CCS ' 254 IOPTW = 1 255 ELSE IF (CC2) THEN 256 MODELW = 'CC2 ' 257 IOPTW = 3 258 ELSE IF (CCSD) THEN 259 MODELW = 'CCSD ' 260 IOPTW = 3 261 ELSE 262 CALL QUIT('Unknown coupled cluster model in CC_CMAT.') 263 END IF 264 265 266 267* start of scratch space for the following: 268 KEND0 = 1 269 LEND0 = LWORK 270 271*=====================================================================* 272* Calculate J term and N^5 terms E1 and E2 in a loop over all required 273* C matrix transformations 274* 275* the N^5 terms H, G, I vanish for the C matrix 276* 277* All N^5 terms drop out for CCS, the E1, E2 vanish also for CC2. 278*=====================================================================* 279 280* memory allocation: 281 KLIAJB = KEND0 282 KENDE2 = KLIAJB + NT2AM(ISYOVOV) 283 LENDE2 = LWORK - KENDE2 284 285 IF (LENDE2.LT.0) THEN 286 CALL QUIT('Insufficient memory in CC_CMAT.(1)') 287 END IF 288 289* read packed (ia|jb) integrals and calculate L(ia|jb) in place: 290 Call CCG_RDIAJB(WORK(KLIAJB),NT2AM(ISYOVOV)) 291 292 IOPTTCME = 1 293 Call CCSD_TCMEPK(WORK(KLIAJB),ONE,ISYOVOV,IOPTTCME) 294 295*---------------------------------------------------------------------* 296* start loop over all C matrix transformations 297*---------------------------------------------------------------------* 298 IADRTH = 1 299 300 DO ITRAN = 1, NCTRAN 301 IDLSTA = ICTRAN(1,ITRAN) 302 IDLSTB = ICTRAN(2,ITRAN) 303 IDLSTC = ICTRAN(3,ITRAN) 304 305 IF (LOCDBG) THEN 306 WRITE (LUPRI,*) MSGDBG, 'ITRAN :',ITRAN 307 WRITE (LUPRI,*) MSGDBG, 'LISTA/IDLSTA:',LISTA(1:2),IDLSTA 308 WRITE (LUPRI,*) MSGDBG, 'LISTB/IDLSTB:',LISTB(1:2),IDLSTB 309 WRITE (LUPRI,*) MSGDBG, 'LISTC/IDLSTC:',LISTC(1:2),IDLSTC 310 CALL FLSHFO(LUPRI) 311 END IF 312*---------------------------------------------------------------------* 313* set symmetries and allocate memory: 314*---------------------------------------------------------------------* 315 ISYMTA = ILSTSYM(LISTA,IDLSTA) 316 ISYMTB = ILSTSYM(LISTB,IDLSTB) 317 ISYMTC = ILSTSYM(LISTC,IDLSTC) 318 319 ISYMFA = MULD2H(ISYOVOV,ISYMTA) 320 ISYMFB = MULD2H(ISYOVOV,ISYMTB) 321 ISYMFC = MULD2H(ISYOVOV,ISYMTC) 322 323 ISYFAB = MULD2H(ISYMFA,ISYMTB) 324 ISYFBC = MULD2H(ISYMFB,ISYMTC) 325 ISYFCA = MULD2H(ISYMFC,ISYMTA) 326 327 ISYRES = MULD2H(ISYFAB,ISYMTC) 328 329* allocate memory for double excitation response vectors: 330* (overlaps with memory for two-electron integrals) 331 KENDE3 = KLIAJB + NT2AM(ISYOVOV) 332 IF (.NOT.CCS) THEN 333 KT2AMPA = KLIAJB 334 KT2AMPB = KLIAJB 335 KT2AMPC = KLIAJB 336 337 KENDE3 = MAX(KENDE3,KLIAJB + NT2SQ(ISYMTA)) 338 KENDE3 = MAX(KENDE3,KLIAJB + NT2SQ(ISYMTB)) 339 KENDE3 = MAX(KENDE3,KLIAJB + NT2SQ(ISYMTC)) 340 END IF 341 342* allocate memory for the result vector: 343 KTHETA1 = KENDE3 344 KENDE3 = KTHETA1 + NT1AM(ISYRES) 345 346 IF (.NOT.CCS) THEN 347 KTHETA2 = KENDE3 348 KENDE3 = KTHETA2 + NT2AM(ISYRES) 349 END IF 350 351 KT1AMPA = KENDE3 352 KT1AMPB = KT1AMPA + NT1AM(ISYMTA) 353 KT1AMPC = KT1AMPB + NT1AM(ISYMTB) 354 KFOCKA = KT1AMPC + NT1AM(ISYMTC) 355 KFOCKB = KFOCKA + NT1AM(ISYMFA) 356 KFOCKC = KFOCKB + NT1AM(ISYMFB) 357 KFABVV = KFOCKC + NT1AM(ISYMFC) 358 KFABOO = KFABVV + NMATAB(ISYFAB) 359 KFBCVV = KFABOO + NMATIJ(ISYFAB) 360 KFBCOO = KFBCVV + NMATAB(ISYFBC) 361 KFCAVV = KFBCOO + NMATIJ(ISYFBC) 362 KFCAOO = KFCAVV + NMATAB(ISYFCA) 363 KENDE3 = KFCAOO + NMATIJ(ISYFCA) 364 LENDE3 = LWORK - KENDE3 365 366 IF (LENDE3 .LT. 0) THEN 367 CALL QUIT('Insufficient work space in CC_CMAT.(2)') 368 END IF 369 370*---------------------------------------------------------------------* 371* read single excitation part of the response vectors: 372*---------------------------------------------------------------------* 373* read A response amplitudes: 374 IOPT = 1 375 Call CC_RDRSP(LISTA,IDLSTA,ISYMTA,IOPT,MODEL, 376 & WORK(KT1AMPA),WORK(KDUM)) 377 378* read B response amplitudes: 379 IOPT = 1 380 Call CC_RDRSP(LISTB,IDLSTB,ISYMTB,IOPT,MODEL, 381 & WORK(KT1AMPB),WORK(KDUM)) 382 383* read C response amplitudes: 384 IOPT = 1 385 Call CC_RDRSP(LISTC,IDLSTC,ISYMTC,IOPT,MODEL, 386 & WORK(KT1AMPC),WORK(KDUM)) 387 388*---------------------------------------------------------------------* 389* calculate occupied/virtual part of first order Fock matrices: 390*---------------------------------------------------------------------* 391* calculate tF^A_ia: 392 Call CCG_LXD(WORK(KFOCKA), ISYMFA, WORK(KT1AMPA), ISYMTA, 393 & WORK(KLIAJB), ISYOVOV, 0) 394 395* calculate tF^B_ia: 396 Call CCG_LXD(WORK(KFOCKB), ISYMFB, WORK(KT1AMPB), ISYMTB, 397 & WORK(KLIAJB), ISYOVOV, 0) 398 399* calculate tF^C_ia: 400 Call CCG_LXD(WORK(KFOCKC), ISYMFC, WORK(KT1AMPC), ISYMTC, 401 & WORK(KLIAJB), ISYOVOV, 0) 402 403*---------------------------------------------------------------------* 404* calculate double one-index transformed Fock matrices: 405*---------------------------------------------------------------------* 406* calculate occ/occ block for FAB: 407 Call CCG_1ITROO(WORK(KFABOO),ISYFAB, 408 & WORK(KFOCKB),ISYMFB, 409 & WORK(KT1AMPA),ISYMTA ) 410 411 IF (LENDE3.LT.NMATIJ(ISYFAB)) 412 & CALL QUIT('Out of memory in CC_CMAT.') 413 414 Call CCG_1ITROO(WORK(KENDE3),ISYFAB, 415 & WORK(KFOCKA),ISYMFA, 416 & WORK(KT1AMPB),ISYMTB ) 417 418 Call DAXPY(NMATIJ(ISYFAB),ONE,WORK(KENDE3),1,WORK(KFABOO),1) 419 420* calculate vir/vir block for FAB: 421 Call CCG_1ITRVV(WORK(KFABVV),ISYFAB, 422 & WORK(KFOCKB),ISYMFB, 423 & WORK(KT1AMPA),ISYMTA ) 424 425 IF (LENDE3.LT.NMATAB(ISYFAB)) 426 & CALL QUIT('Out of memory in CC_CMAT.') 427 428 Call CCG_1ITRVV(WORK(KENDE3),ISYFAB, 429 & WORK(KFOCKA),ISYMFA, 430 & WORK(KT1AMPB),ISYMTB ) 431 432 Call DAXPY(NMATAB(ISYFAB),ONE,WORK(KENDE3),1,WORK(KFABVV),1) 433 434* calculate occ/occ block for FBC: 435 Call CCG_1ITROO(WORK(KFBCOO),ISYFBC, 436 & WORK(KFOCKB),ISYMFB, 437 & WORK(KT1AMPC),ISYMTC ) 438 439 IF (LENDE3.LT.NMATIJ(ISYFBC)) 440 & CALL QUIT('Out of memory in CC_CMAT.') 441 442 Call CCG_1ITROO(WORK(KENDE3),ISYFBC, 443 & WORK(KFOCKC),ISYMFC, 444 & WORK(KT1AMPB),ISYMTB ) 445 446 Call DAXPY(NMATIJ(ISYFBC),ONE,WORK(KENDE3),1,WORK(KFBCOO),1) 447 448* calculate vir/vir block for FBC: 449 Call CCG_1ITRVV(WORK(KFBCVV),ISYFBC, 450 & WORK(KFOCKB),ISYMFB, 451 & WORK(KT1AMPC),ISYMTC ) 452 453 IF (LENDE3.LT.NMATAB(ISYFBC)) 454 & CALL QUIT('Out of memory in CC_CMAT.') 455 456 Call CCG_1ITRVV(WORK(KENDE3),ISYFBC, 457 & WORK(KFOCKC),ISYMFC, 458 & WORK(KT1AMPB),ISYMTB ) 459 460 Call DAXPY(NMATAB(ISYFBC),ONE,WORK(KENDE3),1,WORK(KFBCVV),1) 461 462* calculate occ/occ block for FCA: 463 Call CCG_1ITROO(WORK(KFCAOO),ISYFCA, 464 & WORK(KFOCKC),ISYMFC, 465 & WORK(KT1AMPA),ISYMTA ) 466 467 IF (LENDE3.LT.NMATIJ(ISYFCA)) 468 & CALL QUIT('Out of memory in CC_CMAT.') 469 470 Call CCG_1ITROO(WORK(KENDE3),ISYFCA, 471 & WORK(KFOCKA),ISYMFA, 472 & WORK(KT1AMPC),ISYMTC ) 473 474 Call DAXPY(NMATIJ(ISYFCA),ONE,WORK(KENDE3),1,WORK(KFCAOO),1) 475 476* calculate vir/vir block for FCA: 477 Call CCG_1ITRVV(WORK(KFCAVV),ISYFCA, 478 & WORK(KFOCKC),ISYMFC, 479 & WORK(KT1AMPA),ISYMTA ) 480 481 IF (LENDE3.LT.NMATAB(ISYFCA)) 482 & CALL QUIT('Out of memory in CC_CMAT.') 483 484 Call CCG_1ITRVV(WORK(KENDE3),ISYFCA, 485 & WORK(KFOCKA),ISYMFA, 486 & WORK(KT1AMPC),ISYMTC ) 487 488 Call DAXPY(NMATAB(ISYFCA),ONE,WORK(KENDE3),1,WORK(KFCAVV),1) 489 490*---------------------------------------------------------------------* 491* J term 492*---------------------------------------------------------------------* 493 IF (LENDE3.LT.NT1AM(ISYRES)) 494 & CALL QUIT('Out of memory in CC_CMAT.') 495 496* initialize single-excitation part of the result vector: 497 CALL DZERO(WORK(KTHETA1),NT1AM(ISYRES)) 498 499* calculate triple one-index transformed contribution T^C x F^AB: 500 IOPT = 1 501 Call CCG_1ITRVO( WORK(KENDE3), ISYRES, 502 & WORK(KFABOO), WORK(KFABVV), ISYFAB, 503 & WORK(KT1AMPC), ISYMTC, IOPT ) 504 505 CALL DAXPY(NT1AM(ISYRES),HALF,WORK(KENDE3),1,WORK(KTHETA1),1) 506 507 508* calculate triple one-index transformed contribution T^A x F^BC: 509 IOPT = 1 510 Call CCG_1ITRVO( WORK(KENDE3), ISYRES, 511 & WORK(KFBCOO), WORK(KFBCVV), ISYFBC, 512 & WORK(KT1AMPA), ISYMTA, IOPT ) 513 514 CALL DAXPY(NT1AM(ISYRES),HALF,WORK(KENDE3),1,WORK(KTHETA1),1) 515 516 517* calculate triple one-index transformed contribution T^B x F^CA: 518 IOPT = 1 519 Call CCG_1ITRVO( WORK(KENDE3), ISYRES, 520 & WORK(KFCAOO), WORK(KFCAVV), ISYFCA, 521 & WORK(KT1AMPB), ISYMTB, IOPT ) 522 523 CALL DAXPY(NT1AM(ISYRES),HALF,WORK(KENDE3),1,WORK(KTHETA1),1) 524 525 IF (LOCDBG) THEN 526 WRITE (LUPRI,*) MSGDBG,'THETA1 after J term:' 527 CALL CC_PRP(WORK(KTHETA1),WORK(KDUM),ISYRES,1,0) 528 CALL FLSHFO(LUPRI) 529 END IF 530 531*---------------------------------------------------------------------* 532* initialize double-excitation part of the result vector: 533*---------------------------------------------------------------------* 534 IF (.NOT.CCS) THEN 535 CALL DZERO(WORK(KTHETA2),NT2AM(ISYRES)) 536 END IF 537 538*---------------------------------------------------------------------* 539* E1 and E2 terms 540*---------------------------------------------------------------------* 541 IF (.NOT. (CCS.OR.CC2.OR.CCSTST) ) THEN 542 543* check available scratch space: 544 IF ( LENDE3 .LT. NT2AM(ISYMTA) 545 & .OR. LENDE3 .LT. NT2AM(ISYMTB) 546 & .OR. LENDE3 .LT. NT2AM(ISYMTC) ) THEN 547 CALL QUIT('Insufficient work space in CC_BMAT.(3)') 548 END IF 549 550* read double excitation part of T^C vector and unpack to full matrix: 551 IOPT = 2 552 CALL CC_RDRSP(LISTC,IDLSTC,ISYMTC,IOPT,MODEL, 553 & WORK(KDUM),WORK(KENDE3)) 554 CALL CCLR_DIASCL(WORK(KENDE3),TWO,ISYMTC) 555 CALL CC_T2SQ(WORK(KENDE3),WORK(KT2AMPC),ISYMTC) 556 557* calculate the contribution to THETA2: 558 CALL CCRHS_E(WORK(KTHETA2),WORK(KT2AMPC),WORK(KFABVV), 559 & WORK(KFABOO),WORK(KENDE3),LENDE3,ISYMTC,ISYFAB) 560 561 562* read double excitation part of T^A vector and unpack to full matrix: 563 IOPT = 2 564 CALL CC_RDRSP(LISTA,IDLSTA,ISYMTA,IOPT,MODEL, 565 & WORK(KDUM),WORK(KENDE3)) 566 CALL CCLR_DIASCL(WORK(KENDE3),TWO,ISYMTA) 567 CALL CC_T2SQ(WORK(KENDE3),WORK(KT2AMPA),ISYMTA) 568 569* calculate the contribution to THETA2: 570 CALL CCRHS_E(WORK(KTHETA2),WORK(KT2AMPA),WORK(KFBCVV), 571 & WORK(KFBCOO),WORK(KENDE3),LENDE3,ISYMTA,ISYFBC) 572 573 574* read double excitation part of T^B vector and unpack to full matrix: 575 IOPT = 2 576 CALL CC_RDRSP(LISTB,IDLSTB,ISYMTB,IOPT,MODEL, 577 & WORK(KDUM),WORK(KENDE3)) 578 CALL CCLR_DIASCL(WORK(KENDE3),TWO,ISYMTB) 579 CALL CC_T2SQ(WORK(KENDE3),WORK(KT2AMPB),ISYMTB) 580 581* calculate the contribution to THETA2: 582 CALL CCRHS_E(WORK(KTHETA2),WORK(KT2AMPB),WORK(KFCAVV), 583 & WORK(KFCAOO),WORK(KENDE3),LENDE3,ISYMTB,ISYFCA) 584 585 586* restore L(ia|jb) integrals: 587 IF (ITRAN.LT.NCTRAN) THEN 588 Call CCG_RDIAJB(WORK(KLIAJB),NT2AM(ISYOVOV)) 589 IOPTTCME = 1 590 Call CCSD_TCMEPK(WORK(KLIAJB),ONE,ISYOVOV,IOPTTCME) 591 END IF 592 593 END IF 594 595*---------------------------------------------------------------------* 596* save result vector: 597*---------------------------------------------------------------------* 598 IF (LOCDBG) THEN 599 WRITE (LUPRI,*) MSGDBG,'THETA after E term:' 600 CALL CC_PRP(WORK(KTHETA1),WORK(KTHETA2),ISYRES,1,1) 601 WRITE (LUPRI,*) MSGDBG,'IOPTRES:',IOPTRES 602 WRITE (LUPRI,*) MSGDBG,'FILCMA:',FILCMA 603 WRITE (LUPRI,*) MSGDBG,'LUCMAT:',LUCMAT 604 WRITE (LUPRI,*) MSGDBG,'ICTRAN(4,ITRAN):',ICTRAN(4,ITRAN) 605 CALL FLSHFO(LUPRI) 606 END IF 607 608 KTHETA0 = -999999 609 610 611 IF (IOPTRES .EQ. 0 .OR. IOPTRES .EQ. 1) THEN 612 613* write to a common direct acces file, 614* store start address in ICTRAN(4,ITRAN) 615 616 ICTRAN(4,ITRAN) = IADRTH 617 618 CALL PUTWA2(LUCMAT,FILCMA,WORK(KTHETA1),IADRTH,NT1AM(ISYRES)) 619 IADRTH = IADRTH + NT1AM(ISYRES) 620 621 IF (.NOT.CCS) THEN 622 CALL PUTWA2(LUCMAT,FILCMA,WORK(KTHETA2),IADRTH,NT2AM(ISYRES)) 623 IADRTH = IADRTH + NT2AM(ISYRES) 624 END IF 625 626 627 ELSE IF (IOPTRES .EQ. 3) THEN 628 629* write to a sequential file by call to CC_WRRSP, 630* use FILCMA as LIST type and ICTRAN(4,ITRAN) as index 631 632 CALL CC_WRRSP(FILCMA,ICTRAN(4,ITRAN),ISYRES,IOPTW,MODELW, 633 & WORK(KTHETA0),WORK(KTHETA1),WORK(KTHETA2), 634 & WORK(KENDE3),LENDE3) 635 636 ELSE IF (IOPTRES .EQ. 4) THEN 637 638* add to a vector on a sequential file by call to CC_WARSP, 639* use FILCMA as LIST type and ICTRAN(4,ITRAN) as index 640 641 CALL CC_WARSP(FILCMA,ICTRAN(4,ITRAN),ISYRES,IOPTW,MODELW, 642 & WORK(KTHETA0),WORK(KTHETA1),WORK(KTHETA2), 643 & WORK(KENDE3),LENDE3) 644 645 ELSE IF (IOPTRES .EQ. 5) THEN 646 647* calculate required dot products 648 649 CALL CCDOTRSP(ICDOTS,CCONS,IOPTW,FILCMA,ITRAN,NCTRAN,MXVEC, 650 & WORK(KTHETA1),WORK(KTHETA2),ISYRES, 651 & WORK(KENDE3),LENDE3) 652 653 ELSE 654 CALL QUIT('illegal value for IOPT in CC_CMAT.') 655 END IF 656 657* end of loop over C matrix transformations: 658 END DO 659 660*=====================================================================* 661* end of J, E1 and E2 terms. 662*=====================================================================* 663 664 IF (LOCDBG) THEN 665 WRITE (LUPRI,*) MSGDBG,'J and E (CC2/CCSD) term finished...' 666 CALL FLSHFO(LUPRI) 667 END IF 668 669* that's it for CCS: 670 IF (CCS.OR.CCSTST) GOTO 8888 671 672*=====================================================================* 673* F term: requires AO integrals... 674* 675* Loop over AO integral shells 676* Loop over C matrix transformations 677* Loop over AO integral distributions 678* 679* Calculation of F term contributions 680* 681* End loop 682* End loop 683* End loop 684* 685* F term drops out for CCS. 686* 687*=====================================================================* 688 IF (.NOT. (CCS.OR.CCSTST) ) THEN 689 690 691*---------------------------------------------------------------------* 692* initialize integral calculation 693*---------------------------------------------------------------------* 694 695 KEND = KEND0 696 LEND = LEND0 697 698 IF (DIRECT) THEN 699 NTOSYM = 1 700 701 IF (HERDIR) THEN 702 CALL HERDI1(WORK(KEND),LEND,IPRERI) 703 ELSE 704 KCCFB1 = KEND 705 KINDXB = KCCFB1 + MXPRIM*MXCONT 706 KEND = KINDXB + (8*MXSHEL*MXCONT + 1)/IRAT 707 LEND = LWORK - KEND 708 CALL ERIDI1(KODCL1,KODCL2,KODBC1,KODBC2,KRDBC1,KRDBC2, 709 * KODPP1,KODPP2,KRDPP1,KRDPP2, 710 * KFREE,LFREE,KEND,WORK(KCCFB1),WORK(KINDXB), 711 * WORK(KEND),LEND,IPRERI) 712 713 KEND = KFREE 714 LEND = LFREE 715 END IF 716 717 KENDSV = KEND 718 LWRKSV = LEND 719 ELSE 720 NTOSYM = NSYM 721 END IF 722 723*---------------------------------------------------------------------* 724* start loop over AO integrals shells: 725*---------------------------------------------------------------------* 726 DO ISYMD1 = 1, NTOSYM 727 728 IF (DIRECT) THEN 729 IF (HERDIR) THEN 730 NTOT = MAXSHL 731 ELSE 732 NTOT = MXCALL 733 ENDIF 734 ELSE 735 NTOT = NBAS(ISYMD1) 736 END IF 737 738 DO ILLL = 1, NTOT 739 740 IF (DIRECT) THEN 741 KEND = KENDSV 742 LEND = LWRKSV 743 744 IF (HERDIR) THEN 745 CALL HERDI2(WORK(KEND),LEND,INDEXA,ILLL,NUMDIS, 746 & IPRINT) 747 ELSE 748 CALL ERIDI2(ILLL,INDEXA,NUMDIS,0,0, 749 * WORK(KODCL1),WORK(KODCL2), 750 * WORK(KODBC1),WORK(KODBC2), 751 * WORK(KRDBC1),WORK(KRDBC2), 752 * WORK(KODPP1),WORK(KODPP2), 753 * WORK(KRDPP1),WORK(KRDPP2), 754 * WORK(KCCFB1),WORK(KINDXB), 755 * WORK(KEND), LEND,IPRERI) 756 END IF 757 758 KRECNR = KEND 759 KEND = KRECNR + (NBUFX(0) - 1)/IRAT + 1 760 LEND = LWORK - KEND 761 762 IF (LEND .LT. 0) THEN 763 CALL QUIT('Insufficient work space in CC_CMAT. (1a)') 764 END IF 765 766 ELSE 767 NUMDIS = 1 768 END IF 769 770 771*---------------------------------------------------------------------* 772* start loop over C matrix transformations: 773*---------------------------------------------------------------------* 774 DO ITRAN = 1, NCTRAN 775 776 IDLSTA = ICTRAN(1,ITRAN) 777 IDLSTB = ICTRAN(2,ITRAN) 778 IDLSTC = ICTRAN(3,ITRAN) 779 780 ISYMTA = ILSTSYM(LISTA,IDLSTA) 781 ISYMTB = ILSTSYM(LISTB,IDLSTB) 782 ISYMTC = ILSTSYM(LISTC,IDLSTC) 783 784 ISYRES = MULD2H(MULD2H(ISYMTA,ISYMTB),ISYMTC) 785 786* single excitation parts of coupled cluster vectors: 787 KT1AMP0 = KEND 788 KT1AMPA = KT1AMP0 + NT1AM(ISYM0) 789 KT1AMPB = KT1AMPA + NT1AM(ISYMTA) 790 KT1AMPC = KT1AMPB + NT1AM(ISYMTB) 791 KENDF1 = KT1AMPC + NT1AM(ISYMTC) 792 793* Lambda-hole matrices: 794 KLAMDH0 = KENDF1 795 KLAMDHA = KLAMDH0 + NGLMDT(ISYM0) 796 KLAMDHB = KLAMDHA + NGLMDT(ISYMTA) 797 KLAMDHC = KLAMDHB + NGLMDT(ISYMTB) 798 KENDF1 = KLAMDHC + NGLMDT(ISYMTC) 799 800* Lambda-particle matrices: 801 KLAMDP0 = KENDF1 802 KLAMDPA = KLAMDP0 + NGLMDT(ISYM0) 803 KLAMDPB = KLAMDPA + NGLMDT(ISYMTA) 804 KLAMDPC = KLAMDPB + NGLMDT(ISYMTB) 805 KENDF1 = KLAMDPC + NGLMDT(ISYMTC) 806 807* the result vector: 808 KTHETA1 = KENDF1 809 KTHETA2 = KTHETA1 + NT1AM(ISYRES) 810 KENDF1 = KTHETA2 + NT2AM(ISYRES) 811 LENDF1 = LWORK - KENDF1 812 813 IF (LENDF1 .LT. 0) THEN 814 CALL QUIT('Insufficient memory in CC_CMAT.(1F)') 815 END IF 816 817 818* read coupled cluster vectors: 819 IOPT = 1 820 CALL CC_RDRSP('R0',0,ISYM0,IOPT,MODEL, 821 & WORK(KT1AMP0),WORK(KDUM)) 822 823 IOPT = 1 824 CALL CC_RDRSP(LISTA,IDLSTA,ISYMTA,IOPT,MODEL, 825 & WORK(KT1AMPA),WORK(KDUM)) 826 827 IOPT = 1 828 CALL CC_RDRSP(LISTB,IDLSTB,ISYMTB,IOPT,MODEL, 829 & WORK(KT1AMPB),WORK(KDUM)) 830 831 IOPT = 1 832 CALL CC_RDRSP(LISTC,IDLSTC,ISYMTC,IOPT,MODEL, 833 & WORK(KT1AMPC),WORK(KDUM)) 834 835* calculate unperturbed Lambda matrices: 836 CALL LAMMAT(WORK(KLAMDP0),WORK(KLAMDH0),WORK(KT1AMP0), 837 & WORK(KENDF1),LENDF1) 838 839* calculate response Lambda matrices: 840 CALL CCLR_LAMTRA(WORK(KLAMDP0),WORK(KLAMDPA),WORK(KLAMDH0), 841 & WORK(KLAMDHA),WORK(KT1AMPA),ISYMTA) 842 843 CALL CCLR_LAMTRA(WORK(KLAMDP0),WORK(KLAMDPB),WORK(KLAMDH0), 844 & WORK(KLAMDHB),WORK(KT1AMPB),ISYMTB) 845 846 CALL CCLR_LAMTRA(WORK(KLAMDP0),WORK(KLAMDPC),WORK(KLAMDH0), 847 & WORK(KLAMDHC),WORK(KT1AMPC),ISYMTC) 848 849 850* read result vector: 851 IF (IOPTRES.EQ.0 .OR. IOPTRES.EQ.1) THEN 852 CALL GETWA2(LUCMAT,FILCMA,WORK(KTHETA1), 853 & ICTRAN(4,ITRAN),NT1AM(ISYRES)+NT2AM(ISYRES)) 854 ELSE IF (IOPTRES.EQ.3) THEN 855 IOPT = 3 856 CALL CC_RDRSP(FILCMA,ICTRAN(4,ITRAN),ISYRES,IOPT,MODEL, 857 & WORK(KTHETA1),WORK(KTHETA2)) 858 ELSE IF (IOPTRES.EQ.4) THEN 859 CALL DZERO( WORK(KTHETA1), NT1AM(ISYRES) ) 860 CALL DZERO( WORK(KTHETA2), NT2AM(ISYRES) ) 861 ELSE IF (IOPTRES.EQ.5) THEN 862 CALL DZERO( WORK(KTHETA1), NT1AM(ISYRES) ) 863 CALL DZERO( WORK(KTHETA2), NT2AM(ISYRES) ) 864 ELSE 865 CALL QUIT('Error in CC_CMAT.') 866 END IF 867 868*---------------------------------------------------------------------* 869* loop over number of distributions on the disk: 870*---------------------------------------------------------------------* 871 DO IDEL2 = 1, NUMDIS 872 873 IF (DIRECT) THEN 874 IDEL = INDEXA(IDEL2) 875 IF (NOAUXB) THEN 876 IDUM = 1 877 CALL IJKAUX(IDEL,IDUM,IDUM,IDUM) 878 END IF 879 ISYDEL = ISAO(IDEL) 880 ELSE 881 IDEL = IBAS(ISYMD1) + ILLL 882 ISYDEL = ISYMD1 883 END IF 884 885* read AO integral distribution and calculate integrals with 886* one index transformed to occupied MO (particle): 887 888 KXINT = KENDF1 889 KENDF2 = KXINT + NDISAO(ISYDEL) 890 LENDF2 = LWORK - KENDF2 891 892 IF (LENDF2 .LT. 0) THEN 893 CALL QUIT('Insufficient work space in CC_CMAT. (3)') 894 END IF 895 896 CALL CCRDAO(WORK(KXINT),IDEL,IDEL2,WORK(KENDF2),LENDF2, 897 & WORK(KRECNR),DIRECT) 898*.....................................................................* 899 900* set option for CC_MOFCON routine: 901 IOPT = 3 902 903*.....................................................................* 904 CALL CC_MOFCON2(WORK(KXINT),WORK(KTHETA2), 905 & WORK(KLAMDPA),WORK(KLAMDHA), 906 & WORK(KLAMDPB),WORK(KLAMDHB), 907 & WORK(KLAMDPC),WORK(KLAMDH0), 908 & ISYMTA,ISYMTB,ISYMTC,ISYM0, 909 & WORK(KENDF2),LENDF2,IDEL, 910 & ISYDEL,ISYRES,ISYM0,IOPT) 911 912 CALL CC_MOFCON2(WORK(KXINT),WORK(KTHETA2), 913 & WORK(KLAMDPA),WORK(KLAMDHA), 914 & WORK(KLAMDPB),WORK(KLAMDHB), 915 & WORK(KLAMDP0),WORK(KLAMDHC), 916 & ISYMTA,ISYMTB,ISYM0,ISYMTC, 917 & WORK(KENDF2),LENDF2,IDEL, 918 & ISYDEL,ISYRES,ISYM0,IOPT) 919 920*.....................................................................* 921 CALL CC_MOFCON2(WORK(KXINT),WORK(KTHETA2), 922 & WORK(KLAMDPB),WORK(KLAMDHB), 923 & WORK(KLAMDPC),WORK(KLAMDHC), 924 & WORK(KLAMDP0),WORK(KLAMDHA), 925 & ISYMTB,ISYMTC,ISYM0,ISYMTA, 926 & WORK(KENDF2),LENDF2,IDEL, 927 & ISYDEL,ISYRES,ISYM0,IOPT) 928 929 CALL CC_MOFCON2(WORK(KXINT),WORK(KTHETA2), 930 & WORK(KLAMDPB),WORK(KLAMDHB), 931 & WORK(KLAMDP0),WORK(KLAMDH0), 932 & WORK(KLAMDPC),WORK(KLAMDHA), 933 & ISYMTB,ISYM0,ISYMTC,ISYMTA, 934 & WORK(KENDF2),LENDF2,IDEL, 935 & ISYDEL,ISYRES,ISYM0,IOPT) 936 937*.....................................................................* 938 CALL CC_MOFCON2(WORK(KXINT),WORK(KTHETA2), 939 & WORK(KLAMDPB),WORK(KLAMDHB), 940 & WORK(KLAMDPC),WORK(KLAMDHC), 941 & WORK(KLAMDPA),WORK(KLAMDH0), 942 & ISYMTB,ISYMTC,ISYMTA,ISYM0, 943 & WORK(KENDF2),LENDF2,IDEL, 944 & ISYDEL,ISYRES,ISYM0,IOPT) 945 946 CALL CC_MOFCON2(WORK(KXINT),WORK(KTHETA2), 947 & WORK(KLAMDPB),WORK(KLAMDHB), 948 & WORK(KLAMDP0),WORK(KLAMDH0), 949 & WORK(KLAMDPA),WORK(KLAMDHC), 950 & ISYMTB,ISYM0,ISYMTA,ISYMTC, 951 & WORK(KENDF2),LENDF2,IDEL, 952 & ISYDEL,ISYRES,ISYM0,IOPT) 953 954*.....................................................................* 955 956 957 END DO ! IDEL2 958*---------------------------------------------------------------------* 959* end of the loop over integral distributions: 960*---------------------------------------------------------------------* 961 IF (LOCDBG) THEN 962 WRITE (LUPRI,*) MSGDBG,'THETA after F term:' 963 WRITE (LUPRI,*) MSGDBG, 'ITRAN :',ITRAN 964 WRITE (LUPRI,*) MSGDBG, 'LISTA/IDLSTA:',LISTA(1:2), 965 & ICTRAN(1,ITRAN) 966 WRITE (LUPRI,*) MSGDBG, 'LISTB/IDLSTB:',LISTB(1:2), 967 & ICTRAN(2,ITRAN) 968 WRITE (LUPRI,*) MSGDBG, 'LISTC/IDLSTC:',LISTC(1:2), 969 & ICTRAN(3,ITRAN) 970 CALL CC_PRP(WORK(KTHETA1),WORK(KTHETA2),ISYRES,1,1) 971 WRITE (LUPRI,*) MSGDBG,'IOPTRES:',IOPTRES 972 WRITE (LUPRI,*) MSGDBG,'FILCMA:',FILCMA 973 WRITE (LUPRI,*) MSGDBG,'LUCMAT:',LUCMAT 974 WRITE (LUPRI,*) MSGDBG,'ICTRAN(4,ITRAN):',ICTRAN(4,ITRAN) 975 CALL FLSHFO(LUPRI) 976 END IF 977 978 KTHETA0 = -999999 979 980 IF (IOPTRES.EQ.0 .OR. IOPTRES.EQ.1) THEN 981 CALL PUTWA2(LUCMAT,FILCMA,WORK(KTHETA1), 982 & ICTRAN(4,ITRAN),NT1AM(ISYRES)+NT2AM(ISYRES)) 983c ELSE IF (IOPTRES.EQ.3) THEN 984c CALL CC_WRRSP(FILCMA,ICTRAN(4,ITRAN),ISYRES,IOPTW,MODELW, 985c & WORK(KTHETA0),WORK(KTHETA1),WORK(KTHETA2), 986c & WORK(KENDF1),LENDF1) 987 ELSE IF (IOPTRES.EQ.3 .OR. IOPTRES.EQ.4) THEN 988 CALL CC_WARSP(FILCMA,ICTRAN(4,ITRAN),ISYRES,IOPTW,MODELW, 989 & WORK(KTHETA0),WORK(KTHETA1),WORK(KTHETA2), 990 & WORK(KENDF1),LENDF1) 991 ELSE IF (IOPTRES.EQ.5) THEN 992 CALL CCDOTRSP(ICDOTS,CCONS,IOPTW,FILCMA,ITRAN,NCTRAN,MXVEC, 993 & WORK(KTHETA1),WORK(KTHETA2),ISYRES, 994 & WORK(KENDF1),LENDF1) 995 ELSE 996 CALL QUIT('Error in CC_CMAT.') 997 END IF 998 999*---------------------------------------------------------------------* 1000* end of the loop over C matrix transformations: 1001*---------------------------------------------------------------------* 1002 END DO ! ITRAN 1003 END DO ! ILLL 1004 END DO ! ISYMD1 1005*=====================================================================* 1006* End of Loop over AO-integrals 1007*=====================================================================* 1008 1009 IF (LOCDBG) THEN 1010 WRITE (LUPRI,*) MSGDBG,'F term section finished...' 1011 CALL FLSHFO(LUPRI) 1012 END IF 1013 1014 END IF ! (.NOT.CCS) 1015*=====================================================================* 1016* end of F term section 1017*=====================================================================* 1018 1019* that's it for CC2 1020 IF (CC2) GOTO 8888 1021 1022 1023*=====================================================================* 1024* Calculate the N^6 terms A, B, C and D in a loop over all required 1025* C matrix transformations. 1026* 1027* All N^6 terms drop out for CCS and CC2. 1028*=====================================================================* 1029 IF (.NOT. (CCS .OR. CC2 .OR. CCSTST) ) THEN 1030 1031 1032*---------------------------------------------------------------------* 1033* loop over all C matrix transformations: 1034*---------------------------------------------------------------------* 1035 DO ITRAN = 1, NCTRAN 1036 IF (LOCDBG) THEN 1037 WRITE (LUPRI,*) MSGDBG, 'N^6 term section:' 1038 WRITE (LUPRI,*) MSGDBG, 'ITRAN :',ITRAN 1039 WRITE (LUPRI,*) MSGDBG, 'LISTA/IDLSTA:',LISTA(1:2), 1040 & ICTRAN(1,ITRAN) 1041 WRITE (LUPRI,*) MSGDBG, 'LISTB/IDLSTB:',LISTB(1:2), 1042 & ICTRAN(2,ITRAN) 1043 WRITE (LUPRI,*) MSGDBG, 'LISTC/IDLSTC:',LISTC(1:2), 1044 & ICTRAN(3,ITRAN) 1045 CALL FLSHFO(LUPRI) 1046 END IF 1047 1048 ISYMTA = ILSTSYM(LISTA,ICTRAN(1,ITRAN)) 1049 ISYMTB = ILSTSYM(LISTC,ICTRAN(2,ITRAN)) 1050 ISYMTC = ILSTSYM(LISTB,ICTRAN(3,ITRAN)) 1051 ISYRES = MULD2H(MULD2H(ISYMTA,ISYMTB),MULD2H(ISYMTC,ISYOVOV)) 1052 1053* read result vector: 1054 KTHETA1 = KEND0 1055 KTHETA2 = KTHETA1 + NT1AM(ISYRES) 1056 KEND1 = KTHETA2 + NT2AM(ISYRES) 1057 LEND1 = LWORK - KEND1 1058 1059 IF (LEND1.LT.0) THEN 1060 CALL QUIT('Insufficient memory in CC_CMAT. (1b)') 1061 END IF 1062 1063 IF (IOPTRES.EQ.0 .OR. IOPTRES.EQ.1) THEN 1064 CALL GETWA2(LUCMAT,FILCMA,WORK(KTHETA1), 1065 & ICTRAN(4,ITRAN),NT1AM(ISYRES)+NT2AM(ISYRES)) 1066 ELSE IF (IOPTRES.EQ.3) THEN 1067 IOPT = 3 1068 CALL CC_RDRSP(FILCMA,ICTRAN(4,ITRAN),ISYRES,IOPT,MODEL, 1069 & WORK(KTHETA1),WORK(KTHETA2)) 1070 ELSE IF (IOPTRES.EQ.4) THEN 1071 CALL DZERO( WORK(KTHETA1), NT1AM(ISYRES) ) 1072 CALL DZERO( WORK(KTHETA2), NT2AM(ISYRES) ) 1073 ELSE IF (IOPTRES.EQ.5) THEN 1074 CALL DZERO( WORK(KTHETA1), NT1AM(ISYRES) ) 1075 CALL DZERO( WORK(KTHETA2), NT2AM(ISYRES) ) 1076 ELSE 1077 CALL QUIT('Error in CC_CMAT.') 1078 END IF 1079 1080* loop over cyclic permutations of A, B and C: 1081 DO ICYCLE = 1, 3 1082 IF (ICYCLE.EQ.1) THEN 1083 IDLSTA = ICTRAN(1,ITRAN) 1084 IDLSTB = ICTRAN(2,ITRAN) 1085 IDLSTC = ICTRAN(3,ITRAN) 1086 1087 LISTXA = LISTA 1088 LISTXB = LISTB 1089 LISTXC = LISTC 1090 ELSE IF (ICYCLE.EQ.2) THEN 1091 IDLSTB = ICTRAN(1,ITRAN) 1092 IDLSTC = ICTRAN(2,ITRAN) 1093 IDLSTA = ICTRAN(3,ITRAN) 1094 1095 LISTXB = LISTA 1096 LISTXC = LISTB 1097 LISTXA = LISTC 1098 ELSE IF (ICYCLE.EQ.3) THEN 1099 IDLSTC = ICTRAN(1,ITRAN) 1100 IDLSTA = ICTRAN(2,ITRAN) 1101 IDLSTB = ICTRAN(3,ITRAN) 1102 1103 LISTXC = LISTA 1104 LISTXA = LISTB 1105 LISTXB = LISTC 1106 ELSE 1107 CALL QUIT('Error in CC_CMAT.') 1108 END IF 1109 1110 ISYMTA = ILSTSYM(LISTXA,IDLSTA) 1111 ISYMTB = ILSTSYM(LISTXB,IDLSTB) 1112 ISYMTC = ILSTSYM(LISTXC,IDLSTC) 1113 ISYMAB = MULD2H(ISYMTA,ISYMTB) 1114 1115 1116*---------------------------------------------------------------------* 1117* read single excitation parts of T^A and T^B and keep them in 1118* memory during the loop: 1119*---------------------------------------------------------------------* 1120 KT1AMPA = KEND1 1121 KT1AMPB = KT1AMPA + NT1AM(ISYMTA) 1122 KEND2 = KT1AMPB + NT1AM(ISYMTB) 1123 LEND2 = LWORK - KEND2 1124 1125 IF (LEND2 .LE. 0) THEN 1126 CALL QUIT('Insufficient work space in CC_CMAT. (2b)') 1127 END IF 1128 1129* read single excitation part of T^A vector: 1130 IOPT = 1 1131 CALL CC_RDRSP(LISTXA,IDLSTA,ISYMTA,IOPT,MODEL, 1132 & WORK(KT1AMPA),WORK(KDUM)) 1133 1134* read single excitation part of T^B vector: 1135 IOPT = 1 1136 CALL CC_RDRSP(LISTXB,IDLSTB,ISYMTB,IOPT,MODEL, 1137 & WORK(KT1AMPB),WORK(KDUM)) 1138 1139 1140*=====================================================================* 1141* A term: 1142*=====================================================================* 1143 1144* calculate Gamma^AB: intermediate: 1145 ISYX4O = MULD2H(ISYOVOV,ISYMAB) 1146 1147 KX4O = KEND2 1148 KXIAJB = KX4O + NGAMMA(ISYX4O) 1149 KENDA3 = KXIAJB + NT2AM(ISYOVOV) 1150 LENDA3 = LWORK - KENDA3 1151 1152 IF (LENDA3 .LE. 0) THEN 1153 CALL QUIT('Insufficient work space in CC_CMAT. (A)') 1154 END IF 1155 1156* read (ia|jb) integrals from file: 1157 Call CCG_RDIAJB (WORK(KXIAJB),NT2AM(ISYOVOV)) 1158 1159* calculate double one-index transformed (ik|jl) integrals: 1160 IOPT = 1 1161 CALL CCG_4O(WORK(KX4O),ISYX4O,WORK(KXIAJB),ISYOVOV, 1162 & WORK(KT1AMPA),ISYMTA,WORK(KT1AMPB),ISYMTB, 1163 & WORK(KENDA3),LENDA3,IOPT) 1164 1165* read double excitation part of T^C vector and unpack to full matrix: 1166 KT2AMPC = KX4O + NGAMMA(ISYX4O) 1167 KENDA3 = KT2AMPC + NT2SQ(ISYMTC) 1168 LENDA3 = LWORK - KENDA3 1169 1170 IF (LENDA3 .LE. NT2AM(ISYMTC)) THEN 1171 CALL QUIT('Insufficient work space in CC_CMAT. (A2)') 1172 END IF 1173 1174 IOPT = 2 1175 CALL CC_RDRSP(LISTXC,IDLSTC,ISYMTC,IOPT,MODEL, 1176 & WORK(KDUM),WORK(KENDA3)) 1177 1178 CALL CCLR_DIASCL(WORK(KENDA3),TWO,ISYMTC) 1179 1180 CALL CC_T2SQ(WORK(KENDA3),WORK(KT2AMPC),ISYMTC) 1181 1182 1183* contract T^C with Gamma^AB to A term contribution: 1184 IOPT = 1 1185 CALL CCRHS_A(WORK(KTHETA2),WORK(KT2AMPC),WORK(KX4O), 1186 & WORK(KENDA3),LENDA3,ISYX4O,ISYMTC,IOPT) 1187 1188 IF (LOCDBG) THEN 1189 XNORM=DDOT(NGAMMA(ISYX4O),WORK(KX4O),1,WORK(KX4O),1) 1190 WRITE (LUPRI,*) 'Norm of X4O intermediate:',XNORM 1191 WRITE (LUPRI,*) MSGDBG,'THETA after A term:' 1192 WRITE (LUPRI,*) MSGDBG, 'ITRAN/ICYCLE :',ITRAN,ICYCLE 1193 WRITE (LUPRI,*) MSGDBG, 'LISTA/IDLSTA:',LISTA(1:2), 1194 & ICTRAN(1,ITRAN) 1195 WRITE (LUPRI,*) MSGDBG, 'LISTB/IDLSTB:',LISTB(1:2), 1196 & ICTRAN(2,ITRAN) 1197 WRITE (LUPRI,*) MSGDBG, 'LISTC/IDLSTC:',LISTC(1:2), 1198 & ICTRAN(3,ITRAN) 1199 XNORM=DDOT(NT2AM(ISYRES),WORK(KTHETA2),1,WORK(KTHETA2),1) 1200 WRITE (LUPRI,*) MSGDBG, 'Norm of THETA2:',XNORM 1201 CALL CC_PRP(WORK(KTHETA1),WORK(KTHETA2),ISYRES,1,1) 1202 CALL FLSHFO(LUPRI) 1203 END IF 1204 1205 1206*=====================================================================* 1207* B term: 1208*=====================================================================* 1209 ISYMKC = MULD2H(ISYOVOV,ISYMTC) 1210 1211 KKINTC = KEND2 1212 KXIAJB = KKINTC + N3ORHF(ISYMKC) 1213 KT2AMPC = KXIAJB + NT2SQ(ISYOVOV) 1214 KENDB3 = KT2AMPC + MAX(NT2AM(ISYMTC),NT2AM(ISYOVOV)) 1215 LENDB3 = LWORK - KENDB3 1216 1217 IF (LENDB3 .LT. 0) THEN 1218 CALL QUIT('Insufficient memory in CC_CMAT. (B)') 1219 END IF 1220 1221* read (ia|jb) integrals from file and unpack them to a full matrix: 1222 Call CCG_RDIAJB (WORK(KT2AMPC),NT2AM(ISYOVOV)) 1223 1224 CALL CC_T2SQ(WORK(KT2AMPC),WORK(KXIAJB),ISYOVOV) 1225 1226 1227* read (ia|jb) integrals from file: 1228 IOPT = 2 1229 CALL CC_RDRSP(LISTXC,IDLSTC,ISYMTC,IOPT,MODEL, 1230 & WORK(KDUM),WORK(KT2AMPC)) 1231 CALL CCLR_DIASCL(WORK(KT2AMPC),TWO,ISYMTC) 1232 1233* construct PHI^C intermediate: 1234 CALL CC_MI(WORK(KKINTC),WORK(KXIAJB),ISYOVOV, 1235 & WORK(KT2AMPC),ISYMTC,WORK(KENDB3),LENDB3) 1236 1237* for CCSD(R12) add correction term: 1238 IF (CCR12.AND..NOT.(CCS.OR.CC2)) THEN 1239 KTR12AM = KXIAJB 1240 KVINT = KTR12AM + NTR12AM(ISYMTC) 1241 KENDB3 = KVINT + NTR12AM(1) 1242 LENDB3 = LWORK - KENDB3 1243 IF (LENDB3.LT.0) THEN 1244 CALL QUIT('Insufficient work space in CC_CMAT (B-R12)') 1245 END IF 1246 1247 IOPT = 32 1248 CALL CC_RDRSP(LISTC,IDLSTC,ISYMTC,IOPT,MODEL,WORK(KDUM), 1249 & WORK(KTR12AM)) 1250 1251 LUNIT = -1 1252 CALL GPOPEN(LUNIT,FCCR12V,'OLD',' ','UNFORMATTED', 1253 & KDUM,.FALSE.) 1254 9999 READ(LUNIT) IAN 1255 READ(LUNIT) (WORK(KVINT-1+I), I=1, NTR12AM(1)) 1256 IF (IAN.NE.IANR12) GOTO 9999 1257 CALL GPCLOSE(LUNIT,'KEEP') 1258 1259 IOPT = 2 1260 CALL CC_R12CV(WORK(KKINTC),WORK(KTR12AM),ISYMTC,WORK(KVINT), 1261 & 1,IOPT,WORK(KENDB3),LENDB3) 1262 1263 END IF 1264 1265* double oneindex-transform PHI^C intermediate to result vector: 1266 CALL CCC_OVOV(WORK(KTHETA2),ISYRES,WORK(KKINTC),ISYMKC, 1267 & WORK(KT1AMPA),ISYMTA,WORK(KT1AMPB),ISYMTB, 1268 & WORK(KENDB3),LENDB3) 1269 1270 1271 1272 IF (LOCDBG) THEN 1273 WRITE (LUPRI,*) MSGDBG,'THETA after B term:' 1274 WRITE (LUPRI,*) MSGDBG, 'ITRAN/ICYCLE :',ITRAN,ICYCLE 1275 WRITE (LUPRI,*) MSGDBG, 'LISTA/IDLSTA:',LISTA(1:2), 1276 & ICTRAN(1,ITRAN) 1277 WRITE (LUPRI,*) MSGDBG, 'LISTB/IDLSTB:',LISTB(1:2), 1278 & ICTRAN(2,ITRAN) 1279 WRITE (LUPRI,*) MSGDBG, 'LISTC/IDLSTC:',LISTC(1:2), 1280 & ICTRAN(3,ITRAN) 1281 XNORM=DDOT(N3ORHF(ISYMKC),WORK(KKINTC),1,WORK(KKINTC),1) 1282 WRITE (LUPRI,*) MSGDBG, 'Norm of KINTC:',XNORM 1283 XNORM=DDOT(NT2AM(ISYRES),WORK(KTHETA2),1,WORK(KTHETA2),1) 1284 WRITE (LUPRI,*) MSGDBG, 'Norm of THETA2:',XNORM 1285 CALL CC_PRP(WORK(KTHETA1),WORK(KTHETA2),ISYRES,1,1) 1286 CALL FLSHFO(LUPRI) 1287 END IF 1288 1289*=====================================================================* 1290* C term: requires 5 x 1/2 V^2 O^2 !!! 1291* D term: requires 5 x 1/2 V^2 O^2 !!! 1292*=====================================================================* 1293 ISYCBAR = MULD2H(ISYMTC,ISYOVOV) 1294 ISYDBAR = MULD2H(ISYMTC,ISYOVOV) 1295 1296 KXIAJB = KEND2 1297 KT2AMPC = KXIAJB + NT2SQ(ISYOVOV) 1298 KCBAR = KT2AMPC + MAX(NT2AM(ISYMTC),NT2AM(ISYOVOV)) 1299 KEND3 = KCBAR + NT2SQ(ISYCBAR) 1300 LEND3 = LWORK - KEND3 1301 1302 KDBAR = KCBAR 1303 1304 IF (LEND3 .LT. 0) THEN 1305 CALL QUIT('Insufficient work space in CC_BMAT. (C)') 1306 END IF 1307 1308* read (ia|jb) and square them: 1309 Call CCG_RDIAJB (WORK(KT2AMPC),NT2AM(ISYOVOV)) 1310 1311 Call CC_T2SQ (WORK(KT2AMPC), WORK(KXIAJB), ISYOVOV) 1312 1313* read double excitation part of T^C vector: 1314 IOPT = 2 1315 CALL CC_RDRSP(LISTXC,IDLSTC,ISYMTC,IOPT,MODEL, 1316 & WORK(KDUM),WORK(KT2AMPC) ) 1317 1318 CALL CCLR_DIASCL(WORK(KT2AMPC),TWO,ISYMTC) 1319 1320* calculate CBAR^C intermediate: 1321 IOPT = 1 1322 CBAFIL = '--------' 1323 LUCBAR = -1 1324 IOFFCD = -1 1325 CALL CCB_CDBAR('C',WORK(KXIAJB),ISYOVOV, WORK(KT2AMPC),ISYMTC, 1326 & WORK(KCBAR), ISYCBAR, WORK(KEND3),LEND3, 1327 & CBAFIL, LUCBAR, IOFFCD, IOPT) 1328 1329* double transform CBAR^C intermediate to C term of the result vector: 1330 CALL CCB_22CD(WORK(KTHETA2),ISYRES,WORK(KCBAR),ISYCBAR, 1331 & WORK(KT1AMPA),ISYMTA,WORK(KT1AMPB),ISYMTB, 1332 & 'C', WORK(KEND3),LEND3) 1333 1334 1335* calculate DBAR^D intermediate: 1336 IOPT = 1 1337 DBAFIL = '--------' 1338 LUDBAR = -1 1339 IOFFCD = -1 1340 CALL CCB_CDBAR('D',WORK(KXIAJB),ISYOVOV, WORK(KT2AMPC),ISYMTC, 1341 & WORK(KDBAR), ISYDBAR, WORK(KEND3),LEND3, 1342 & DBAFIL, LUDBAR, IOFFCD, IOPT) 1343 1344* double transform DBAR^C intermediate to D term of the result vector: 1345 CALL CCB_22CD(WORK(KTHETA2),ISYRES,WORK(KDBAR),ISYDBAR, 1346 & WORK(KT1AMPA),ISYMTA,WORK(KT1AMPB),ISYMTB, 1347 & 'D', WORK(KEND3),LEND3) 1348 1349 1350 IF (LOCDBG) THEN 1351 WRITE (LUPRI,*) MSGDBG,'THETA after C and D terms:' 1352 WRITE (LUPRI,*) MSGDBG, 'ITRAN/ICYCLE :',ITRAN,ICYCLE 1353 WRITE (LUPRI,*) MSGDBG, 'LISTA/IDLSTA:',LISTA(1:2), 1354 & ICTRAN(1,ITRAN) 1355 WRITE (LUPRI,*) MSGDBG, 'LISTB/IDLSTB:',LISTB(1:2), 1356 & ICTRAN(2,ITRAN) 1357 WRITE (LUPRI,*) MSGDBG, 'LISTC/IDLSTC:',LISTC(1:2), 1358 & ICTRAN(3,ITRAN) 1359 XNORM=DDOT(NT2AM(ISYRES),WORK(KTHETA2),1,WORK(KTHETA2),1) 1360 WRITE (LUPRI,*) MSGDBG, 'Norm of THETA2:',XNORM 1361 CALL CC_PRP(WORK(KTHETA1),WORK(KTHETA2),ISYRES,1,1) 1362 CALL FLSHFO(LUPRI) 1363 END IF 1364 1365*---------------------------------------------------------------------* 1366* End loop over cyclic permutations, save result vector 1367*---------------------------------------------------------------------* 1368 END DO ! ICYCLE 1369 1370 KTHETA0 = -999999 1371 1372 IF (IOPTRES.EQ.0 .OR. IOPTRES.EQ.1) THEN 1373 CALL PUTWA2(LUCMAT,FILCMA,WORK(KTHETA1), 1374 & ICTRAN(4,ITRAN),NT1AM(ISYRES)+NT2AM(ISYRES)) 1375c ELSE IF (IOPTRES.EQ.3) THEN 1376c CALL CC_WRRSP(FILCMA,ICTRAN(4,ITRAN),ISYRES,IOPTW,MODELW, 1377c & WORK(KTHETA0),WORK(KTHETA1),WORK(KTHETA2), 1378c & WORK(KEND1),LEND1) 1379 ELSE IF (IOPTRES.EQ.3 .OR. IOPTRES.EQ.4) THEN 1380 CALL CC_WARSP(FILCMA,ICTRAN(4,ITRAN),ISYRES,IOPTW,MODELW, 1381 & WORK(KTHETA0),WORK(KTHETA1),WORK(KTHETA2), 1382 & WORK(KEND1),LEND1) 1383 ELSE IF (IOPTRES.EQ.5) THEN 1384 CALL CCDOTRSP(ICDOTS,CCONS,IOPTW,FILCMA,ITRAN,NCTRAN,MXVEC, 1385 & WORK(KTHETA1),WORK(KTHETA2),ISYRES, 1386 & WORK(KEND1),LEND1) 1387 ELSE 1388 CALL QUIT('Error in CC_CMAT.') 1389 END IF 1390*---------------------------------------------------------------------* 1391* End loop over C matrix transformations 1392*---------------------------------------------------------------------* 1393 END DO 1394 1395 IF (LOCDBG) THEN 1396 WRITE (LUPRI,*) MSGDBG,'N^6 term section finished...' 1397 CALL FLSHFO(LUPRI) 1398 END IF 1399 1400 END IF ! (.NOT. (CCS .OR. CC2)) 1401*=====================================================================* 1402* End of the N^6 term section 1403*=====================================================================* 1404 1405*=====================================================================* 1406* restore result vectors and clean up and return: 1407*=====================================================================* 14088888 CONTINUE 1409 1410*---------------------------------------------------------------------* 1411* if IOPTRES=1 and enough work space available, read result 1412* vectors back into memory: 1413*---------------------------------------------------------------------* 1414 1415* check size of work space: 1416 IF (IOPTRES .EQ. 1) THEN 1417 LENALL = IADRTH-1 1418 IF (LENALL .GT. LWORK) IOPTRES = 0 1419 END IF 1420 1421* read the result vectors back into memory: 1422 IF (IOPTRES .EQ. 1) THEN 1423 1424 CALL GETWA2(LUCMAT,FILCMA,WORK(1),1,LENALL) 1425 1426 IF (LOCDBG) THEN 1427 DO ITRAN = 1, NCTRAN 1428 IF (ITRAN.LT.NCTRAN) THEN 1429 LEN = ICTRAN(4,ITRAN+1)-ICTRAN(4,ITRAN) 1430 ELSE 1431 LEN = IADRTH-ICTRAN(4,NCTRAN) 1432 END IF 1433 KTHETA1 = ICTRAN(4,ITRAN) 1434 XNORM = DDOT(LEN, WORK(KTHETA1),1, WORK(KTHETA1),1) 1435 WRITE (LUPRI,*) 'Read C matrix transformation nb. ',ITRAN 1436 WRITE (LUPRI,*) 'Adress, length, NORM:',ICTRAN(4,NCTRAN), 1437 & LEN,XNORM 1438 END DO 1439 CALL FLSHFO(LUPRI) 1440 END IF 1441 END IF 1442 1443*---------------------------------------------------------------------* 1444* close C matrix file & return 1445*---------------------------------------------------------------------* 1446* check return option for the result vectors: 1447 IF (IOPTRES .EQ. 0 .OR. IOPTRES .EQ. 1) THEN 1448 CALL WCLOSE2(LUCMAT,FILCMA,'KEEP') 1449 1450 ELSE IF (IOPTRES.EQ.3 .OR. IOPTRES.EQ.4 .OR. IOPTRES.EQ.5) THEN 1451 CONTINUE 1452 ELSE 1453 CALL QUIT('Illegal value of IOPTRES in CC_CMAT.') 1454 END IF 1455 1456 1457*=====================================================================* 1458 1459 CALL QEXIT('CC_CMAT') 1460 1461 RETURN 1462 END 1463*=====================================================================* 1464* END OF SUBROUTINE CC_CMAT 1465*=====================================================================* 1466 1467*---------------------------------------------------------------------* 1468c/* Deck CC_CTST */ 1469*=====================================================================* 1470 SUBROUTINE CC_CTST(WORK,LWORK) 1471#if defined (IMPLICIT_NONE) 1472 IMPLICIT NONE 1473#else 1474# include "implicit.h" 1475#endif 1476#include "priunit.h" 1477#include "ccsdinp.h" 1478#include "ccsdsym.h" 1479#include "ccorb.h" 1480 1481* local parameters: 1482 CHARACTER MSGDBG*(18) 1483 PARAMETER (MSGDBG='[debug] CC_CTST> ') 1484 1485 LOGICAL LOCDBG 1486 PARAMETER (LOCDBG = .FALSE.) 1487 INTEGER MXCTRAN 1488 PARAMETER (MXCTRAN = 2) 1489 1490 INTEGER LWORK 1491#if defined (SYS_CRAY) 1492 REAL WORK(LWORK) 1493 REAL DDOT, RDUM 1494#else 1495 DOUBLE PRECISION WORK(LWORK) 1496 DOUBLE PRECISION DDOT, RDUM 1497#endif 1498 1499 CHARACTER*(3) LISTA, LISTB, LISTC, LISTD 1500 CHARACTER*(8) FILCMA, LABELA 1501 CHARACTER*(10) MODEL 1502 INTEGER IOPTRES, IDUM 1503 INTEGER ICTRAN(4,MXCTRAN), NCTRAN 1504 INTEGER IDLSTA, IDLSTB, IDLSTC, ISYMA, ISYMB, ISYMC, ISYMABC 1505 INTEGER KTHETA1, KTHETA2, KT1AMPC, KT2AMPC, KRESLT1, KRESLT2 1506 INTEGER KEND1, LEND1, IOPT, ISYMD, KT1AMPD, KT2AMPD, IDLSTD 1507 1508* external function: 1509 INTEGER IR1TAMP 1510 INTEGER IL1ZETA 1511 INTEGER ILSTSYM 1512 1513 1514 CALL QENTER('CC_CTST') 1515 1516 1517*---------------------------------------------------------------------* 1518* call C matrix transformation: 1519*---------------------------------------------------------------------* 1520 LISTA = 'R1' 1521 LISTB = 'R1' 1522 LISTC = 'R1' 1523 LISTD = 'L1' 1524 IDLSTA = IR1TAMP('ZDIPLEN ',.FALSE.,0.0D0,1) 1525 IDLSTB = IR1TAMP('ZDIPLEN ',.FALSE.,0.0D0,1) 1526 IDLSTC = IR1TAMP('ZDIPLEN ',.FALSE.,0.0D0,1) 1527 IDLSTD = IL1ZETA('ZDIPLEN ',.FALSE.,0.0D0,1) 1528 ICTRAN(1,1) = IDLSTA 1529 ICTRAN(2,1) = IDLSTB 1530 ICTRAN(3,1) = IDLSTC 1531 NCTRAN = 1 1532 1533 IOPTRES = 1 1534 FILCMA = 'CCCMAT' 1535 1536 WRITE (LUPRI,*) 'CCTST: call CC_CMAT:' 1537 1538 CALL CC_CMAT(ICTRAN, NCTRAN, LISTA, LISTB, LISTC, 1539 & IOPTRES, FILCMA, IDUM, RDUM, 0, WORK, LWORK ) 1540 1541 1542 ISYMA = ILSTSYM(LISTA,IDLSTA) 1543 ISYMB = ILSTSYM(LISTB,IDLSTB) 1544 ISYMC = ILSTSYM(LISTC,IDLSTC) 1545 ISYMD = ILSTSYM(LISTD,IDLSTD) 1546 ISYMABC = MULD2H(MULD2H(ISYMA,ISYMB),ISYMC) 1547 1548 KTHETA1 = ICTRAN(4,1) 1549 KTHETA2 = KTHETA1 + NT1AM(ISYMABC) 1550 KEND1 = KTHETA2 + NT2AM(ISYMABC) 1551 LEND1 = LWORK - KEND1 1552 1553 IF (NSYM.EQ.1 .AND. LOCDBG) THEN 1554 KT1AMPC = KTHETA2 + NT2AM(ISYMABC) 1555 KT2AMPC = KT1AMPC + NT1AM(ISYMC) 1556 KRESLT1 = KT2AMPC + NT2AM(ISYMC) 1557 KRESLT2 = KRESLT1 + NT1AM(ISYMABC) 1558 KEND1 = KRESLT2 + NT2AM(ISYMABC) 1559 LEND1 = LWORK - KEND1 1560 1561 IF (LEND1 .LT. 0) THEN 1562 CALL QUIT('Insufficient work space in CC_CTST.') 1563 END IF 1564 1565 IOPT = 3 1566 Call CC_RDRSP(LISTC,IDLSTC,ISYMC,IOPT,MODEL, 1567 & WORK(KT1AMPC),WORK(KT2AMPC)) 1568 1569 ! zero singles or doubles C vector: 1570C CALL DZERO(WORK(KT1AMPC),NT1AM(ISYMC)) 1571C CALL DZERO(WORK(KT2AMPC),NT2AM(ISYMC)) 1572 CALL DZERO(WORK(KRESLT1),NT1AM(ISYMABC)+NT2AM(ISYMABC)) 1573 IPRINT = 5 1574 1575 CALL CC_FDC(NT1AM(ISYMABC),NT2AM(ISYMABC), 1576 > LISTA,IDLSTA,LISTB,IDLSTB, 1577 > WORK(KT1AMPC), WORK(KRESLT1), 1578 > WORK(KEND1), LEND1) 1579 1580 IPRINT = 0 1581 1582 IF (.TRUE.) THEN 1583 WRITE (LUPRI,*) 1584 * 'LISTA, IDLSTA, ISYMA:',LISTA(1:2),IDLSTA,ISYMA 1585 WRITE (LUPRI,*) 1586 * 'LISTB, IDLSTB, ISYMB:',LISTB(1:2),IDLSTB,ISYMB 1587 WRITE (LUPRI,*) 1588 * 'LISTC, IDLSTC, ISYMC:',LISTC(1:2),IDLSTC,ISYMC 1589 WRITE (LUPRI,*) 'ISYMABC:',ISYMABC 1590 WRITE (LUPRI,*) 1591 WRITE (LUPRI,*) 'finite difference Theta vector:' 1592 Call CC_PRP(WORK(KRESLT1),WORK(KRESLT2),ISYMABC,1,1) 1593 WRITE (LUPRI,*) 'analytical Theta vector:' 1594 Call CC_PRP(WORK(KTHETA1),WORK(KTHETA2),ISYMABC,1,1) 1595 END IF 1596 1597 KT1AMPD = KEND1 1598 KT2AMPD = KT1AMPD + NT1AM(ISYMABC) 1599 KEND1 = KT2AMPD + NT2AM(ISYMABC) 1600 LEND1 = LWORK - KEND1 1601 IF (LEND1 .LT. 0) THEN 1602 CALL QUIT('Insufficient work space in CC_CTST.') 1603 END IF 1604 1605 IOPT = 3 1606 Call CC_RDRSP(LISTD,IDLSTD,ISYMD,IOPT,MODEL, 1607 & WORK(KT1AMPD),WORK(KT2AMPD)) 1608 CALL CCLR_DIASCL(WORK(KT2AMPD),0.5d0,ISYMD) 1609 1610 WRITE (LUPRI,*) 'Dot product with T^D for analytical C matrix:', 1611 > DDOT(NT1AM(ISYMD),WORK(KT1AMPD),1,WORK(KTHETA1),1)+ 1612 > DDOT(NT2AM(ISYMD),WORK(KT2AMPD),1,WORK(KTHETA2),1) 1613 1614 WRITE (LUPRI,*) 'Dot product with T^D for numerical C matrix:', 1615 > DDOT(NT1AM(ISYMD),WORK(KT1AMPD),1,WORK(KRESLT1),1)+ 1616 > DDOT(NT2AM(ISYMD),WORK(KT2AMPD),1,WORK(KRESLT2),1) 1617 1618 1619 1620 Call DAXPY(NT1AM(ISYMABC),-1.0d0,WORK(KTHETA1),1, 1621 & WORK(KRESLT1),1) 1622 IF (.NOT.CCS) THEN 1623 Call DAXPY(NT2AM(ISYMABC),-1.0d0,WORK(KTHETA2),1, 1624 & WORK(KRESLT2),1) 1625 ELSE 1626 Call DZERO(WORK(KRESLT2),NT2AM(ISYMABC)) 1627 END IF 1628 1629 WRITE (LUPRI,*) 'Norm of difference between analytical THETA ' 1630 > // 'vector and the numerical result:' 1631 WRITE (LUPRI,*) 'singles excitation part:', 1632 > DSQRT(DDOT(NT1AM(ISYMA),WORK(KRESLT1),1,WORK(KRESLT1),1)) 1633 WRITE (LUPRI,*) 'double excitation part: ', 1634 > DSQRT(DDOT(NT2AM(ISYMA),WORK(KRESLT2),1,WORK(KRESLT2),1)) 1635 1636 WRITE (LUPRI,*) 'difference vector:' 1637 Call CC_PRP(WORK(KRESLT1),WORK(KRESLT2),ISYMABC,1,1) 1638 1639 CALL FLSHFO(LUPRI) 1640 1641 1642 ELSE IF (NSYM.NE.1 .AND. LOCDBG) THEN 1643 WRITE (LUPRI,*) 'analytical Theta vector:' 1644 Call CC_PRP(WORK(KTHETA1),WORK(KTHETA2),ISYMABC,1,1) 1645 1646 WRITE (LUPRI,*) 'CC_CTST> can not calculate finite '// 1647 & 'difference C matrix' 1648 WRITE (LUPRI,*) 'CC_CTST> with symmetry.' 1649 1650 KT1AMPD = KEND1 1651 KT2AMPD = KT1AMPD + NT1AM(ISYMABC) 1652 KEND1 = KT2AMPD + NT2AM(ISYMABC) 1653 LEND1 = LWORK - KEND1 1654 IF (LEND1 .LT. 0) THEN 1655 CALL QUIT('Insufficient work space in CC_CTST.') 1656 END IF 1657 1658 IOPT = 3 1659 Call CC_RDRSP(LISTD,IDLSTD,ISYMD,IOPT,MODEL, 1660 & WORK(KT1AMPD),WORK(KT2AMPD)) 1661 CALL CCLR_DIASCL(WORK(KT2AMPD),0.5d0,ISYMD) 1662 1663 WRITE (LUPRI,*) 'Dot product with T^D for analytical C matrix:', 1664 > DDOT(NT1AM(ISYMD),WORK(KT1AMPD),1,WORK(KTHETA1),1)+ 1665 > DDOT(NT2AM(ISYMD),WORK(KT2AMPD),1,WORK(KTHETA2),1) 1666 1667 END IF 1668 1669 CALL QEXIT('CC_CTST') 1670 1671 RETURN 1672 END 1673*=====================================================================* 1674*---------------------------------------------------------------------* 1675c/* Deck CCC_OVOV */ 1676*=====================================================================* 1677 SUBROUTINE CCC_OVOV(THETA2, ISYRES, XKINT, ISYXKI, 1678 & T1AMPA, ISYMTA, T1AMPB, ISYMTB, 1679 & WORK, LWORK ) 1680*---------------------------------------------------------------------* 1681* 1682* Purpose: double one-index transform XKINT_{iljk} intermediate 1683* with T^A and T^B single-excitation amplitudes 1684* 1685* THETA2(ai,bj) += P^{AB} t^A_{ak} t^B_{bl} XKINT_{iljk} 1686* 1687* needed for CCSD part of C matrix transformation 1688* 1689* Christof Haettig, Maj 1997 1690*=====================================================================* 1691#if defined (IMPLICIT_NONE) 1692 IMPLICIT NONE 1693#else 1694# include "implicit.h" 1695#endif 1696 1697#include "priunit.h" 1698#include "ccsdsym.h" 1699#include "ccorb.h" 1700 1701#if defined (SYS_CRAY) 1702 REAL ZERO, ONE 1703#else 1704 DOUBLE PRECISION ZERO, ONE, TWO 1705#endif 1706 PARAMETER (ZERO = 0.0d0, ONE = 1.0d0, TWO = 2.0d0) 1707 1708 LOGICAL LOCDBG 1709 PARAMETER (LOCDBG = .FALSE.) 1710 1711 INTEGER ISYMTA, ISYMTB, ISYXKI, ISYRES 1712 INTEGER LWORK 1713 1714#if defined (SYS_CRAY) 1715 REAL THETA2(*) ! dimension (NT2AM(ISYRES)) 1716 REAL XKINT(*) ! dimension (N3ORHF(ISYXKI)) 1717 REAL T1AMPA(*) ! dimension (NT1AM(ISYMTA)) 1718 REAL T1AMPB(*) ! dimension (NT1AM(ISYMTB)) 1719 REAL WORK(LWORK) 1720 REAL SUM1, DDOT, SUM2 1721#else 1722 DOUBLE PRECISION THETA2(*) ! dimension (NT2AM(ISYRES)) 1723 DOUBLE PRECISION XKINT(*) ! dimension (N3ORHF(ISYXKI)) 1724 DOUBLE PRECISION T1AMPA(*) ! dimension (NT1AM(ISYMTA)) 1725 DOUBLE PRECISION T1AMPB(*) ! dimension (NT1AM(ISYMTB)) 1726 DOUBLE PRECISION WORK(LWORK) 1727 DOUBLE PRECISION SUM1, DDOT, SUM2 1728#endif 1729 1730 INTEGER ISYMA, ISYMI, ISYMB, ISYMJ, ISYMK, ISYML 1731 INTEGER ISYMAI, ISYMBJ, ISYMJL, ISYJLI, NAIBJ 1732 INTEGER KSCR1, KSCR2, KEND1, LEND1, KEND2, LEND2, KOFF1, KOFF2 1733 INTEGER KOFFX, KOFFT, NTOTJLI, NTOTA, NTOTB, NTOTJ, NAI, NBJ 1734 1735 1736* statement function: 1737 INTEGER INDEX 1738 INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 1739 1740 CALL QENTER('CCC_OVOV') 1741 1742*---------------------------------------------------------------------* 1743* begin: 1744*---------------------------------------------------------------------* 1745* check symmetries: 1746 IF (MULD2H(ISYMTA,ISYMTB) .NE. MULD2H(ISYRES,ISYXKI)) THEN 1747 CALL QUIT('Symmetry mismatch in CCC_OVOV.') 1748 END IF 1749 1750 IF (LOCDBG) THEN 1751 WRITE (LUPRI,*) 'CCC_OVOV> ISYRES:', ISYRES 1752 WRITE (LUPRI,*) 'CCC_OVOV> ISYXKI:', ISYXKI 1753 WRITE (LUPRI,*) 'CCC_OVOV> ISYMTA:', ISYMTA 1754 WRITE (LUPRI,*) 'CCC_OVOV> ISYMTB:', ISYMTB 1755 WRITE (LUPRI,*) 'CCC_OVOV> XKINT:' 1756C WRITE (LUPRI,'(5F14.6)') XKINT 1757 END IF 1758 1759 SUM1 = 0.0d0 1760 1761*---------------------------------------------------------------------* 1762* loop over A indeces: 1763*---------------------------------------------------------------------* 1764 DO ISYMA = 1, NSYM 1765 ISYMK = MULD2H(ISYMA,ISYMTA) 1766 ISYJLI = MULD2H(ISYXKI,ISYMK) 1767 1768 IF (NRHF(ISYMK).NE.0) THEN 1769 1770 1771 KSCR1 = 1 1772 KEND1 = KSCR1 + NMAIJK(ISYJLI) 1773 LEND1 = LWORK - KEND1 1774 1775 IF (LEND1.LT.0) THEN 1776 CALL QUIT('Insufficient memory in CCC_OVOV.') 1777 END IF 1778 1779 DO A = 1, NVIR(ISYMA) 1780 1781*---------------------------------------------------------------------* 1782* transform K index: sum_k XKINT_{jli;k} * t^A(ka) 1783*---------------------------------------------------------------------* 1784 KOFFX = I3ORHF(ISYJLI,ISYMK) + 1 1785 KOFFT = IT1AM(ISYMA,ISYMK) + A 1786 1787 NTOTJLI = MAX(1,NMAIJK(ISYJLI)) 1788 NTOTA = MAX(1,NVIR(ISYMA)) 1789 1790 Call DGEMV('N', NMAIJK(ISYJLI), NRHF(ISYMK), 1791 & ONE, XKINT(KOFFX), NTOTJLI, 1792 & T1AMPA(KOFFT), NTOTA, 1793 & ZERO, WORK(KSCR1), 1 ) 1794 1795*---------------------------------------------------------------------* 1796* loop over I indeces: 1797*---------------------------------------------------------------------* 1798 DO ISYMI = 1, NSYM 1799 ISYMJL = MULD2H(ISYJLI,ISYMI) 1800 ISYMAI = MULD2H(ISYMA,ISYMI) 1801 ISYMBJ = MULD2H(ISYRES,ISYMAI) 1802 1803 KSCR2 = KEND1 1804 KEND2 = KSCR2 + NT1AM(ISYMBJ) 1805 LEND2 = LWORK - KEND2 1806 1807 IF (LEND2.LT.0) THEN 1808 CALL QUIT('Insufficient memory in CCC_OVOV.') 1809 END IF 1810 1811 1812 DO I = 1, NRHF(ISYMI) 1813 1814*---------------------------------------------------------------------* 1815* transform L index: sum_l t^B(bl) * SCR_{jl;i} 1816*---------------------------------------------------------------------* 1817 DO ISYMB = 1, NSYM 1818 ISYMJ = MULD2H(ISYMBJ,ISYMB) 1819 ISYML = MULD2H(ISYMTB,ISYMB) 1820 IF (MULD2H(ISYMJ,ISYML).NE.ISYMJL) 1821 & CALL QUIT('Error in CCC_OVOV (1).') 1822 1823 KOFF1 = KSCR1 + IMAIJK(ISYMJL,ISYMI) + 1824 & NMATIJ(ISYMJL)*(I-1) + IMATIJ(ISYMJ,ISYML) 1825 KOFF2 = KSCR2 + IT1AM(ISYMB,ISYMJ) 1826 KOFFT = IT1AM(ISYMB,ISYML) + 1 1827 1828 NTOTJ = MAX(1,NRHF(ISYMJ)) 1829 NTOTB = MAX(1,NVIR(ISYMB)) 1830 1831 Call DGEMM('N','T',NVIR(ISYMB),NRHF(ISYMJ),NRHF(ISYML), 1832 & ONE, T1AMPB(KOFFT), NTOTB, WORK(KOFF1), NTOTJ, 1833 & ZERO, WORK(KOFF2), NTOTB ) 1834 END DO 1835 1836*---------------------------------------------------------------------* 1837* store result in THETA2 vector: 1838*---------------------------------------------------------------------* 1839 NAI = IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I-1) + A 1840 1841 IF (ISYMAI .EQ. ISYMBJ) THEN 1842 1843 WORK(KSCR2-1+NAI) = TWO * WORK(KSCR2-1+NAI) 1844 1845 DO NBJ = 1, NT1AM(ISYMBJ) 1846 NAIBJ = IT2AM(ISYMAI,ISYMBJ) + INDEX(NAI,NBJ) 1847 THETA2(NAIBJ) = THETA2(NAIBJ) + WORK(KSCR2-1+NBJ) 1848 END DO 1849 1850 ELSE IF (ISYMAI .LT. ISYMBJ) THEN 1851 1852 DO NBJ = 1, NT1AM(ISYMBJ) 1853 NAIBJ = IT2AM(ISYMAI,ISYMBJ) + NT1AM(ISYMAI)*(NBJ-1)+NAI 1854 THETA2(NAIBJ) = THETA2(NAIBJ) + WORK(KSCR2-1+NBJ) 1855 END DO 1856 1857 ELSE IF (ISYMBJ .LT. ISYMAI) THEN 1858 1859 DO NBJ = 1, NT1AM(ISYMBJ) 1860 NAIBJ = IT2AM(ISYMAI,ISYMBJ) + NT1AM(ISYMBJ)*(NAI-1)+NBJ 1861 THETA2(NAIBJ) = THETA2(NAIBJ) + WORK(KSCR2-1+NBJ) 1862 END DO 1863 1864 ELSE 1865 CALL QUIT('Error in CCC_OVOV (2).') 1866 END IF 1867 1868 1869 1870*---------------------------------------------------------------------* 1871* end loops over I and A: 1872*---------------------------------------------------------------------* 1873 END DO ! I 1874 END DO ! ISYMI 1875 1876 END DO ! A 1877 END IF ! NRHF(ISYMK).NE.0 1878 END DO ! ISYMA 1879 1880 CALL QEXIT('CCC_OVOV') 1881 1882 RETURN 1883 END 1884*---------------------------------------------------------------------* 1885* END OF SUBROUTINE CCC_OVOV * 1886*---------------------------------------------------------------------* 1887