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_BAMAT */ 21*=====================================================================* 22 SUBROUTINE CC_BAMAT(IBATRAN,NBATRAN,LISTO,LISTA,LISTB,IOPTRES, 23 & FILBAM, IBDOTS, BCONS,MXVEC,WORK,LWORK ) 24*---------------------------------------------------------------------* 25* 26* Purpose: linear transformation of two CC amplitude vectors, 27* T^A and T^B, with the CC B{O} matrix 28* 29* The linear transformations are calculated for a list 30* of operators and T^A and T^B vectors: 31* 32* LISTO -- 'o1' 33* LISTA -- type of T^A vectors 34* LISTB -- type of T^B vectors 35* IBATRAN(1,*) -- indeces of the operators 36* IBATRAN(2,*) -- indeces of T^A vectors 37* IBATRAN(3,*) -- indeces of T^B vectors 38* IBATRAN(4,*) -- indeces or addresses of result vectors 39* NBATRAN -- number of requested transformations 40* FILBAM -- file name / list type of result vectors 41* or list type of vectors to be dotted on 42* IBDOTS -- indeces of vectors to be dotted on 43* BCONS -- contains the dot products on return 44* 45* return of the result vectors: 46* 47* IOPTRES = 0 : all result vectors are written to a direct 48* access file, FILBAM is used as file name 49* the start addresses of the vectors are 50* returned in IBATRAN(4,*) 51* 52* IOPTRES = 1 : the vectors are kept and returned in WORK 53* if possible, start addresses returned in 54* IBATRAN(4,*). N.B.: if WORK is not large 55* enough iopt is automatically reset to 0!!! 56* 57* IOPTRES = 3 : each result vector is written to its own 58* file by a call to CC_WRRSP, FILBAM is used 59* as list type and IBATRAN(4,*) as list index 60* NOTE that IBATRAN(4,*) is in this case input! 61* 62* IOPTRES = 4 : each result vector is written/added to a 63* file by a call to CC_WARSP, FILBAM is used 64* as list type and IBATRAN(4,*) as list index 65* NOTE that IBATRAN(4,*) is in this case input! 66* 67* IOPTRES = 5 : the result vectors are dotted on a array 68* of vectors, the type of the arrays given 69* by FILBAM and the indeces from IBDOTS 70* the result of the dot products is returned 71* in the BCONS array 72* 73* Written by Christof Haettig, Maj 1997. 74* 75* Adapted for CC-R12: Chrsitian Neiss, July 2005 76* 77*=====================================================================* 78#if defined (IMPLICIT_NONE) 79 IMPLICIT NONE 80#else 81# include "implicit.h" 82#endif 83#include "ccsdinp.h" 84#include "priunit.h" 85#include "ccsdsym.h" 86#include "ccorb.h" 87#include "ccroper.h" 88#include "cclists.h" 89#include "r12int.h" 90#include "ccr12int.h" 91 92* local parameters: 93 CHARACTER MSGDBG*(18) 94 PARAMETER (MSGDBG='[debug] CC_BAMAT> ') 95 96 LOGICAL LOCDBG 97 PARAMETER (LOCDBG = .FALSE.) 98 99 INTEGER KDUM 100 PARAMETER ( KDUM = +99 999 999 ) ! dummy address for work space 101 102 INTEGER ISYMT0 103 PARAMETER ( ISYMT0 = 1 ) ! symmetry of the reference state 104 105 INTEGER LUBMAT 106 107 CHARACTER*(*) LISTO, LISTA, LISTB, FILBAM 108 INTEGER IOPTRES 109 INTEGER NBATRAN, MXVEC, LWORK 110 INTEGER IBATRAN(MXDIM_BATRAN,NBATRAN) 111 INTEGER IBDOTS(MXVEC,NBATRAN) 112 113#if defined (SYS_CRAY) 114 REAL WORK(LWORK) 115 REAL BCONS(MXVEC,NBATRAN) 116 REAL ZERO, ONE, TWO, HALF 117 REAL XNORM 118#else 119 DOUBLE PRECISION WORK(LWORK) 120 DOUBLE PRECISION BCONS(MXVEC,NBATRAN) 121 DOUBLE PRECISION ZERO, ONE, TWO, HALF 122 DOUBLE PRECISION XNORM 123#endif 124 PARAMETER (ZERO = 0.0d0, ONE = 1.0d0, TWO = 2.0d0, HALF = 0.5d0) 125 126 CHARACTER*(10) MODEL, MODELW 127 CHARACTER*(8) LABELO 128 129 INTEGER KEND1, KEND2, LEND2, LENALL, KEND3, LEND3 130 INTEGER ISYMTA, ISYMTB, ISYMO, ISYRES, ISYMAO, ISYMBO, IOPTW 131 INTEGER IDLSTO, IDLSTA, IDLSTB, IOPTWE 132 INTEGER KPERT, KOO, KOV, KVV, KAOO, KBOO, KAVV, KBVV, KSCR 133 INTEGER KT1AMPA, KT1AMPB, KT1AMP0, KLAMDP0, KLAMDH0 134 INTEGER KTHETA0, KTHETA1, KTHETA2, KT2AMPA, KT2AMPB 135 INTEGER IRREP, IERR, IERRB, IOPT, ISYM, LEN, ITRAN, IADRTH, 136 & KTHETA1EFF, KTHETA2EFF 137 INTEGER IOPTWR12,LENMOD,KTHETAR12,KATRANR12 138 INTEGER KLAMDPB,KLAMDHB,KLAMDPA,KLAMDHA,KT12AMP,KXINTTRI,KXINTSQ 139 INTEGER LUNIT,IAN,DUMMY,ISYM1,ISYM2,IDLST1,IDLST2 140 CHARACTER APROXR12*3,LIST1*3,LIST2*3 141 142* external functions: 143 INTEGER ILSTSYM 144#if defined (SYS_CRAY) 145 REAL DDOT, FREQLST, FREQB 146#else 147 DOUBLE PRECISION DDOT, FREQLST, FREQB 148#endif 149 150*---------------------------------------------------------------------* 151* begin: 152*---------------------------------------------------------------------* 153 IF (LOCDBG) THEN 154 Call AROUND('ENTERED CC_BAMAT') 155 WRITE (LUPRI,*) 'LISTO : ',LISTO 156 WRITE (LUPRI,*) 'LISTA : ',LISTA 157 WRITE (LUPRI,*) 'LISTB : ',LISTB 158 WRITE (LUPRI,*) 'FILBAM: ',FILBAM 159 WRITE (LUPRI,*) 'NBATRAN: ',NBATRAN 160 WRITE (LUPRI,*) 'IOPTRES:',IOPTRES 161 CALL FLSHFO(LUPRI) 162 END IF 163 164 ! well, this is no longer true, since CC3 is implemented 165 ! but it is not yet debugged... though, be aware!! 166 IF (CCSDT) THEN 167 WRITE(LUPRI,'(/1x,a)') 'B{O} matrix transformations not ' 168 & //'implemented for triples yet...' 169 CALL QUIT('Triples not implemented for B '// 170 & 'matrix transformations') 171 END IF 172 173 IF ( .NOT. (CCS .OR. CC2 .OR. CCSD .OR. CC3) ) THEN 174 WRITE(LUPRI,'(/1x,a)') 'CC_BAMAT called for a Coupled Cluster ' 175 & //'method not implemented in CC_BAMAT...' 176 CALL QUIT('Unknown CC method in CC_BAMAT.') 177 END IF 178 179 IF (LISTO(1:2).NE.'o1') THEN 180 WRITE (LUPRI,*) 'LISTO must refer to operator list '// 181 & 'o1 in CC_BAMAT.' 182 CALL QUIT('Illegal LISTO in CC_BAMAT.') 183 END IF 184 185 IF (LISTA(1:1).NE.'R' .OR. LISTB(1:1).NE.'R') THEN 186 WRITE(LUPRI,*) 'LISTA and LISTB must refer to t-amplitude', 187 & ' vectors in CC_BAMAT.' 188 CALL QUIT('Illegal LISTA or LISTB in CC_BAMAT.') 189 END IF 190 191 IF (CCS) THEN 192 MODELW = 'CCS ' 193 IOPTW = 1 194 ELSE IF (CC2) THEN 195 MODELW = 'CC2 ' 196 IOPTW = 3 197 ELSE IF (CCSD) THEN 198 MODELW = 'CCSD ' 199 IOPTW = 3 200 ELSE IF (CC3) THEN 201 MODELW = 'CC3 ' 202 IOPTW = 3 203 IOPTWE = 24 204 ELSE 205 CALL QUIT('Unknown coupled cluster model in CC_BAMAT.') 206 END IF 207 IF (CCR12) THEN 208 APROXR12 = ' ' 209 CALL CCSD_MODEL(MODELW,LENMOD,10,MODELW,10,APROXR12) 210 IOPTWR12 = 32 211 END IF 212 213* check return option for the result vectors: 214 IF (IOPTRES .EQ. 0 .OR. IOPTRES .EQ. 1) THEN 215 216 LUBMAT = -1 217 CALL WOPEN2(LUBMAT,FILBAM,64,0) 218 219 ELSE IF (IOPTRES.EQ.3 .OR. IOPTRES.EQ.4) THEN 220 CONTINUE 221 ELSE IF (IOPTRES.EQ.5) THEN 222 IF (MXVEC*NBATRAN.NE.0) CALL DZERO(BCONS,MXVEC*NBATRAN) 223 ELSE 224 CALL QUIT('Illegal value of IOPTRES in CC_BAMAT.') 225 END IF 226 227*=====================================================================* 228* calculate B matrix transformations: 229*=====================================================================* 230 231 KEND1 = 1 232 233 IADRTH = 1 234 235 DO ITRAN = 1, NBATRAN 236 237 IDLSTO = IBATRAN(1,ITRAN) 238 IDLSTA = IBATRAN(2,ITRAN) 239 IDLSTB = IBATRAN(3,ITRAN) 240 241 LABELO = LBLOPR(IDLSTO) 242 243 ISYMO = ILSTSYM(LISTO,IDLSTO) 244 ISYMTA = ILSTSYM(LISTA,IDLSTA) 245 ISYMTB = ILSTSYM(LISTB,IDLSTB) 246 ISYMAO = MULD2H(ISYMTA,ISYMO) 247 ISYMBO = MULD2H(ISYMTB,ISYMO) 248 ISYRES = MULD2H(MULD2H(ISYMTA,ISYMTB),MULD2H(ISYMT0,ISYMO)) 249 250*---------------------------------------------------------------------* 251* allocate work space for the result vector(s): 252*---------------------------------------------------------------------* 253 IF (CCS) THEN 254 KTHETA1 = KEND1 255 KTHETA2 = KDUM 256 KEND2 = KTHETA1 + NT1AM(ISYRES) 257 ELSE 258 KTHETA1 = KEND1 259 KTHETA2 = KTHETA1 + NT1AM(ISYRES) 260 KEND2 = KTHETA2 + NT2AM(ISYRES) 261 IF (CCR12) THEN 262 KTHETAR12 = KTHETA2 + NT2AM(ISYRES) 263 KEND2 = KTHETAR12 + NTR12AM(ISYRES) 264 END IF 265 END IF 266 IF (CCSDT .AND. IOPTRES.NE.5) THEN 267 KTHETA1EFF = KEND2 268 KTHETA2EFF = KTHETA1EFF + NT1AM(ISYRES) 269 KEND2 = KTHETA2EFF + NT2AM(ISYRES) 270 END IF 271 272 IF (LOCDBG) THEN 273 WRITE (LUPRI,*) 'B{O} matrix transformation for ITRAN,',ITRAN 274 WRITE (LUPRI,*) 'IADRTH:',IADRTH 275 WRITE (LUPRI,*) 'LISTO,IDLSTO:',LISTO,IDLSTO 276 WRITE (LUPRI,*) 'LISTA,IDLSTA:',LISTA,IDLSTA 277 WRITE (LUPRI,*) 'LISTB,IDLSTB:',LISTB,IDLSTB 278 WRITE (LUPRI,*) 'ISYMO,ISYMTA,ISYMTB:',ISYMO,ISYMTA,ISYMTB 279 CALL FLSHFO(LUPRI) 280 END IF 281 282*---------------------------------------------------------------------* 283* allocate memory for property integrals and response vectors: 284* 1) AO/MO perturbation integrals (complete matrix) 285* 2) MO transformed perturbation integrals (occ/occ block) 286* 3) MO transformed perturbation integrals (vir/vir block) 287* 4) MO transformed perturbation integrals (occ/vir block) 288* 5) singles excitation part of T^A 289* 6) singles excitation part of T^B 290* 7) singles excitation part of T^0 291* 8) zeroth-order lambda particle matrix 292* 9) zeroth-order lambda hole matrix 293* 10) a scratch vector with the size of the singles result vector 294*---------------------------------------------------------------------* 295 KPERT = KEND2 296 KOO = KPERT + N2BST(ISYMO) 297 KOV = KOO + NMATIJ(ISYMO) 298 KVV = KOV + NT1AM(ISYMO) 299 KAOO = KVV + NMATAB(ISYMO) 300 KBOO = KAOO + NMATIJ(ISYMAO) 301 KAVV = KBOO + NMATIJ(ISYMBO) 302 KBVV = KAVV + NMATAB(ISYMAO) 303 KT1AMPA = KBVV + NMATAB(ISYMBO) 304 KT1AMPB = KT1AMPA + NT1AM(ISYMTA) 305 KT1AMP0 = KT1AMPB + NT1AM(ISYMTB) 306 KLAMDP0 = KT1AMP0 + NT1AM(ISYMT0) 307 KLAMDH0 = KLAMDP0 + NGLMDT(ISYMT0) 308 KSCR = KLAMDH0 + NGLMDT(ISYMT0) 309 KEND2 = KSCR + NT1AM(ISYRES) 310 LEND2 = LWORK - KEND2 311 312 IF (LEND2 .LE. 0) THEN 313 CALL QUIT('Insufficient work space in CC_BAMAT. (8)') 314 END IF 315 316* read single excitation part of T^A: 317 IOPT = 1 318 CALL CC_RDRSP(LISTA,IDLSTA,ISYMTA,IOPT,MODEL, 319 & WORK(KT1AMPA),WORK(KDUM)) 320 321* read single excitation part of T^B: 322 IOPT = 1 323 CALL CC_RDRSP(LISTB,IDLSTB,ISYMTB,IOPT,MODEL, 324 & WORK(KT1AMPB),WORK(KDUM)) 325 326* read single excitation part of zeroth-order coupled cluster vector: 327 IOPT = 1 328 CALL CC_RDRSP('R0',0,ISYMT0,IOPT,MODEL, 329 & WORK(KT1AMP0),WORK(KDUM)) 330 331* get zeroth-order Lambda matrices: 332 CALL LAMMAT(WORK(KLAMDP0),WORK(KLAMDH0),WORK(KT1AMP0), 333 & WORK(KEND2),LEND2) 334 335* read the AO integrals: 336 CALL CCPRPAO(LABELO,.TRUE.,WORK(KPERT),IRREP,ISYM,IERR, 337 & WORK(KEND2),LEND2) 338 IF (IERR.NE.0 .OR. IRREP.NE.ISYMO) THEN 339 CALL QUIT('CC_BAMAT: error while reading operator '//LABELO) 340 END IF 341 342* transform one-electron integrals in place: 343 CALL CC_FCKMO(WORK(KPERT),WORK(KLAMDP0),WORK(KLAMDH0), 344 & WORK(KEND2),LEND2,ISYMO,1,1) 345 346* gather occ/occ, occ/vir, and vir/vir blocks: 347 CALL CC_GATHEROO(WORK(KPERT),WORK(KOO),ISYMO) 348 CALL CC_GATHEROV(WORK(KPERT),WORK(KOV),ISYMO) 349 CALL CC_GATHERVV(WORK(KPERT),WORK(KVV),ISYMO) 350 351*---------------------------------------------------------------------* 352* calculate O^{A} = [O, T1^A] 353*---------------------------------------------------------------------* 354* Ftilde^{A}, occupied/occupied blocks: 355 CALL CCG_1ITROO(WORK(KAOO),ISYMAO, 356 & WORK(KOV), ISYMO, WORK(KT1AMPA),ISYMTA ) 357 358* Ftilde^{B}, occupied/occupied blocks: 359 CALL CCG_1ITROO(WORK(KBOO),ISYMBO, 360 & WORK(KOV), ISYMO, WORK(KT1AMPB),ISYMTB ) 361 362* Ftilde^{A}, virtual/virtual blocks: 363 CALL CCG_1ITRVV(WORK(KAVV),ISYMAO, 364 & WORK(KOV), ISYMO, WORK(KT1AMPA),ISYMTA ) 365 366* Ftilde^{B}, virtual/virtual blocks: 367 CALL CCG_1ITRVV(WORK(KBVV),ISYMBO, 368 & WORK(KOV), ISYMO, WORK(KT1AMPB),ISYMTB ) 369 370 371 IF (LOCDBG) THEN 372 XNORM=DDOT(NMATIJ(ISYMTA),WORK(KAOO),1,WORK(KAOO),1) 373 WRITE (LUPRI,*) 'Norm of O^AOO:',XNORM 374 XNORM=DDOT(NMATAB(ISYMTA),WORK(KAVV),1,WORK(KAVV),1) 375 WRITE (LUPRI,*) 'Norm of O^AVV:',XNORM 376 XNORM=DDOT(NMATIJ(ISYMTB),WORK(KBOO),1,WORK(KBOO),1) 377 WRITE (LUPRI,*) 'Norm of O^BOO:',XNORM 378 XNORM=DDOT(NMATAB(ISYMTB),WORK(KBVV),1,WORK(KBVV),1) 379 WRITE (LUPRI,*) 'Norm of O^BVV:',XNORM 380 WRITE (LUPRI,*) 'T^A (singles only):' 381 Call CC_PRP(WORK(KT1AMPA),WORK,ISYMTA,1,0) 382 WRITE (LUPRI,*) 'T^B (singles only):' 383 Call CC_PRP(WORK(KT1AMPB),WORK,ISYMTB,1,0) 384 CALL FLSHFO(LUPRI) 385 END IF 386 387*---------------------------------------------------------------------* 388* initialize the singles part of the result vector THETA: 389*---------------------------------------------------------------------* 390 CALL DZERO(WORK(KTHETA1),NT1AM(ISYRES)) 391 392*---------------------------------------------------------------------* 393* J contribution: vir/occ block of double transformed integrals: 394*---------------------------------------------------------------------* 395 396* 1. contribution [O^A,T^B]: 397 IOPT = 1 398 CALL CCG_1ITRVO(WORK(KSCR),ISYRES, WORK(KBOO), 399 & WORK(KBVV),ISYMBO, 400 & WORK(KT1AMPA),ISYMTA,IOPT ) 401 402 CALL DAXPY(NT1AM(ISYRES),HALF,WORK(KSCR),1,WORK(KTHETA1),1) 403 404* 2. contribution [O^B,T^A]: 405 IOPT = 1 406 CALL CCG_1ITRVO(WORK(KSCR),ISYRES, WORK(KAOO), 407 & WORK(KAVV),ISYMAO, 408 & WORK(KT1AMPB),ISYMTB,IOPT ) 409 410 CALL DAXPY(NT1AM(ISYRES),HALF,WORK(KSCR),1,WORK(KTHETA1),1) 411 412 413 IF (LOCDBG) THEN 414 XNORM=DDOT(NT1AM(ISYRES),WORK(KTHETA1),1,WORK(KTHETA1),1) 415 WRITE (LUPRI,*) 'Norm of O^ABOV (J contribution):',XNORM 416C WRITE (LUPRI,'(/5X,A)') 'AVV / AOO:' 417C CALL CC_PREI(WORK(KAVV),WORK(KAOO),ISYMAO,1) 418C WRITE (LUPRI,'(/5X,A)') 'BVV / BOO:' 419C CALL CC_PREI(WORK(KBVV),WORK(KBOO),ISYMBO,1) 420 CALL FLSHFO(LUPRI) 421 END IF 422 423 424*---------------------------------------------------------------------* 425* initialize the doubles part of the result vector THETA: 426*---------------------------------------------------------------------* 427 IF (.NOT. CCS ) THEN 428 CALL DZERO( WORK(KTHETA2), NT2AM(ISYRES) ) 429 END IF 430 431*---------------------------------------- 432* first E1 & E2 contributions, T^B x F^A: 433*---------------------------------------- 434C IF (.NOT. (CCS .OR. CC2 .OR. CCSTST) ) THEN 435 IF (.NOT. (CCS .OR. CCSTST) ) THEN 436 KT2AMPB = KEND2 437 KEND3 = KT2AMPB + NT2SQ(ISYMTB) 438 LEND3 = LWORK - KEND3 439 440 IF (LEND3 .LT. NT2AM(ISYMTB)) THEN 441 CALL QUIT('Insufficient work space in CC_BAMAT. (15)') 442 END IF 443 444 IOPT = 2 445 CALL CC_RDRSP(LISTB,IDLSTB,ISYMTB,IOPT, 446 & MODEL,WORK(KDUM),WORK(KEND3)) 447 448 CAll CCLR_DIASCL(WORK(KEND3),TWO,ISYMTB) 449 450 CALL CC_T2SQ(WORK(KEND3),WORK(KT2AMPB),ISYMTB) 451 452* calculate the contribution to THETA2: 453 CALL CCRHS_E(WORK(KTHETA2),WORK(KT2AMPB),WORK(KAVV), 454 & WORK(KAOO),WORK(KEND3),LEND3,ISYMTB,ISYMAO) 455 456 IF (LOCDBG) THEN 457 XNORM=DDOT(NMATIJ(ISYMAO),WORK(KAOO),1,WORK(KAOO),1) 458 WRITE (LUPRI,*) 'Norm of KAOO intermediate:',XNORM 459 XNORM=DDOT(NMATAB(ISYMAO),WORK(KAVV),1,WORK(KAVV),1) 460 WRITE (LUPRI,*) 'Norm of KAVV intermediate:',XNORM 461 XNORM=DDOT(NT2AM(ISYRES),WORK(KTHETA2),1,WORK(KTHETA2),1) 462 WRITE (LUPRI,*) 'Norm of THETA2 after first E contribution:', 463 & XNORM 464 Call AROUND('final result for B{O} matrix transformation:') 465 Call CC_PRP(WORK(KTHETA1),WORK(KTHETA2),ISYRES,1,1) 466 CALL FLSHFO(LUPRI) 467 END IF 468 469 END IF ! (.NOT. (CCS .OR. CC2 .OR. CCSTST)) 470 471*----------------------------------------- 472* second E1 & E2 contributions, T^A x F^B: 473*----------------------------------------- 474 IF (.NOT. (CCS .OR. CCSTST) ) THEN 475 KT2AMPA = KEND2 476 KEND3 = KT2AMPA + NT2SQ(ISYMTA) 477 LEND3 = LWORK - KEND3 478 479 IF (LEND3 .LT. NT2AM(ISYMTA)) THEN 480 CALL QUIT('Insufficient work space in CC_BAMAT. (16)') 481 END IF 482 483 IOPT = 2 484 CALL CC_RDRSP(LISTA,IDLSTA,ISYMTA,IOPT, 485 & MODEL,WORK(KDUM),WORK(KEND3)) 486 487 CAll CCLR_DIASCL(WORK(KEND3),TWO,ISYMTA) 488 489 CALL CC_T2SQ(WORK(KEND3),WORK(KT2AMPA),ISYMTA) 490 491* calculate the contribution to THETA2: 492 CALL CCRHS_E(WORK(KTHETA2),WORK(KT2AMPA),WORK(KBVV), 493 & WORK(KBOO),WORK(KEND3),LEND3,ISYMTA,ISYMBO) 494 495 IF (LOCDBG) THEN 496 XNORM=DDOT(NMATIJ(ISYMBO),WORK(KBOO),1,WORK(KBOO),1) 497 WRITE (LUPRI,*) 'Norm of KBOO intermediate:',XNORM 498 XNORM=DDOT(NMATAB(ISYMBO),WORK(KBVV),1,WORK(KBVV),1) 499 WRITE (LUPRI,*) 'Norm of KBVV intermediate:',XNORM 500 XNORM=DDOT(NT2AM(ISYRES),WORK(KTHETA2),1,WORK(KTHETA2),1) 501 WRITE (LUPRI,*) 'Norm of THETA2 after second E contribution:', 502 & XNORM 503 Call AROUND('final result for B{O} matrix transformation:') 504 Call CC_PRP(WORK(KTHETA1),WORK(KTHETA2),ISYRES,1,1) 505 CALL FLSHFO(LUPRI) 506 END IF 507 508 END IF ! (.NOT. (CCS .OR. CCSTST)) 509 510*---------------------------------------------------------------------* 511* initialize the R12 doubles part of the result vector THETA: 512*---------------------------------------------------------------------* 513 IF ( CCR12 ) THEN 514 CALL DZERO(WORK(KTHETAR12),NTR12AM(ISYRES)) 515 516 KLAMDPB = KEND2 517 KLAMDHB = KLAMDPB + NGLMDT(ISYMTB) 518 KLAMDPA = KLAMDHB + NGLMDT(ISYMTB) 519 KLAMDHA = KLAMDPA + NGLMDT(ISYMTA) 520 KT12AMP = KLAMDHA + NGLMDT(ISYMTA) 521 KXINTTRI= KT12AMP + MAX(NTR12SQ(ISYMTA),NTR12SQ(ISYMTB)) 522 KXINTSQ = KXINTTRI + MAX(NR12R12P(1),NTR12SQ(ISYRES)) 523 KSCR = KXINTSQ + NR12R12SQ(1) 524 KEND3 = KSCR + NTR12SQ(ISYRES) 525 LEND3 = LWORK - KEND3 526 IF (LEND3 .LT. 0) THEN 527 CALL QUIT('Insufficient work space in CC_BAMAT') 528 END IF 529 530 CALL CCPRPAO(LABELO,.TRUE.,WORK(KPERT),IRREP,ISYM,IERR, 531 & WORK(KEND2),LEND2) 532 IF (IERR.NE.0 .OR. IRREP.NE.ISYMO) THEN 533 CALL QUIT('CC_BAMAT: error while reading operator '//LABELO) 534 END IF 535 536 CALL CCLR_LAMTRA(WORK(KLAMDP0),WORK(KLAMDPA), WORK(KLAMDH0), 537 & WORK(KLAMDHA),WORK(KT1AMPA),ISYMTA) 538 539 CALL CCLR_LAMTRA(WORK(KLAMDP0),WORK(KLAMDPB), WORK(KLAMDH0), 540 & WORK(KLAMDHB),WORK(KT1AMPB),ISYMTB) 541 542 LUNIT = -1 543 CALL GPOPEN(LUNIT,FCCR12X,'OLD',' ','UNFORMATTED', 544 & DUMMY,.FALSE.) 545 REWIND(LUNIT) 546 8888 READ(LUNIT) IAN 547 READ(LUNIT) (WORK(KXINTTRI-1+I), I=1, NR12R12P(1)) 548 IF (IAN.NE.IANR12) GOTO 8888 549 CALL GPCLOSE(LUNIT,'KEEP') 550 IOPT = 2 551 CALL CCR12UNPCK2(WORK(KXINTTRI),1,WORK(KXINTSQ),'N',IOPT) 552C 553 DO I = 1, 2 554 IF (I.EQ.1) THEN 555 ISYM1 = ISYMTA 556 ISYM2 = ISYMTB 557 LIST1 = LISTA 558 IDLST1 = IDLSTA 559 ELSE IF (I.EQ.2) THEN 560 ISYM1 = ISYMTB 561 ISYM2 = ISYMTA 562 LIST1 = LISTB 563 IDLST1 = IDLSTB 564 END IF 565 !read R12 response amplitudes 566 CALL CC_R12GETCT(WORK(KT12AMP),ISYM1,2,KETSCL,.FALSE.,'N', 567 & DUMMY,DUMMY,DUMMY,LIST1,IDLST1,WORK(KEND3),LEND3) 568 !calculate.... 569 IF (I.EQ.1) THEN 570 CALL CC_R12XI2A(WORK(KSCR),ISYRES,WORK(KT12AMP),ISYM1, 571 & WORK(KPERT),ISYMO,WORK(KLAMDP0), 572 & WORK(KLAMDHB),ISYM2,'N', 573 & WORK(KEND3),LEND3) 574 CALL DCOPY(NTR12SQ(ISYRES),WORK(KSCR),1,WORK(KXINTTRI),1) 575 ELSE IF (I.EQ.2) THEN 576 CALL CC_R12XI2A(WORK(KSCR),ISYRES,WORK(KT12AMP),ISYM1, 577 & WORK(KPERT),ISYMO,WORK(KLAMDP0), 578 & WORK(KLAMDHA),ISYM2,'N', 579 & WORK(KEND3),LEND3) 580 END IF 581 END DO 582C 583 CALL DAXPY(NTR12SQ(ISYRES),ONE,WORK(KXINTTRI),1,WORK(KSCR),1) 584 CALL DZERO(WORK(KXINTTRI),NTR12SQ(ISYRES)) 585 CALL CC_R12XI2B(WORK(KXINTTRI),'N',WORK(KXINTSQ),1,'N', 586 & WORK(KSCR),ISYRES,'N',-ONE) 587C 588 IOPT = 1 589 CALL CCR12PCK2(WORK(KTHETAR12),ISYRES,.FALSE.,WORK(KXINTTRI), 590 & 'N',IOPT) 591 CALL CCLR_DIASCLR12(WORK(KTHETAR12),BRASCL,ISYRES) 592C 593 END IF 594 595*---------------------------------------------------------------------* 596* compute CC3 contributions: 597*---------------------------------------------------------------------* 598 IF (CCSDT) THEN 599 600 IF (IOPTRES.EQ.5) THEN 601 FREQB = 0.0D0 602 ELSE 603 FREQB = FREQLST(FILBAM,IBATRAN(4,ITRAN)) 604 END IF 605 606 CALL CCSDT_BAMAT_NODDY(IOPTRES,FREQB,LABELO,ISYMO, 607 & LISTA,IDLSTA, 608 & LISTB,IDLSTB, 609 & WORK(KTHETA1),WORK(KTHETA2), 610 & WORK(KTHETA1EFF),WORK(KTHETA2EFF), 611 & IBDOTS,BCONS,FILBAM,ITRAN, 612 & NBATRAN,MXVEC,WORK(KEND2),LEND2) 613 614 END IF 615 616*---------------------------------------------------------------------* 617* write result vector to output: 618*---------------------------------------------------------------------* 619 IF (IOPTRES .EQ. 0 .OR. IOPTRES .EQ. 1) THEN 620 621* write to a common direct access file, 622* store start address in IBATRAN(4,ITRAN) 623 624 IBATRAN(4,ITRAN) = IADRTH 625 626 CALL PUTWA2(LUBMAT,FILBAM,WORK(KTHETA1),IADRTH,NT1AM(ISYRES)) 627 IADRTH = IADRTH + NT1AM(ISYRES) 628 629 IF (.NOT.CCS) THEN 630 CALL PUTWA2(LUBMAT,FILBAM,WORK(KTHETA2),IADRTH,NT2AM(ISYRES)) 631 IADRTH = IADRTH + NT2AM(ISYRES) 632 END IF 633 634 IF (CCR12) THEN 635 CALL PUTWA2(LUBMAT,FILBAM,WORK(KTHETAR12),IADRTH, 636 & NTR12AM(ISYRES)) 637 IADRTH = IADRTH + NTR12AM(ISYRES) 638 END IF 639 640 IF (LOCDBG) THEN 641 WRITE (LUPRI,*) 'B{O} matrix transform. nb. ',ITRAN, 642 & ' saved on file.' 643 WRITE (LUPRI,*) 'ADRESS, LENGTH:', 644 & IBATRAN(4,ITRAN),IADRTH-IBATRAN(4,ITRAN) 645 XNORM = DDOT(NT1AM(ISYRES),WORK(KTHETA1),1,WORK(KTHETA1),1) 646 IF (.NOT.CCS) XNORM = XNORM + 647 & DDOT(NT2AM(ISYRES),WORK(KTHETA2),1,WORK(KTHETA2),1) 648 IF (CCR12) XNORM = XNORM + 649 & DDOT(NTR12AM(ISYRES),WORK(KTHETAR12),1,WORK(KTHETAR12),1) 650 WRITE (LUPRI,*) 'Norm:', XNORM 651 652 Call AROUND('B{O} matrix transformation written to file:') 653 Call CC_PRP(WORK(KTHETA1),WORK(KTHETA2),ISYRES,1,1) 654 IF (CCR12) CALL CC_PRPR12(WORK(KTHETAR12),ISYRES,1,.TRUE.) 655 END IF 656 657 ELSE IF ( IOPTRES .EQ. 3 .OR. IOPTRES .EQ. 4 ) THEN 658 659* write to a sequential file by a call to CC_WRRSP/CC_WARSP, 660* use FILBAM as LIST type and IBATRAN(4,ITRAN) as index 661 662 KTHETA0 = -999999 663 664 IF (IOPTRES.EQ.3) THEN 665 CALL CC_WRRSP(FILBAM,IBATRAN(4,ITRAN),ISYRES,IOPTW,MODELW, 666 & WORK(KTHETA0),WORK(KTHETA1),WORK(KTHETA2), 667 & WORK(KEND2),LEND2) 668 IF (CCR12) THEN 669 CALL CC_WRRSP(FILBAM,IBATRAN(4,ITRAN),ISYRES,IOPTWR12, 670 & MODELW,DUMMY,DUMMY,WORK(KTHETAR12), 671 & WORK(KEND2),LEND2) 672 END IF 673 IF (CCSDT) THEN 674 CALL CC_WRRSP(FILBAM,IBATRAN(4,ITRAN),ISYRES,IOPTWE,MODEL, 675 & WORK(KTHETA0),WORK(KTHETA1EFF), 676 & WORK(KTHETA2EFF),WORK(KEND2),LEND2) 677 END IF 678 ELSE IF (IOPTRES.EQ.4) THEN 679 CALL CC_WARSP(FILBAM,IBATRAN(4,ITRAN),ISYRES,IOPTW,MODELW, 680 & WORK(KTHETA0),WORK(KTHETA1),WORK(KTHETA2), 681 & WORK(KEND2),LEND2) 682 IF (CCR12) THEN 683 CALL CC_WARSP(FILBAM,IBATRAN(4,ITRAN),ISYRES,IOPTWR12, 684 & MODELW,DUMMY,DUMMY,WORK(KTHETAR12), 685 & WORK(KEND2),LEND2) 686 END IF 687 IF (CCSDT) THEN 688 CALL CC_WARSP(FILBAM,IBATRAN(4,ITRAN),ISYRES,IOPTWE,MODELW, 689 & WORK(KTHETA0),WORK(KTHETA1EFF), 690 & WORK(KTHETA2EFF),WORK(KEND2),LEND2) 691 END IF 692 END IF 693 694 IF (LOCDBG) THEN 695 WRITE (LUPRI,*) 'Write B{O} * ',LISTA,' * ',LISTB, 696 & ' transformation', 697 & ' as ',FILBAM,' type vector to file.' 698 WRITE (LUPRI,*) 'index of inp. O integr:',IBATRAN(1,ITRAN) 699 WRITE (LUPRI,*) 'index of inp. A vector:',IBATRAN(2,ITRAN) 700 WRITE (LUPRI,*) 'index of inp. B vector:',IBATRAN(3,ITRAN) 701 WRITE (LUPRI,*) 'index of result vector:',IBATRAN(4,ITRAN) 702 XNORM = DDOT(NT1AM(ISYRES),WORK(KTHETA1),1,WORK(KTHETA1),1) 703 IF (.NOT.CCS) XNORM = XNORM + 704 & DDOT(NT2AM(ISYRES),WORK(KTHETA2),1,WORK(KTHETA2),1) 705 IF (CCR12) XNORM = XNORM + 706 & DDOT(NTR12AM(ISYRES),WORK(KTHETAR12),1,WORK(KTHETAR12),1) 707 WRITE (LUPRI,*) 'norm^2 of result vector:',XNORM 708 END IF 709 710 ELSE IF ( IOPTRES .EQ. 5 ) THEN 711 712 CALL CCDOTRSP(IBDOTS,BCONS,IOPTW,FILBAM,ITRAN,NBATRAN,MXVEC, 713 & WORK(KTHETA1),WORK(KTHETA2),ISYRES, 714 & WORK(KEND2),LEND2) 715 IF (CCR12) THEN 716 CALL CCDOTRSP(IBDOTS,BCONS,IOPTWR12,FILBAM,ITRAN,NBATRAN, 717 & MXVEC,DUMMY,WORK(KTHETAR12),ISYRES, 718 & WORK(KEND2),LEND2) 719 END IF 720 721 ELSE 722 CALL QUIT('Illegal value for IOPTRES in CC_BAMAT.') 723 END IF 724*---------------------------------------------------------------------* 725* End of loop over B matrix transformations 726*---------------------------------------------------------------------* 727 END DO 728 729*---------------------------------------------------------------------* 730* if IOPTRES=1 and enough work space available, read result 731* vectors back into memory: 732*---------------------------------------------------------------------* 733 734* check size of work space: 735 IF (IOPTRES .EQ. 1) THEN 736 LENALL = IADRTH-1 737 IF (LENALL .GT. LWORK) IOPTRES = 0 738 END IF 739 740* read the result vectors back into memory: 741 IF (IOPTRES .EQ. 1) THEN 742 743 CALL GETWA2(LUBMAT,FILBAM,WORK(1),1,LENALL) 744 745 IF (LOCDBG) THEN 746 DO ITRAN = 1, NBATRAN 747 IF (ITRAN.LT.NBATRAN) THEN 748 LEN = IBATRAN(4,ITRAN+1)-IBATRAN(4,ITRAN) 749 ELSE 750 LEN = IADRTH-IBATRAN(4,NBATRAN) 751 END IF 752 KTHETA1 = IBATRAN(4,ITRAN) 753 XNORM = DDOT(LEN, WORK(KTHETA1),1, WORK(KTHETA1),1) 754 WRITE (LUPRI,*) 'Read B matrix transformation nb. ',NBATRAN 755 WRITE (LUPRI,*) 'Adress, length, NORM:',IBATRAN(4,NBATRAN), 756 & LEN,XNORM 757 END DO 758 CALL FLSHFO(LUPRI) 759 END IF 760 END IF 761 762*---------------------------------------------------------------------* 763* close B matrix file & return 764*---------------------------------------------------------------------* 765* check return option for the result vectors: 766 IF (IOPTRES .EQ. 0 .OR. IOPTRES .EQ. 1) THEN 767 CALL WCLOSE2(LUBMAT,FILBAM,'DELETE') 768 769 ELSE IF (IOPTRES.EQ.3 .OR. IOPTRES.EQ.4 .OR. IOPTRES.EQ.5) THEN 770 CONTINUE 771 ELSE 772 CALL QUIT('Illegal value of IOPTRES in CC_BAMAT.') 773 END IF 774 775 776*=====================================================================* 777 778 RETURN 779 END 780*=====================================================================* 781* END OF SUBROUTINE CC_BAMAT 782*=====================================================================* 783