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_FMATRIX */ 21*=====================================================================* 22 SUBROUTINE CC_FMATRIX(IFTRAN, NFTRAN, LISTL, LISTR, IOPTRES, 23 & FILFMA, IFDOTS, FCONS, MXVEC, WORK, LWORK) 24*---------------------------------------------------------------------* 25* 26* Purpose: batched loop over F matrix transformations 27* (needed if the number of transformations exceeds the 28* limit MAXSIM defined on ccsdio.h ) 29* 30* for more detailed documentation see: CC_FMAT 31* 32* Written by Christof Haettig, November 1998, based on CC_BMATRIX. 33* 34*=====================================================================* 35#if defined (IMPLICIT_NONE) 36 IMPLICIT NONE 37#else 38# include "implicit.h" 39#endif 40#include "priunit.h" 41#include "ccsdinp.h" 42#include "maxorb.h" 43Cholesky 44#include "ccdeco.h" 45Cholesky 46#include "ccsdio.h" 47#include "ccnoddy.h" 48 49 LOGICAL LOCDBG 50 PARAMETER (LOCDBG = .FALSE.) 51 52 CHARACTER*(*) LISTL, LISTR, FILFMA 53 CHARACTER*8 FN3VI2, FNTOC 54 CHARACTER*4 FNDKBC 55 CHARACTER*6 FNDELD, FNCKJD, FN3VI 56 INTEGER IOPTRES 57 INTEGER NFTRAN, MXVEC, LWORK 58 INTEGER IFTRAN(3,NFTRAN) 59 INTEGER IFDOTS(MXVEC,NFTRAN) 60 61#if defined (SYS_CRAY) 62 REAL WORK(LWORK) 63 REAL FCONS(MXVEC,NFTRAN) 64#else 65 DOUBLE PRECISION WORK(LWORK) 66 DOUBLE PRECISION FCONS(MXVEC,NFTRAN) 67COMMENT COMMENT 68 DOUBLE PRECISION DDOT 69COMMENT COMMENT 70#endif 71 72 LOGICAL LADD 73 INTEGER MAXFTRAN, NTRAN, ISTART, IBATCH, NBATCH 74 INTEGER KIADRBFD, KIADRZ0, KIADRPQ0, KIADRPQR, KIADRPQMO 75 INTEGER KIADRBFI, KRHO1, KRHO2, KFIM, KMGDL, KRIM, KZDEN, KGIM 76 INTEGER KLAMPA, KLAMHA, KCTR2, KLAMP, KLAMH, KDENS, KFOCK 77 INTEGER KEND, LEND 78 INTEGER MXRAMP, MXLAMP, MXBTRAN, KB1TRAN, KB1DOTS, KB1CON, 79 & KB2TRAN, KB2DOTS, KB2CON, NB1TRAN, NB2TRAN 80 81 CALL QENTER('CC_FMATRIX') 82 83 IF (IOPTRES.EQ.1) 84 & CALL QUIT('IOPTRES cannnot be used with CC_FMATRIX') 85 86*---------------------------------------------------------------------* 87* Main section: CC_FMATNEW, driven by a loop over transformed vectors 88* ------------- 89* singles and doubles models: 90* calculation of the single and double excitation 91* parts of the transformed vectors and respective 92* dot products (IOPTRES=5). 93* 94* triple models: 95* calculation of the effective transformed vectors, 96* for IOPTRES=5 only singles and doubles contributions 97* and contributions from <L3|...|HF> are computed here 98* the remaining contributions <L2|[...,R3]|HF> are 99* calculated in CCSDT_FMAT_CONT (see below). 100*---------------------------------------------------------------------* 101* 102Cholesky 103C 104 IF (CHOINT .AND. CC2) THEN 105 106 CALL CC_CHOFMAT(IFTRAN, NFTRAN, LISTL, LISTR, IOPTRES, 107 & FILFMA, IFDOTS, FCONS, MXVEC, WORK, LWORK) 108 109 GOTO 9999 110 END IF 111C 112Cholesky 113* 114 MAXFTRAN = MAXSIM 115 116 NBATCH = (NFTRAN+MAXFTRAN-1)/MAXFTRAN 117 118 IF (LOCDBG) THEN 119 WRITE (LUPRI,*) 'Batching over F matrix transformations:' 120 WRITE (LUPRI,*) 'nb. of batches needed:', NBATCH 121 END IF 122 123 DO IBATCH = 1, NBATCH 124 ISTART = (IBATCH-1) * MAXFTRAN + 1 125 NTRAN = MIN(NFTRAN-(ISTART-1),MAXFTRAN) 126 127 KIADRBFD = 1 128 KIADRZ0 = KIADRBFD + MXCORB_CC*NTRAN 129 KIADRPQ0 = KIADRZ0 + MXCORB_CC*NTRAN 130 KIADRPQR = KIADRPQ0 + MXCORB_CC*NTRAN 131 KIADRPQMO = KIADRPQR + MXCORB_CC*NTRAN 132 KIADRBFI = KIADRPQMO + MXCORB_CC*NTRAN 133 KRHO1 = KIADRBFI + MXCORB_CC*NTRAN 134 KRHO2 = KRHO1 + NTRAN 135 KFIM = KRHO2 + NTRAN 136 KMGDL = KFIM + NTRAN 137 KRIM = KMGDL + NTRAN 138 KZDEN = KRIM + NTRAN 139 KGIM = KZDEN + NTRAN 140 KLAMPA = KGIM + NTRAN 141 KLAMHA = KLAMPA + NTRAN 142 KCTR2 = KLAMHA + NTRAN 143 KLAMP = KCTR2 + NTRAN 144 KLAMH = KLAMP + NTRAN 145 KDENS = KLAMH + NTRAN 146 KFOCK = KDENS + NTRAN 147 KEND = KFOCK + NTRAN 148 LEND = LWORK - KEND 149 150 IF (LEND .LT. 0) THEN 151 WRITE (LUPRI,*) 'Insufficient work space in CC_FMATRIX.' 152 WRITE (LUPRI,*) 'Available :',LWORK,' words,' 153 WRITE (LUPRI,*) 'Need at least:',KEND, ' words.' 154 CALL QUIT('Insufficient work space in CC_FMATRIX.') 155 END IF 156 157 IF (LOCDBG) THEN 158 WRITE (LUPRI,*) 'Batch No.:',IBATCH 159 WRITE (LUPRI,*) 'start at :',ISTART 160 WRITE (LUPRI,*) '# transf.:',NTRAN 161 END IF 162 163 CALL CC_FMATNEW(IFTRAN(1,ISTART), NTRAN, 164 & LISTL, LISTR, IOPTRES, FILFMA, 165 & IFDOTS(1,ISTART), FCONS(1,ISTART), 166 & WORK(KIADRBFD), WORK(KIADRZ0), WORK(KIADRPQ0), 167 & WORK(KIADRPQR),WORK(KIADRPQMO),WORK(KIADRBFI), 168 & WORK(KRHO1),WORK(KRHO2),WORK(KFIM),WORK(KMGDL), 169 & WORK(KRIM),WORK(KZDEN),WORK(KGIM),WORK(KLAMPA), 170 & WORK(KLAMHA), WORK(KCTR2), WORK(KLAMP), 171 & WORK(KLAMH), WORK(KDENS), WORK(KFOCK), 172 & MXVEC, WORK(KEND), LEND ) 173 174 END DO 175*---------------------------------------------------------------------* 176* special triples section: 177* ------------------------ 178* compute the contributions <L2|[...,R3]|HF> to 179* F matrix contractions 180*---------------------------------------------------------------------* 181 IF (IOPTRES.EQ.5 .AND. CCSDT .AND. CCSDT_F_ALTER) THEN 182 MXRAMP = MAX(NFTRAN,MXVEC) 183 MXLAMP = NFTRAN 184 MXBTRAN = MXRAMP*MXRAMP 185 186 KB1TRAN = 1 187 KB1DOTS = KB1TRAN + MXBTRAN * 3 188 KB1CON = KB1DOTS + MXBTRAN * MXLAMP 189 KEND = KB1CON + MXBTRAN * MXLAMP 190 191 IF (LISTR(1:3).NE.FILFMA(1:3)) THEN 192 KB2TRAN = KEND 193 KB2DOTS = KB2TRAN + MXBTRAN * 3 194 KB2CON = KB2DOTS + MXBTRAN * MXLAMP 195 KEND = KB2CON + MXBTRAN * MXLAMP 196 ELSE 197 KB2TRAN = -999999 198 KB2DOTS = -999999 199 KB2CON = -999999 200 END IF 201 202 LEND = LWORK - KEND 203 IF (LEND.LE.0) CALL QUIT('Out of memory in CC_FMATRIX.') 204 205 CALL DZERO(WORK(KB1TRAN),MXBTRAN*3) 206 CALL DZERO(WORK(KB1DOTS),MXBTRAN*MXLAMP) 207 CALL DZERO(WORK(KB1CON), MXBTRAN*MXLAMP) 208 IF (LISTR(1:3).NE.FILFMA(1:3)) THEN 209 CALL DZERO(WORK(KB2TRAN),MXBTRAN*3) 210 CALL DZERO(WORK(KB2DOTS),MXBTRAN*MXLAMP) 211 CALL DZERO(WORK(KB2CON), MXBTRAN*MXLAMP) 212 END IF 213 214 LADD = .FALSE. 215 CALL CCSDT_F2B_SETUP(IFTRAN,IFDOTS,FCONS,NFTRAN,MXVEC,LADD, 216 & WORK(KB1TRAN),WORK(KB1DOTS), 217 & WORK(KB1CON),NB1TRAN,MXLAMP, 218 & WORK(KB2TRAN),WORK(KB2DOTS), 219 & WORK(KB2CON),NB2TRAN,MXLAMP, 220 & LISTL,LISTR,FILFMA,MXBTRAN,MXBTRAN) 221 222 IF (NODDY_FMAT) THEN 223 CALL CCSDT_FBC_NODDY(CCSDT_F_ALTER, 224 & WORK(KB1TRAN),WORK(KB1DOTS), 225 & WORK(KB1CON),NB1TRAN,MXLAMP, 226 & LISTL,LISTR,FILFMA,WORK(KEND),LEND) 227 ELSE 228 FNCKJD = 'CKJDEL' 229 FNDKBC = 'DKBC' 230 FNDELD = 'CKDELD' 231 FNTOC = 'CCSDT_OC' 232 FN3VI = 'CC3_VI' 233 FN3VI2 = 'CC3_VI12' 234 235 IF (LISTR(1:3).EQ.'R1 '.OR.LISTR(1:3).EQ.'RE ' 236 * .OR.LISTR(1:3).EQ.'RC ') THEN 237 CALL CCSDT_FBMAT(WORK(KB1TRAN),WORK(KB1DOTS), 238 * WORK(KB1CON),NB1TRAN,MXLAMP, 239 * LISTL,LISTR,FILFMA,FNCKJD,FNDKBC, 240 * FNDELD,FNTOC,FN3VI,FN3VI2, 241 * WORK(KEND),LEND) 242 ELSE IF (LISTR(1:3).EQ.'R2 '.OR.LISTR(1:3).EQ.'ER1') THEN 243 CALL CC3_FBMATT3ZU(WORK(KB1TRAN),WORK(KB1DOTS), 244 * WORK(KB1CON),NB1TRAN,MXLAMP, 245 * LISTL,LISTR,FILFMA, 246 * WORK(KEND),LEND) 247 ELSE 248 CALL QUIT('Unknown LISTR in CC_FMAT.') 249 END IF 250 251 END IF 252 253 IF (LISTR(1:3).NE.FILFMA(1:3)) THEN 254 IF (NODDY_FMAT) THEN 255 CALL CCSDT_FBC_NODDY(CCSDT_F_ALTER, 256 & WORK(KB2TRAN),WORK(KB2DOTS), 257 & WORK(KB2CON),NB2TRAN,MXLAMP, 258 & LISTL,FILFMA,LISTR,WORK(KEND),LEND) 259 ELSE 260 FNCKJD = 'CKJDEL' 261 FNDKBC = 'DKBC' 262 FNDELD = 'CKDELD' 263 FNTOC = 'CCSDT_OC' 264 FN3VI = 'CC3_VI' 265 FN3VI2 = 'CC3_VI12' 266 267 IF (FILFMA(1:3).EQ.'R1 '.OR.FILFMA(1:3).EQ.'RE ' 268 * .OR.FILFMA(1:3).EQ.'RC ') THEN 269 CALL CCSDT_FBMAT(WORK(KB2TRAN),WORK(KB2DOTS), 270 * WORK(KB2CON),NB2TRAN,MXLAMP, 271 * LISTL,FILFMA,LISTR,FNCKJD,FNDKBC, 272 * FNDELD,FNTOC,FN3VI,FN3VI2, 273 * WORK(KEND),LEND) 274 ELSE IF (FILFMA(1:3).EQ.'R2 '.OR.FILFMA(1:3).EQ.'ER1') THEN 275 CALL CC3_FBMATT3ZU(WORK(KB2TRAN),WORK(KB2DOTS), 276 * WORK(KB2CON),NB2TRAN,MXLAMP, 277 * LISTL,FILFMA,LISTR, 278 * WORK(KEND),LEND) 279 ELSE 280 CALL QUIT('Unknown FILFMA in CC_FMAT.') 281 END IF 282 283 END IF 284 END IF 285 286 LADD = .TRUE. 287 CALL CCSDT_F2B_SETUP(IFTRAN,IFDOTS,FCONS,NFTRAN,MXVEC,LADD, 288 & WORK(KB1TRAN),WORK(KB1DOTS), 289 & WORK(KB1CON),NB1TRAN,MXLAMP, 290 & WORK(KB2TRAN),WORK(KB2DOTS), 291 & WORK(KB2CON),NB2TRAN,MXLAMP, 292 & LISTL,LISTR,FILFMA,MXBTRAN,MXBTRAN) 293 294 END IF 295 296 IF (LOCDBG) WRITE(LUPRI,*) 'Leave CC_FMATRIX...' 297 298 9999 CONTINUE 299 CALL QEXIT('CC_FMATRIX') 300 301 RETURN 302 END 303 304*---------------------------------------------------------------------* 305* END OF SUBROUTINE CC_FMATRIX * 306*---------------------------------------------------------------------* 307*---------------------------------------------------------------------* 308c/* Deck CC_FTRAN */ 309*=====================================================================* 310 SUBROUTINE CC_FTRAN(LISTL, IDLSTL, LISTR, IDLSTR, WORK, LWORK) 311*---------------------------------------------------------------------* 312* 313* Purpose: shortcut to F matrix transformation for a single 314* transformation (should be avoided and maybe completely 315* eliminated) 316* 317* Written by Christof Haettig, Januar 1999, based on CC_FMATRIX. 318* 319*=====================================================================* 320#if defined (IMPLICIT_NONE) 321 IMPLICIT NONE 322#else 323# include "implicit.h" 324#endif 325#include "priunit.h" 326#include "maxorb.h" 327#include "ccorb.h" 328#include "ccsdsym.h" 329#include "ccsdinp.h" 330 331 INTEGER LUFMAT 332 333 CHARACTER*(*) LISTL, LISTR 334 CHARACTER*(8) FILFMA 335 INTEGER IOPTRES, IDLSTL, IDLSTR, IDUM, NFTRAN, LWORK 336 INTEGER IFTRAN(3,1) 337 INTEGER KIADRBFD, KIADRZ0, KIADRPQ0, KIADRPQR, KIADRPQMO 338 INTEGER KIADRBFI, KRHO1, KRHO2, KFIM, KMGDL, KRIM, KZDEN, KGIM 339 INTEGER KLAMPA, KLAMHA, KCTR2, KLAMP, KLAMH, KDENS, KFOCK 340 INTEGER KEND, LEND, ISYML, ISYMR, ISYRES, LENALL 341 342 INTEGER ILSTSYM 343 344 REAL*8 WORK(LWORK), RDUM 345 346 CALL QENTER('CC_FTRAN') 347 348 IF (CCSDT.OR.CC3) CALL 349 & QUIT('Triples not carried through in CC_FTRAN!') 350 351 352 NFTRAN = 1 353 IFTRAN(1,1) = IDLSTL 354 IFTRAN(2,1) = IDLSTR 355 IOPTRES = 0 356 FILFMA = 'CC__FMAT' 357 358 KIADRBFD = 1 359 KIADRZ0 = KIADRBFD + MXCORB_CC*NFTRAN 360 KIADRPQ0 = KIADRZ0 + MXCORB_CC*NFTRAN 361 KIADRPQR = KIADRPQ0 + MXCORB_CC*NFTRAN 362 KIADRPQMO = KIADRPQR + MXCORB_CC*NFTRAN 363 KIADRBFI = KIADRPQMO + MXCORB_CC*NFTRAN 364 KRHO1 = KIADRBFI + MXCORB_CC*NFTRAN 365 KRHO2 = KRHO1 + NFTRAN 366 KFIM = KRHO2 + NFTRAN 367 KMGDL = KFIM + NFTRAN 368 KRIM = KMGDL + NFTRAN 369 KZDEN = KRIM + NFTRAN 370 KGIM = KZDEN + NFTRAN 371 KLAMPA = KGIM + NFTRAN 372 KLAMHA = KLAMPA + NFTRAN 373 KCTR2 = KLAMHA + NFTRAN 374 KLAMP = KCTR2 + NFTRAN 375 KLAMH = KLAMP + NFTRAN 376 KDENS = KLAMH + NFTRAN 377 KFOCK = KDENS + NFTRAN 378 KEND = KFOCK + NFTRAN 379 LEND = LWORK - KEND 380 381 IF (LEND .LT. 0) THEN 382 WRITE (LUPRI,*) 'Insufficient work space in CC_FTRAN.' 383 WRITE (LUPRI,*) 'Available :',LWORK,' words,' 384 WRITE (LUPRI,*) 'Need at least:',KEND, ' words.' 385 CALL QUIT('Insufficient work space in CC_FTRAN.') 386 END IF 387 388 CALL CC_FMATNEW( IFTRAN, NFTRAN, LISTL, LISTR, IOPTRES, FILFMA, 389 & IDUM, RDUM, 390 & WORK(KIADRBFD), WORK(KIADRZ0), WORK(KIADRPQ0), 391 & WORK(KIADRPQR),WORK(KIADRPQMO),WORK(KIADRBFI), 392 & WORK(KRHO1),WORK(KRHO2),WORK(KFIM),WORK(KMGDL), 393 & WORK(KRIM),WORK(KZDEN),WORK(KGIM),WORK(KLAMPA), 394 & WORK(KLAMHA), WORK(KCTR2), WORK(KLAMP), 395 & WORK(KLAMH), WORK(KDENS), WORK(KFOCK), 396 & 0, WORK(KEND), LEND ) 397 398 ISYML = ILSTSYM(LISTL,IDLSTL) 399 ISYMR = ILSTSYM(LISTR,IDLSTR) 400 ISYRES = MULD2H(ISYML,ISYMR) 401 402 LENALL = NT1AM(ISYRES) + NT2AM(ISYRES) 403 IF (CCS) LENALL = NT1AM(ISYRES) 404 IF (LENALL.GT.LWORK) THEN 405 WRITE (LUPRI,*) 'Insufficient work space in CC_FTRAN.' 406 WRITE (LUPRI,*) 'Available :',LWORK,' words,' 407 WRITE (LUPRI,*) 'Need at least:',LENALL, ' words.' 408 CALL QUIT('Insufficient work space in CC_FTRAN.') 409 END IF 410 411 LUFMAT = -1 412 CALL WOPEN2(LUFMAT, FILFMA, 64, 0) 413 CALL GETWA2(LUFMAT,FILFMA,WORK(1),1,LENALL) 414 CALL WCLOSE2(LUFMAT, FILFMA, 'DELETE') 415 416 CALL QEXIT('CC_FTRAN') 417 418 RETURN 419 END 420 421*---------------------------------------------------------------------* 422* END OF SUBROUTINE CC_FTRAN * 423*---------------------------------------------------------------------* 424*---------------------------------------------------------------------* 425c/* Deck CC_FMATNEW */ 426*=====================================================================* 427 SUBROUTINE CC_FMATNEW(IFTRAN, NFTRAN, LISTL, LISTR, IOPTRES, 428 & FILFMA, IFDOTS, FCONS, 429 & IADRBFD, IADRZ0, IADRPQ0, IADRPQR, 430 & IADRPQMO, IADRBFI, KRHO1, KRHO2, KFIM, 431 & KMGDL, KRIM, KZDEN, KGIM, KLAMPA, 432 & KLAMHA, KCTR2, KLAMP, KLAMH, 433 & KDENS, KFOCK, MXVEC, WORK, LWORK ) 434*---------------------------------------------------------------------* 435* 436* Purpose: AO-direct calculation of a linear transformation of two 437* CC amplitude vectors, L and R, with the CC F matrix 438* (or B matrix for L different from zeroth-order multipl.) 439* 440* L.EQ.'L0' ---> F * R 441* L.NE.'L0' ---> L * B * R 442* 443* the linear transformations are calculated for a list 444* of L vectors and a list of R vectors: 445* 446* LISTL -- type of L vectors 447* LISTR -- type of R vectors 448* IFTRAN(1,*) -- indeces of L vectors 449* IFTRAN(2,*) -- indeces of R vectors 450* IFTRAN(3,*) -- indeces or addresses of result vectors 451* NFTRAN -- number of requested transformations 452* FILFMA -- file name / list type of result vectors 453* or list type of vectors to be dotted on 454* IFDOTS -- indeces of vectors to be dotted on 455* FCONS -- contains the dot products on return 456* 457* return of the result vectors: 458* 459* IOPTRES = 0 : all result vectors are written to a direct 460* access file, FILFMA is used as file name 461* the start addresses of the vectors are 462* returned in IFTRAN(3,*) 463* 464* IOPTRES = 1 : the vectors are kept and returned in WORK 465* if possible, start addresses returned in 466* IFTRAN(3,*). N.B.: if WORK is not large 467* enough IOPTRES is automatically reset to 0!! 468* 469* IOPTRES = 3 : each result vector is written to its own 470* file by a call to CC_WRRSP, FILFMA is used 471* as list type and IFTRAN(3,*) as list index 472* NOTE that IFTRAN(3,*) is in this case input! 473* 474* IOPTRES = 4 : each result vector is added to a vector on 475* file by a call to CC_WARSP, FILFMA is used 476* as list type and IFTRAN(3,*) as list index 477* NOTE that IFTRAN(3,*) is in this case input! 478* 479* IOPTRES = 5 : the result vectors are dotted on a array 480* of vectors, the type of the arrays given 481* by FILFMA and the indeces from IFDOTS 482* the result of the dot products is returned 483* in the FCONS array 484* 485* Written by Christof Haettig, November 1998, based on CC_FMAT. 486* 487*=====================================================================* 488 USE PELIB_INTERFACE, ONLY: USE_PELIB, PELIB_IFC_TRANSFORMER, 489 & PELIB_IFC_QRTRANSFORMER 490#if defined (IMPLICIT_NONE) 491 IMPLICIT NONE 492#else 493# include "implicit.h" 494#endif 495#include "priunit.h" 496#include "ccsdinp.h" 497#include "ccsections.h" 498#include "ccsdsym.h" 499#include "maxorb.h" 500#include "mxcent.h" 501#include "ccsdio.h" 502#include "ccorb.h" 503#include "cciccset.h" 504#include "cbieri.h" 505#include "distcl.h" 506#include "iratdef.h" 507#include "eritap.h" 508#include "ccisao.h" 509#include "ccfield.h" 510#include "aovec.h" 511#include "blocks.h" 512#include "ccslvinf.h" 513#include "ccnoddy.h" 514#include "second.h" 515#include "dummy.h" 516#include "r12int.h" 517#include "qm3.h" 518!#include "qmmm.h" 519 520* local parameters: 521 CHARACTER MSGDBG*(17) 522 PARAMETER (MSGDBG='[debug] CC_FMAT> ') 523 524 LOGICAL LOCDBG 525 PARAMETER (LOCDBG = .FALSE.) 526 527 LOGICAL APPEND, NOAPPEND 528 PARAMETER (APPEND = .TRUE., NOAPPEND = .FALSE.) 529 530 INTEGER KDUM 531 PARAMETER( KDUM = +99 999 999 ) ! dummy address for work space 532 INTEGER ISYM0 533 PARAMETER( ISYM0 = 1 ) ! symmetry of the reference state 534 535 INTEGER LUBF, LUBFD, LUC, LUD, LUFK, LURIM, LUBFI 536 INTEGER LURHO, LUFIM, LUMGD, LUXIM, LUYIM, LUBDZ0, LUO3 537 INTEGER LUCBAR, LUDBAR, LUPQRR, LUPQR0, LUPQMO, LUFMAT, LUZDN 538 INTEGER LUGIM, LUMIM 539 INTEGER LUCKJD, LUDKBC, LUTOC, LU3VI, LUDKBC3, LU3FOP, LU3FOP2 540 INTEGER LU3VI2, LU3FOPX, LU3FOP2X 541 542 CHARACTER*(8) BFFIL, FNBFD, CTFIL, DTFIL, FILFCK, FILRIM, FILO3 543 CHARACTER*(8) RHOFIL,FIMFIL,FILMGD,FILXIM,FILYIM,FNBFI,FILZDN 544 CHARACTER*(8) CBAFIL, DBAFIL, FNPQRR, FNPQR0, FNPQMO, FNBDZ0 545 CHARACTER*(8) GIMFIL, FNTOC, FN3VI2, FILMIM 546 CHARACTER*(7) FN3FOP2X 547 CHARACTER*(6) FNCKJD, FN3VI, FN3FOP2, FNDELD, FN3FOPX 548 CHARACTER*(5) FNDKBC3, FN3FOP 549 CHARACTER*(4) FNDKBC 550 CHARACTER*(1) RSPTYP 551 PARAMETER (BFFIL ='CCFM_BFI', FNBFD ='CCBFDENS', 552 & CBAFIL='CCFM_CBA', DBAFIL='CCFM_DBA', 553 & CTFIL ='CCFM_CIM', DTFIL ='CCFM_DIM', 554 & FIMFIL='CCFM_FIM', RHOFIL='CCFM_RHO', 555 & FILRIM='CCFM_RIM', FILMGD='CCFM_MGD', 556 & FILXIM='CCFM_XIM', FILYIM='CCFM_YIM', 557 & FILMIM='CCFM_MIM', 558 & FILFCK='CCFM_FCK', FNPQRR='CCFMPQRR', 559 & FNPQR0='CCFMPQR0', FNPQMO='CCFMPQMO', 560 & FNBDZ0='CCFMBDZ0', FNBFI ='CCFM_BZI', 561 & FILO3 ='CCFM_O3I', FILZDN='CCFM_ZDN', 562 & GIMFIL='CCFM_GIM', FNCKJD='CKJDEL', 563 & FNDKBC='DKBC' , FNTOC ='CCSDT_OC', 564 & FN3VI ='CC3_VI' , FNDKBC3='DKBC3', 565 & FN3FOP='PTFOP' , FN3FOP2='PTFOP2', 566 & FN3FOPX='PTFOPX' , FN3FOP2X='PTFOP2X', 567 & FNDELD='CKDELD' , FN3VI2 = 'CC3_VI12' ) 568 569 570 CHARACTER*(*) LISTL, LISTR, FILFMA 571 INTEGER IOPTRES, IOPTTCME, IDUM 572 INTEGER NFTRAN, MXVEC, LWORK 573 INTEGER IFTRAN(3,NFTRAN) 574 INTEGER IFDOTS(MXVEC,NFTRAN) 575 576 577#if defined (SYS_CRAY) 578 REAL WORK(LWORK) 579 REAL FCONS(MXVEC,NFTRAN) 580 REAL DUM, XNORM, FF, DUMMY 581 REAL DTIME, CONVRT, TIMALL 582 REAL TIMFCK, TIMBF, TIMC, TIMD, TIMIMR, TIMXYM, TIMIO 583 REAL TIMPQ, TIMIML, TIMIM2, TIMPRE, TIMINT, TIMRDAO 584 REAL TIMTRBT, TIMIM0, TIMTRN, TIMGZ, TIMFG 585 REAL ZERO, ONE, TWO, HALF, FACT 586C REAL ECURRSAV 587#else 588 DOUBLE PRECISION WORK(LWORK) 589 DOUBLE PRECISION FCONS(MXVEC,NFTRAN) 590 DOUBLE PRECISION XNORM, FF 591 DOUBLE PRECISION DTIME, CONVRT, TIMALL 592 DOUBLE PRECISION TIMFCK, TIMBF, TIMC, TIMD, TIMIMR, TIMXYM, TIMIO 593 DOUBLE PRECISION TIMPQ, TIMIML, TIMIM2, TIMPRE, TIMINT, TIMRDAO 594 DOUBLE PRECISION TIMTRBT, TIMIM0, TIMTRN, TIMGZ, TIMFG 595 DOUBLE PRECISION ZERO, ONE, TWO, HALF, FACT 596C DOUBLE PRECISION ECURRSAV 597#endif 598 PARAMETER (ZERO = 0.0d0, ONE = 1.0d0, TWO = 2.0d0, HALF = 0.5d0) 599 600 CHARACTER*(10) MODEL, MODELW, CDUMMY 601 CHARACTER*(3) APROXR12 602 LOGICAL LLSAME, CCSDR12, OSQSAV, OORSAV 603 INTEGER INDEXA(MXCORB_CC) 604 INTEGER INTMEDR(2,MAXSIM), NINTR, IINTR 605 INTEGER INTMEDL(4,MAXSIM), NINTL, IINTL 606 INTEGER INTMED2(4,MAXSIM), NINT2, IINT2 607 INTEGER IRHGH(0:MAXSIM), I2HGH(0:MAXSIM), IOFFCD(0:MAXSIM+1) 608 INTEGER IADRBFD(MXCORB_CC,NFTRAN), IADRZ0(MXCORB_CC,NFTRAN) 609 INTEGER IADRPQ0(MXCORB_CC,NFTRAN), IADRPQR(MXCORB_CC,NFTRAN) 610 INTEGER IADRPQMO(MXCORB_CC,NFTRAN),IADRBFI(MXCORB_CC,NFTRAN) 611 INTEGER KRHO1(NFTRAN),KFIM(NFTRAN),KMGDL(NFTRAN),KRIM(NFTRAN) 612 INTEGER KZDEN(NFTRAN),KGIM(NFTRAN) 613 INTEGER KLAMPA(NFTRAN),KLAMHA(NFTRAN),KRHO2(NFTRAN),KCTR2(NFTRAN) 614 INTEGER KLAMP(NFTRAN),KLAMH(NFTRAN),KDENS(NFTRAN),KFOCK(NFTRAN) 615 INTEGER LENX,LENY,LENM,LENMGD,LENR,LENFIM,LENRHO,LENFK,LENBF,LEN 616 INTEGER MT2BGD,MDISAO,MDSRHF,MSCRATCH,MEMAVAIL,NWORK,NSECMAR 617 INTEGER NNWORK, ICDEL2, IBATCH, KEND2, JEND2, LWRK2 618 INTEGER IOPTW, ITRAN, IOPT, ICORE, IF, IOPTRSP, ISTARTBFI, IVEC 619 INTEGER ISYM, ICOUNT, ISYMAK, ISYBET, IBSRHF(8,8), NBSRHF(8) 620 INTEGER ISYML, ISYMR, ISYRES, IDLSTL, IDLSTR, KEND1, JEND1 621 INTEGER KFOCK0, KDENS0, KT1AMP0, KDSRHFA, ISYMXA, KFCKC0, KDNSC0 622 INTEGER KLAMP0, KLAMH0, KTHETA1, KTHETA2, KZETA1, KT1AMPA 623 INTEGER KCCFB1, KINDXB, KODCL1, KODCL2, KODBC1, KODBC2, NBATCH 624 INTEGER KRDBC1, KRDBC2, KODPP1, KODPP2, KRDPP1, KRDPP2, KFREE 625 INTEGER LFREE, KEND, LWRK, KENDSV, LWRKSV, KEND0, LWRK0, LWRK1 626 INTEGER NTOSYM, NTOT, KRECNR, ISYMD1, ILLL, NUMDIS, IDEL2, IDEL 627 INTEGER ISYDEL, KXINT, KEND3, K3OINT, KBFRHF, KDCRHF, KDSRHF 628 INTEGER KFINT, KRIMA, KFOCKA, KTHETA0, KYINT, KYBARA, LENZDN 629 INTEGER LWRK3, IADRTH, KLAMDPA, KLAMDHA, LENALL, KT2AMPA 630 INTEGER ISYCTR, ISYAMP, KLIAJB, KA2IM, KXIDJL, KXJLID 631 INTEGER KXIAJB, NDBLE, LENGIM, KGZETA, KFOCK0AO, KXBARA 632 INTEGER KTHETA1EFF, KTHETA2EFF, IOPTWE 633 INTEGER KTHETAR12,LENMOD,IOPTWR12,KSCR1,KSCR2,ICON,IAMP 634 INTEGER KMINT 635 636* external functions: 637 INTEGER ICCSET1 638 INTEGER ICCSET2 639 INTEGER ILSTSYM 640 641 REAL*8, ALLOCATABLE :: FOCKMAT(:), FOCKTEMP(:) 642 643#if defined (SYS_CRAY) 644 REAL DDOT 645#else 646 DOUBLE PRECISION DDOT 647#endif 648 649 CALL QENTER('CC_FMATNEW') 650 651*---------------------------------------------------------------------* 652* begin: 653*---------------------------------------------------------------------* 654 IF (LOCDBG) THEN 655 Call AROUND('ENTERED CC_FMAT') 656 IF (DIRECT) WRITE(LUPRI,'(/1X,A)') 'AO direct transformation' 657 WRITE (LUPRI,*) MSGDBG, 'LISTL : ',LISTL(1:3) 658 WRITE (LUPRI,*) MSGDBG, 'LISTR : ',LISTR(1:3) 659 WRITE (LUPRI,*) MSGDBG, 'FILFMA: ',FILFMA(1:3) 660 WRITE (LUPRI,*) MSGDBG, 'NFTRAN: ',NFTRAN 661 WRITE (LUPRI,*) MSGDBG, 'IOPTRES:',IOPTRES 662 CALL FLSHFO(LUPRI) 663 END IF 664 665 IF ( .NOT. (CCS .OR. CC2 .OR. CCSD .OR. CC3) ) THEN 666 WRITE (LUPRI,'(/1x,2a)')'CC_FMAT called for a Coupled Cluster ', 667 & 'method not implemented in CC_FMAT...' 668 CALL QUIT('Unknown CC method in CC_FMAT.') 669 END IF 670 671 IF (.NOT. DUMPCD) THEN 672 WRITE (LUPRI,*) 'DUMPCD = ',DUMPCD 673 WRITE (LUPRI,*) 'CC_FMAT requires DUMPCD=.TRUE.' 674 CALL QUIT('DUMPCD=.FALSE. , CC_FMAT requires DUMPCD=.TRUE.') 675 END IF 676 677 IF (NFTRAN .GT. MAXSIM) THEN 678 WRITE (LUPRI,*) 679 * 'Error in CC_FMAT: NFTRAN is larger than MAXSIM.' 680 CALL QUIT('Error in CC_FMAT: NFTRAN is larger than MAXSIM.') 681 END IF 682 683 IF (IPRINT.GT.2) THEN 684 685 WRITE (LUPRI,'(//1X,A1,50("="),A1)') '+','+' 686 687 WRITE (LUPRI,'(1x,A52)') 688 & '| F MATRIX TRANSFORMATION SECTION |' 689 690 IF (IOPTRES.EQ.0) THEN 691 WRITE (LUPRI,'(1X,A52)') 692 & '| (results are written to a direct access file) |' 693 ELSE IF (IOPTRES.EQ.1) THEN 694 WRITE (LUPRI,'(1X,A52)') 695 & '| (results are returned in memory) |' 696 ELSE IF (IOPTRES.EQ.3) THEN 697 WRITE (LUPRI,'(1X,A52)') 698 & '| (result is written to file) |' 699 ELSE IF (IOPTRES.EQ.4) THEN 700 WRITE (LUPRI,'(1X,A52)') 701 & '| (result is added to a vector on file) |' 702 ELSE IF (IOPTRES.EQ.5) THEN 703 WRITE (LUPRI,'(1X,A52)') 704 & '| (result used to calculate dot products) |' 705 END IF 706 707 WRITE (LUPRI,'(1X,A1,50("-"),A1)') '+','+' 708 709 END IF 710 711* initialize timings: 712 TIMALL = SECOND() 713 TIMIM0 = ZERO 714 TIMIML = ZERO 715 TIMIMR = ZERO 716 TIMIM2 = ZERO 717 TIMPRE = ZERO 718 TIMXYM = ZERO 719 TIMPQ = ZERO 720 TIMFCK = ZERO 721 TIMBF = ZERO 722 TIMC = ZERO 723 TIMD = ZERO 724 TIMFG = ZERO 725 TIMGZ = ZERO 726 TIMINT = ZERO 727 TIMRDAO = ZERO 728 TIMTRBT = ZERO 729 TIMIO = ZERO 730 731* set option and model to write vectors to file: 732 IOPTW = 0 733 IOPTWR12 = 0 734 IF (CCS) THEN 735 MODELW = 'CCS ' 736 IOPTW = 1 737 ELSE IF (CC2) THEN 738 MODELW = 'CC2 ' 739 IOPTW = 3 740 ELSE IF (CCSD) THEN 741 MODELW = 'CCSD ' 742 IOPTW = 3 743 ELSE IF (CC3) THEN 744 MODELW = 'CC3 ' 745 IOPTW = 3 746 IOPTWE = 24 747 ELSE 748 CALL QUIT('Unknown coupled cluster model in CC_FMAT.') 749 END IF 750 IF (CCR12) THEN 751 APROXR12 = ' ' 752 CALL CCSD_MODEL(MODELW,LENMOD,10,MODELW,10,APROXR12) 753 IOPTWR12 = 32 754 END IF 755 756 !set CCSDR12 757 CCSDR12 = CCR12 .AND. .NOT.(CC2.OR.CCS) 758 759* check return option for the result vectors: 760 LUFMAT = -1 761 IF (IOPTRES .EQ. 0 .OR. IOPTRES .EQ. 1) THEN 762 CALL WOPEN2(LUFMAT, FILFMA, 64, 0) 763 ELSE IF (IOPTRES .EQ. 3 .OR. IOPTRES .EQ. 4) THEN 764 CONTINUE 765 ELSE IF (IOPTRES .EQ. 5) THEN 766 IF (MXVEC*NFTRAN.NE.0) CALL DZERO(FCONS,MXVEC*NFTRAN) 767 ELSE 768 CALL QUIT('Illegal value of IOPTRES in CC_FMAT.') 769 END IF 770 771* precalculate symmetry array for BSRHF: 772 DO ISYM = 1, NSYM 773 ICOUNT = 0 774 DO ISYMAK = 1, NSYM 775 ISYBET = MULD2H(ISYMAK,ISYM) 776 IBSRHF(ISYMAK,ISYBET) = ICOUNT 777 ICOUNT = ICOUNT + NT1AO(ISYMAK)*NBAS(ISYBET) 778 END DO 779 NBSRHF(ISYM) = ICOUNT 780 END DO 781 782*=====================================================================* 783* build nonredundant arrays of response vectors and pairs of them 784* for which intermediates have to be calculated 785*=====================================================================* 786 DTIME = SECOND() 787 788* array for intermediates which depend on the left vectors: 789 NINTL = 0 790 DO ITRAN = 1, NFTRAN 791 I=ICCSET2(INTMEDL,LISTL,IFTRAN(1,ITRAN), 792 & 'R0 ',0, NINTL,MAXSIM,APPEND) 793 END DO 794 795* array for intermediates which depend on the right vectors: 796 NINTR = 0 797 DO ITRAN = 1, NFTRAN 798 I=ICCSET1(INTMEDR,LISTR,IFTRAN(2,ITRAN),NINTR,MAXSIM,APPEND) 799 END DO 800 801* array for intermediates that depend on left and right vectors: 802 NINT2 = 0 803 DO ITRAN = 1, NFTRAN 804 I=ICCSET2(INTMED2,LISTL,IFTRAN(1,ITRAN), 805 & LISTR,IFTRAN(2,ITRAN),NINT2,MAXSIM,APPEND) 806 END DO 807 808 809 IF (LOCDBG) THEN 810 WRITE (LUPRI,'(/A)') 811 & 'List of response vector for left intermediates:' 812 WRITE (LUPRI,'((/5X,2I5))') ((INTMEDL(I,J),I=1,2),J=1,NINTL) 813 WRITE (LUPRI,'(/A)') 814 & 'List of response vector for right intermediates:' 815 WRITE (LUPRI,'((/5X,2I5))') ((INTMEDR(I,J),I=1,2),J=1,NINTR) 816 WRITE (LUPRI,'(/A)') 817 & 'List of vector pairs for 2. order intermediates:' 818 WRITE (LUPRI,'((/5X,4I5))') ((INTMED2(I,J),I=1,4),J=1,NINT2) 819 END IF 820 821 TIMPRE = TIMPRE + SECOND() - DTIME 822*---------------------------------------------------------------------* 823* estimate scratch space requirements 824*---------------------------------------------------------------------* 825 DTIME = SECOND() 826 827 MT2BGD = 0 828 MDISAO = 0 829 MDSRHF = 0 830 DO ISYM = 1, NSYM 831 MT2BGD = MAX(MT2BGD,NT2BGD(ISYM)) 832 MDISAO = MAX(MDISAO,NDISAO(ISYM)) 833 MDSRHF = MAX(MDSRHF,NDSRHF(ISYM)) 834 END DO 835 836* 5 x a NT2BGD type intermediate 837* + integral arrays + some reserve 838 839 MSCRATCH = 5*MT2BGD + MDISAO + MDSRHF + 10*N2BASX 840 IF (CC2) MSCRATCH = MSCRATCH+MDISAO+2*MDSRHF 841 IF (.NOT.(CCS.OR.CC2)) MSCRATCH = MSCRATCH+MDISAO+4*MDSRHF 842 843 IF (LOCDBG) THEN 844 WRITE (LUPRI,*) MSGDBG, 'scratch space estimate MSCRATCH:', 845 & MSCRATCH 846 CALL FLSHFO(LUPRI) 847 END IF 848 849 TIMPRE = TIMPRE + SECOND() - DTIME 850*---------------------------------------------------------------------* 851* estimate memory for 'in core' version and batched versions: 852*---------------------------------------------------------------------* 853 DTIME = SECOND() 854 855 MEMAVAIL = LWORK - MSCRATCH 856 857 NWORK = 0 858 NBATCH = 1 859 IF (CCS) THEN 860 NSECMAR = 10 * N2BASX 861 ELSE IF (CC2 .OR. CCSD .OR. CCSDT) THEN 862 NSECMAR = 10 * MT2BGD 863 ELSE 864 CALL QUIT('Unknown CC model in CC_FMAT.') 865 END IF 866 867 IRHGH(0) = 0 868 I2HGH(0) = 0 869 870* intermediates that dependent on right response vectors: 871* (see routine CCBPRE1 for details) 872 DO IINTR = 1, NINTR 873 IDLSTR = INTMEDR(1,IINTR) 874 ISYM = ILSTSYM(LISTR,IDLSTR) 875 876 NNWORK = 2*NGLMDT(ISYM) + 2*N2BST(ISYM) 877 IF (CCSD .OR. CCSDT) THEN 878 NNWORK = 2*NGLMDT(ISYM)+NT2AOIJ(ISYM)+NEMAT1(ISYM) 879 END IF 880 881 IF( (NWORK+NNWORK+NSECMAR).GT.MEMAVAIL ) THEN 882 IRHGH(NBATCH) = IINTR - 1 883 I2HGH(NBATCH) = 0 884 885 NBATCH = NBATCH + 1 886 NWORK = 0 887 END IF 888 NWORK = NWORK + NNWORK 889 IF (NWORK .GT. LWORK) THEN 890 WRITE (LUPRI,*) 'Insufficient work space in CC_FMAT. (01)' 891 WRITE (LUPRI,*) 'Need at least:',NNWORK, ' words.' 892 CALL FLSHFO(LUPRI) 893 CALL QUIT('Insufficient work space in CC_FMAT. (01)') 894 END IF 895 END DO 896 897* intermediates that dependent on left and right vectors: 898* (see routine CCFPRE2 for details) 899 DO IINT2 = 1, NINT2 900 IDLSTL = INTMED2(1,IINT2) 901 IDLSTR = INTMED2(3,IINT2) 902 ISYML = ILSTSYM(LISTL,IDLSTL) 903 ISYMR = ILSTSYM(LISTR,IDLSTR) 904 ISYRES = MULD2H(ISYML,ISYMR) 905 906 IF (CCS) THEN 907 NNWORK = 2*N2BST(ISYRES) 908 ELSE IF (CC2) THEN 909 NNWORK = 2*N2BST(ISYRES) + 2*NGLMDT(ISYMR) + 910 & NT2SQ(ISYML) + NT1AM(ISYRES) 911 ELSE IF (CCSD .OR. CCSDT) THEN 912 NNWORK = 2*NGLMDT(ISYMR) + 913 & NT1AO(ISYRES) + NT1AM(ISYRES) + NT2AOIJ(ISYRES) 914 ELSE 915 CALL QUIT('Unknown CC model in CC_FMAT.') 916 END IF 917 918 IF( (NWORK+NNWORK+NSECMAR).GT.MEMAVAIL ) THEN 919 IRHGH(NBATCH) = NINTR 920 I2HGH(NBATCH) = IINT2 - 1 921 922 NBATCH = NBATCH + 1 923 NWORK = 0 924 END IF 925 NWORK = NWORK + NNWORK 926 IF (NWORK .GT. LWORK) THEN 927 WRITE (LUPRI,*) 'Insufficient work space in CC_FMAT. (02)' 928 WRITE (LUPRI,*) 'Need at least:',NNWORK,' words.' 929 CALL FLSHFO(LUPRI) 930 CALL QUIT('Insufficient work space in CC_FMAT. (02)') 931 END IF 932 END DO 933 934 IRHGH(NBATCH) = NINTR 935 I2HGH(NBATCH) = NINT2 936 937 IF (LOCDBG .AND. (NBATCH.EQ.1)) THEN 938 WRITE (LUPRI,*) MSGDBG,'one batch only... will be done in core.' 939 WRITE (LUPRI,*) MSGDBG, 'memory for intermediates: ', NWORK 940 WRITE (LUPRI,*) MSGDBG, 'remaining scratch space: ', LWORK-NWORK 941 CALL FLSHFO(LUPRI) 942 ELSE IF (LOCDBG .AND. (NBATCH.GT.1)) THEN 943 WRITE (LUPRI,*) MSGDBG, 944 & 'more than one batch... choose I/O algorithm.' 945 WRITE (LUPRI,*) MSGDBG, 'max. memory for intermediates: ', 946 & MEMAVAIL 947 WRITE (LUPRI,*) MSGDBG, 'number of batches: ',NBATCH 948 CALL FLSHFO(LUPRI) 949 END IF 950 951 TIMPRE = TIMPRE + SECOND() - DTIME 952*---------------------------------------------------------------------* 953* read zeroth-order singles amplitudes, allocate space for Fock matrix, 954* and prepare zeroth-order lambda matrices and density: 955*---------------------------------------------------------------------* 956 DTIME = SECOND() 957 958 KFOCK0AO = 1 959 KFOCK0 = KFOCK0AO + N2BAST 960 961 KDENS0 = KFOCK0 + N2BAST 962 KT1AMP0 = KDENS0 + N2BAST 963 KLAMP0 = KT1AMP0 + NT1AMX 964 KLAMH0 = KLAMP0 + NLAMDT 965 KEND0 = KLAMH0 + NLAMDT 966 967 IF (FROIMP.OR.FROEXP) THEN 968 KFCKC0 = KEND0 969 KDNSC0 = KFCKC0 + N2BAST 970 KEND0 = KDNSC0 + N2BAST 971 END IF 972 973 LWRK0 = LWORK - KEND0 974 975 IF (LWRK0 .LT. 0) THEN 976 CALL QUIT('Insufficient work space in CC_FMAT. (0)') 977 END IF 978 979* read zeroth order amplitudes: 980 IOPT = 1 981 Call CC_RDRSP('R0',0,ISYM0,IOPT,MODEL,WORK(KT1AMP0),WORK(KDUM)) 982 983* get unperturbed Lambda matrices: 984 Call LAMMAT(WORK(KLAMP0),WORK(KLAMH0),WORK(KT1AMP0), 985 & WORK(KEND0),LWRK0) 986 987* calculate the density matrix: 988 ICORE = 1 ! include core contribution 989 CALL CC_AODENS(WORK(KLAMP0),WORK(KLAMH0),WORK(KDENS0), 990 & ISYM0,ICORE, WORK(KEND0),LWRK0) 991 992* calculate pure core contribution to the density matrix, 993* and initialize core contribution to Fock matrix with zeros 994 IF (FROIMP.OR.FROEXP) THEN 995 ICORE = 0 ! exclude core contribution 996 CALL CC_AODENS(WORK(KLAMP0),WORK(KLAMH0),WORK(KDNSC0), 997 & ISYM0,ICORE, WORK(KEND0),LWRK0) 998 CALL DSCAL(N2BAST,-ONE,WORK(KDNSC0),1) 999 CALL DAXPY(N2BAST,+ONE,WORK(KDENS0),1,WORK(KDNSC0),1) 1000 CALL DZERO(WORK(KFCKC0),N2BAST) 1001 END IF 1002 1003* initialize Fock matrix with the one-electron integrals: 1004 CALL CCRHS_ONEAO(WORK(KFOCK0),WORK(KEND0),LWRK0) 1005 DO IF= 1, NFIELD 1006 FF = EFIELD(IF) 1007 CALL CC_ONEP(WORK(KFOCK0),WORK(KEND0),LWRK0,FF,1,LFIELD(IF) ) 1008 END DO 1009 1010* Solvent contribution to one-electron integrals. 1011* T^g contribution to transformation. SLV98,OC, CCMM JK, NYQMMM KS 1012 IF (CCSLV) THEN 1013 IF (.NOT.CCMM) CALL CCSL_RHSTG(WORK(KFOCK0),WORK(KEND0),LWRK0) 1014 IF (CCMM) THEN 1015 IF (.NOT. NYQMMM) THEN 1016 CALL CCMM_RHSTG(WORK(KFOCK0),WORK(KEND0),LWRK0) 1017 ELSE IF (NYQMMM) THEN 1018 IF (HFFLD) THEN 1019 CALL CCMM_ADDGHF(WORK(KFOCK0),WORK(KEND0),LWRK0) 1020 ELSE 1021 CALL CCMM_ADDG(WORK(KFOCK0),WORK(KEND0),LWRK0) 1022 END IF 1023 END IF 1024 END IF 1025 ENDIF 1026 IF (USE_PELIB()) THEN 1027 ALLOCATE(FOCKMAT(NNBASX),FOCKTEMP(N2BAST)) 1028 IF (HFFLD) THEN 1029 CALL GET_FROM_FILE('FOCKMHF',NNBASX,FOCKMAT) 1030 ELSE 1031 CALL GET_FROM_FILE('FOCKMAT',NNBASX,FOCKMAT) 1032 END IF 1033 CALL DSPTSI(NBAS,FOCKMAT,FOCKTEMP) 1034 CALL DAXPY(N2BAST,1.0d0,FOCKTEMP,1,WORK(KFOCK0),1) 1035 DEALLOCATE(FOCKMAT,FOCKTEMP) 1036 END IF 1037C 1038C==================================================== 1039C 1040 IF (LOCDBG) THEN 1041 WRITE (LUPRI,*) MSGDBG, 'norm of T1AMP0:', 1042 & DDOT(NT1AMX,WORK(KT1AMP0),1,WORK(KT1AMP0),1) 1043 WRITE (LUPRI,*) MSGDBG, 'norm of XLAMP0:', 1044 & DDOT(NLAMDT,WORK(KLAMP0),1,WORK(KLAMP0),1) 1045 WRITE (LUPRI,*) MSGDBG, 'norm of XLAMH0:', 1046 & DDOT(NLAMDT,WORK(KLAMH0),1,WORK(KLAMH0),1) 1047 WRITE (LUPRI,*) MSGDBG, 'norm of DENS0:', 1048 & DDOT(N2BAST,WORK(KDENS0),1,WORK(KDENS0),1) 1049 WRITE (LUPRI,*) MSGDBG, 'norm of FOCK0:', 1050 & DDOT(N2BAST,WORK(KFOCK0),1,WORK(KFOCK0),1) 1051 END IF 1052 1053 TIMPRE = TIMPRE + SECOND() - DTIME 1054 TIMIM0 = TIMIM0 + SECOND() - DTIME 1055 1056*---------------------------------------------------------------------* 1057* precalculate some intermediates which depend only on the left 1058* vectors: X, Y, M, Chi, Yps, P, Q, and the BFZeta density 1059*---------------------------------------------------------------------* 1060 DTIME = SECOND() 1061 1062 IOPTRSP = 1 1063 CALL CCFPREINT(INTMEDL,NINTL,WORK(KLAMP0),WORK(KLAMH0), 1064 & CDUMMY,LENX,CDUMMY,LENY,CDUMMY,LENM, 1065 & CDUMMY,LENMGD,CDUMMY,LENZDN, 1066 & FNPQR0,IADRPQ0,FNPQMO,IADRPQMO,FNBDZ0,IADRZ0, 1067 & TIMXYM,TIMBF,TIMIO,TIMPQ,IOPTRSP,WORK(KEND0),LWRK0) 1068 1069 TIMIML = TIMIML + SECOND() - DTIME 1070 1071*---------------------------------------------------------------------* 1072* precalculate some intermediates which depend only on the right 1073* vectors: CBAR, DBAR, and the BF density 1074*---------------------------------------------------------------------* 1075 DTIME = SECOND() 1076 1077 CALL CCBPRE1INT(INTMEDR,NINTR,IOFFCD,IADRBFD, 1078 & CBAFIL,DBAFIL,FNBFD, 1079 & WORK(KLAMP0),WORK(KLAMH0), 1080 & TIMIO,TIMC,TIMD,TIMBF,WORK(KEND0),LWRK0) 1081 1082 TIMIMR = TIMIMR + SECOND() - DTIME 1083 1084*---------------------------------------------------------------------* 1085* precalculate some intermediates which depend on left and right 1086* vectors: response X, Y, M, Chi, Yps, P, Q, and BFZeta density 1087*---------------------------------------------------------------------* 1088 DTIME = SECOND() 1089 1090 IOPTRSP = 2 1091 CALL CCFPREINT(INTMED2,NINT2,WORK(KLAMP0),WORK(KLAMH0), 1092 & FILXIM,LENX,FILYIM,LENY,FILMIM,LENM, 1093 & FILMGD,LENMGD,FILZDN,LENZDN, 1094 & FNPQRR,IADRPQR,CDUMMY,IDUMMY,CDUMMY,IDUMMY, 1095 & TIMXYM,TIMBF,TIMIO,TIMPQ,IOPTRSP,WORK(KEND0),LWRK0) 1096 1097 TIMIM2 = TIMIM2 + SECOND() - DTIME 1098 1099*---------------------------------------------------------------------* 1100* open files for local intermediates generated in the loop over the 1101* AO integrals and which need to be initialized here: 1102*---------------------------------------------------------------------* 1103 DTIME = SECOND() 1104 1105 CALL CCFOPEN(LUBF, LUFK, LURIM, LUFIM, LUGIM, LURHO, 1106 & BFFIL, FILFCK, FILRIM, FIMFIL, GIMFIL, RHOFIL, 1107 & LENBF, LENFK, LENR, LENFIM, LENGIM, LENRHO, 1108 & NINTR, NINT2, WORK(KEND0), LWRK0 ) 1109 1110* open some other files needed in the loop over integrals, but need 1111* not to be initialized: 1112 LUBFD = -1 1113 LUZDN = -1 1114 LUC = -1 1115 LUD = -1 1116 LUMGD = -1 1117 LUBFI = -1 1118 LUBFD = -1 1119 LUPQR0 = -1 1120 LUPQRR = -1 1121 LUBDZ0 = -1 1122 LUO3 = -1 1123 IF (CCS.OR.CC2) THEN 1124 CALL WOPEN2(LUZDN, FILZDN, 64,0) 1125 ELSE IF (.NOT.(CCS.OR.CC2)) THEN 1126 CALL WOPEN2(LUC, CTFIL, 64,0) 1127 CALL WOPEN2(LUD, DTFIL, 64,0) 1128 CALL WOPEN2(LUMGD, FILMGD,64,0) 1129 CALL WOPEN2(LUBFI, FNBFI, 64,0) 1130 CALL WOPEN2(LUBFD, FNBFD, 64,0) 1131 CALL WOPEN2(LUPQR0,FNPQR0,64,0) 1132 CALL WOPEN2(LUPQRR,FNPQRR,64,0) 1133 CALL WOPEN2(LUBDZ0,FNBDZ0,64,0) 1134 CALL WOPEN2(LUO3 , FILO3, 64,0) 1135 END IF 1136 1137* initialize offsets C and D intermediates: 1138 ICDEL2 = 0 1139 1140* initialize start address for BFI intermediates: 1141 ISTARTBFI = 1 1142 1143 TIMPRE = TIMPRE + SECOND() - DTIME 1144*---------------------------------------------------------------------* 1145* if all vectors and intermediates fit into the memory, read all 1146* response vectors before the loop over AO integral shells: 1147*---------------------------------------------------------------------* 1148 DTIME = SECOND() 1149 1150 IF (NBATCH .EQ. 1) THEN 1151 1152 CALL CCBPRE1(INTMEDR, 1, NINTR, 1153 & KRHO2, KLAMP, KLAMH, KDENS, KFOCK, KRIM, 1154 & LUBF,BFFIL,LENBF,LUFK,FILFCK,LENFK, 1155 & LURIM,FILRIM,LENR, 1156 & WORK(KLAMP0), WORK(KLAMH0), 1157 & WORK, LWORK, KEND0, JEND1 ) 1158 KEND1 = JEND1 1159 1160 CALL CCFPRE2(INTMED2,1,NINT2, 1161 & KFIM, LUFIM,FIMFIL,LENFIM, 1162 & KGIM, LUGIM,GIMFIL,LENGIM, 1163 & KRHO1,LURHO,RHOFIL,LENRHO, 1164 & KMGDL,LUMGD,FILMGD,LENMGD, 1165 & KZDEN,LUZDN,FILZDN,LENZDN, 1166 & KCTR2,KLAMPA,KLAMHA, 1167 & WORK(KLAMP0),WORK(KLAMH0), 1168 & WORK, LWORK, KEND1, JEND1 ) 1169 KEND1 = JEND1 1170 1171 IF (LOCDBG) THEN 1172 WRITE (LUPRI,*) MSGDBG, 1173 & 'allocated work space for intermediates:' 1174 WRITE (LUPRI,*) MSGDBG,'KRHO2 :',(KRHO2(I),I=1,NINTR) 1175 WRITE (LUPRI,*) MSGDBG,'KLAMP :',(KLAMP(I),I=1,NINTR) 1176 WRITE (LUPRI,*) MSGDBG,'KLAMH :',(KLAMH(I),I=1,NINTR) 1177 WRITE (LUPRI,*) MSGDBG,'KDENS :',(KDENS(I),I=1,NINTR) 1178 WRITE (LUPRI,*) MSGDBG,'KFOCK :',(KFOCK(I),I=1,NINTR) 1179 WRITE (LUPRI,*) MSGDBG,'KRIM :',(KRIM(I),I=1,NINTR) 1180 WRITE (LUPRI,*) MSGDBG,'KRHO1 :',(KRHO1(I),I=1,NINT2) 1181 WRITE (LUPRI,*) MSGDBG,'KZDEN :',(KZDEN(I),I=1,NINT2) 1182 WRITE (LUPRI,*) MSGDBG,'KFIM :',(KFIM(I),I=1,NINT2) 1183 WRITE (LUPRI,*) MSGDBG,'KGIM :',(KGIM(I),I=1,NINT2) 1184 WRITE (LUPRI,*) MSGDBG,'KMGDL :',(KMGDL(I),I=1,NINT2) 1185 WRITE (LUPRI,*) MSGDBG,'KLAMPA:',(KLAMPA(I),I=1,NINT2) 1186 WRITE (LUPRI,*) MSGDBG,'KLAMHA:',(KLAMHA(I),I=1,NINT2) 1187 WRITE (LUPRI,*) MSGDBG,'KEND1:',KEND1 1188 CALL FLSHFO(LUPRI) 1189 END IF 1190 1191 ELSE 1192 KEND1 = KEND0 1193 END IF 1194 1195 LWRK1 = LWORK - KEND1 1196 IF (LWRK1 .LT. 0) THEN 1197 CALL QUIT('Insufficient work space in CC_FMAT. (1)') 1198 END IF 1199 1200 TIMPRE = TIMPRE + SECOND() - DTIME 1201*---------------------------------------------------------------------* 1202* initialize integral calculation 1203*---------------------------------------------------------------------* 1204 DTIME = SECOND() 1205 1206 KEND = KEND1 1207 LWRK = LWRK1 1208 1209 IF (DIRECT) THEN 1210 NTOSYM = 1 1211 1212 IF (HERDIR) THEN 1213 CALL HERDI1(WORK(KEND),LWRK,IPRERI) 1214 ELSE 1215 KCCFB1 = KEND 1216 KINDXB = KCCFB1 + MXPRIM*MXCONT 1217 KEND = KINDXB + (8*MXSHEL*MXCONT + 1)/IRAT 1218 LWRK = LWORK - KEND 1219 CALL ERIDI1(KODCL1,KODCL2,KODBC1,KODBC2,KRDBC1,KRDBC2, 1220 * KODPP1,KODPP2,KRDPP1,KRDPP2, 1221 * KFREE,LFREE,KEND,WORK(KCCFB1),WORK(KINDXB), 1222 * WORK(KEND),LWRK,IPRERI) 1223 1224 KEND = KFREE 1225 LWRK = LFREE 1226 END IF 1227 1228 KENDSV = KEND 1229 LWRKSV = LWRK 1230 ELSE 1231 NTOSYM = NSYM 1232 END IF 1233 1234 TIMINT = TIMINT + SECOND() - DTIME 1235*---------------------------------------------------------------------* 1236* start loop over AO integral shells: 1237*---------------------------------------------------------------------* 1238 DO ISYMD1 = 1, NTOSYM 1239 1240 IF (DIRECT) THEN 1241 IF (HERDIR) THEN 1242 NTOT = MAXSHL 1243 ELSE 1244 NTOT = MXCALL 1245 ENDIF 1246 ELSE 1247 NTOT = NBAS(ISYMD1) 1248 END IF 1249 1250 DO ILLL = 1, NTOT 1251 1252 DTIME = SECOND() 1253 1254 IF (DIRECT) THEN 1255 KEND = KENDSV 1256 LWRK = LWRKSV 1257 1258 IF (HERDIR) THEN 1259 CALL HERDI2(WORK(KEND),LWRK,INDEXA,ILLL,NUMDIS, 1260 & IPRINT) 1261 ELSE 1262 CALL ERIDI2(ILLL,INDEXA,NUMDIS,0,0, 1263 * WORK(KODCL1),WORK(KODCL2), 1264 * WORK(KODBC1),WORK(KODBC2), 1265 * WORK(KRDBC1),WORK(KRDBC2), 1266 * WORK(KODPP1),WORK(KODPP2), 1267 * WORK(KRDPP1),WORK(KRDPP2), 1268 * WORK(KCCFB1),WORK(KINDXB), 1269 * WORK(KEND),LWRK,IPRERI) 1270 END IF 1271 1272 KRECNR = KEND 1273 KEND = KRECNR + (NBUFX(0) - 1)/IRAT + 1 1274 LWRK = LWORK - KEND 1275 1276 IF (LWRK .LT. 0) THEN 1277 CALL QUIT('Insufficient work space in CC_FMAT. (1a)') 1278 END IF 1279 1280 ELSE 1281 NUMDIS = 1 1282 END IF 1283 1284 TIMINT = TIMINT + SECOND() - DTIME 1285*---------------------------------------------------------------------* 1286* if out of core: allocate memory and get response vectors: 1287*---------------------------------------------------------------------* 1288 DO IBATCH = 1, NBATCH 1289 KEND2 = KEND ! reset memory for each batch 1290 1291 IF (LOCDBG) THEN 1292 WRITE (LUPRI,*) MSGDBG, IBATCH, 1293 & '-th. batch out of ',NBATCH 1294 WRITE (LUPRI,*) MSGDBG, 'IR:',IRHGH(IBATCH-1)+1,' -- ', 1295 & IRHGH(IBATCH) 1296 WRITE (LUPRI,*) MSGDBG, 'I2:',I2HGH(IBATCH-1)+1,' -- ', 1297 & I2HGH(IBATCH) 1298 END IF 1299 1300 IF (NBATCH.GT.1) THEN 1301 1302 DTIME = SECOND() 1303 1304 CALL CCBPRE1(INTMEDR,IRHGH(IBATCH-1)+1,IRHGH(IBATCH), 1305 & KRHO2, KLAMP, KLAMH, KDENS, KFOCK, KRIM, 1306 & LUBF,BFFIL,LENBF,LUFK,FILFCK,LENFK, 1307 & LURIM,FILRIM,LENR, 1308 & WORK(KLAMP0), WORK(KLAMH0), 1309 & WORK, LWORK, KEND2, JEND2 ) 1310 KEND2 = JEND2 1311 1312 CALL CCFPRE2(INTMED2,I2HGH(IBATCH-1)+1,I2HGH(IBATCH), 1313 & KFIM, LUFIM,FIMFIL,LENFIM, 1314 & KGIM, LUGIM,GIMFIL,LENGIM, 1315 & KRHO1,LURHO,RHOFIL,LENRHO, 1316 & KMGDL,LUMGD,FILMGD,LENMGD, 1317 & KZDEN,LUZDN,FILZDN,LENZDN, 1318 & KCTR2,KLAMPA,KLAMHA, 1319 & WORK(KLAMP0),WORK(KLAMH0), 1320 & WORK, LWORK, KEND2, JEND2 ) 1321 KEND2 = JEND2 1322 1323 TIMPRE = TIMPRE + SECOND() - DTIME 1324 1325 END IF 1326 1327 LWRK2 = LWORK - KEND2 1328 IF (LWRK2 .LT. 0) THEN 1329 CALL QUIT('Insufficient work space in CC_FMAT. (2)') 1330 END IF 1331 1332*---------------------------------------------------------------------* 1333* loop over number of distributions on the disk: 1334*---------------------------------------------------------------------* 1335 DO IDEL2 = 1, NUMDIS 1336 1337 IF (DIRECT) THEN 1338 IDEL = INDEXA(IDEL2) 1339 IF (NOAUXB) THEN 1340 IDUM = 1 1341 CALL IJKAUX(IDEL,IDUM,IDUM,IDUM) 1342 END IF 1343 ISYDEL = ISAO(IDEL) 1344 ELSE 1345 IDEL = IBAS(ISYMD1) + ILLL 1346 ISYDEL = ISYMD1 1347 END IF 1348cch 1349c write(lupri,*) 'NOAUXB,IDEL2,IDEL,ISYDEL:', 1350c & NOAUXB,IDEL2,IDEL,ISYDEL 1351cch 1352 1353* read AO integral distribution and calculate integrals with 1354* one index transformed to occupied MO (particle): 1355 1356 KXINT = KEND2 1357 KEND3 = KXINT + NDISAO(ISYDEL) 1358 1359 IF (CC2) THEN 1360 KDSRHF = KEND3 1361 KDSRHFA = KDSRHF + NDSRHF(ISYDEL) 1362 KEND3 = KDSRHFA + MDSRHF 1363 ELSE IF (CCSD .OR. CCSDT) THEN 1364 K3OINT = KEND3 1365 KBFRHF = K3OINT + NMAIJK(ISYDEL) 1366 KDCRHF = KBFRHF + NBSRHF(ISYDEL) 1367 KDSRHF = KDCRHF + NBSRHF(ISYDEL) 1368 KEND3 = KDSRHF + NDSRHF(ISYDEL) 1369 END IF 1370 1371 LWRK3 = LWORK - KEND3 1372 IF (LWRK3 .LT. 0) THEN 1373 CALL QUIT('Insufficient work space in CC_FMAT. (3)') 1374 END IF 1375 1376 DTIME = SECOND() 1377 CALL CCRDAO(WORK(KXINT),IDEL,IDEL2,WORK(KEND3),LWRK3, 1378 & WORK(KRECNR),DIRECT) 1379 TIMRDAO = TIMRDAO + SECOND() - DTIME 1380 1381C 1382 IF (LOCDBG) THEN 1383 XNORM = DDOT(NDISAO(ISYDEL),WORK(KXINT),1,WORK(KXINT),1) 1384 WRITE (LUPRI,*) 'IDEL,XNORM:',IDEL,XNORM 1385 END IF 1386C 1387 1388 IF (CCSD .OR. CCSDT) THEN 1389C ----------------------------------------------------- 1390C some integral transformations and presortings for 1391C BF, C, and D intermediate and 21G term: 1392C DSRHF -- standard (**|kdel) integrals 1393C BFRHF -- (**|kdel) presorted for B & D interm. 1394C DCRHF -- (**|kdel) presorted for C & D interm. 1395C 3OINT -- (ij|kdel) integrals for 21G term 1396C ----------------------------------------------------- 1397 DTIME = SECOND() 1398 1399 CALL CCTRBT(WORK(KXINT),WORK(KDSRHF),WORK(KLAMP0),ISYM0, 1400 & WORK(KEND3),LWRK3,ISYDEL) 1401 1402 CALL CC_INT3O(WORK(K3OINT),WORK(KDSRHF),WORK(KLAMH0), 1403 & ISYM0,WORK(KLAMP0),WORK(KEND3), 1404 & LWRK3,IDEL,ISYDEL,LUO3,FILO3) 1405 1406 CALL CC_BFBSORT(WORK(KDSRHF),WORK(KBFRHF),ISYDEL) 1407 1408 CALL CCB_CDSORT(WORK(KXINT),ISYDEL,WORK(KDCRHF), 1409 & WORK(KLAMP0),ISYM0,WORK(KEND3),LWRK3) 1410 1411 TIMTRBT = TIMTRBT + SECOND() - DTIME 1412 TIMIM0 = TIMIM0 + SECOND() - DTIME 1413 1414 KEND3 = KDSRHF 1415 LWRK3 = LWORK - KEND3 1416 1417 KDSRHF = KDUM 1418 1419 END IF 1420 1421 1422C ---------------------------------------------------------- 1423C calculate zeroth-order Fock matrix (Fhat), as well as the 1424C pure core contribution to the zeroth-order Fock matrix 1425C ---------------------------------------------------------- 1426 IF (IBATCH .EQ. 1) THEN 1427 DTIME = SECOND() 1428 1429 CALL CC_AOFOCK(WORK(KXINT), WORK(KDENS0), 1430 * WORK(KFOCK0),WORK(KEND3), 1431 * LWRK3,IDEL,ISYDEL,.FALSE.,DUMMY,ISYM0) 1432 1433 IF (FROIMP.OR.FROEXP) THEN 1434 CALL CC_AOFOCK(WORK(KXINT), WORK(KDNSC0), 1435 * WORK(KFCKC0),WORK(KEND3), 1436 * LWRK3,IDEL,ISYDEL,.FALSE.,DUMMY,ISYM0) 1437 END IF 1438 1439 TIMFCK = TIMFCK + SECOND() - DTIME 1440 TIMIM0 = TIMIM0 + SECOND() - DTIME 1441 END IF 1442 1443 1444C ---------------------------------------------------------- 1445C calculate intermediates that depend only on right vectors: 1446C ---------------------------------------------------------- 1447 DO IINTR = IRHGH(IBATCH-1)+1, IRHGH(IBATCH) 1448 IDLSTR = INTMEDR(1,IINTR) 1449 ISYMR = ILSTSYM(LISTR,IDLSTR) 1450 1451* calculate addresses for C & D intermediates: 1452 IT2DLR(IDEL,IINTR) = ICDEL2 1453 ICDEL2 = ICDEL2 + NT2BCD(MULD2H(ISYDEL,ISYMR)) 1454 1455 DTIME = SECOND() 1456 CALL CCBINT1(WORK(KXINT), WORK(KBFRHF), WORK(KDCRHF), 1457 & IDEL, ISYDEL, WORK(KRHO2(IINTR)), 1458 & WORK(KLAMP0), WORK(KLAMH0), 1459 & WORK(KLAMP(IINTR)), WORK(KLAMH(IINTR)), 1460 & ISYMR, IINTR, 1461 & WORK(KDENS(IINTR)), WORK(KFOCK(IINTR)), 1462 & WORK(KRIM(IINTR)), 1463 & LUC, CTFIL, LUD, DTFIL, 1464 & LUBFD, FNBFD, IADRBFD(1,IINTR), 1465 & WORK(KEND3), LWRK3, 1466 & TIMFCK, TIMBF, TIMC, TIMD ) 1467 TIMIMR = TIMIMR + SECOND() - DTIME 1468 1469 END DO 1470 1471 1472C ---------------------------------------------------------- 1473C calculate intermediates that depend on both the left and 1474C the right vectors: 1475C ---------------------------------------------------------- 1476 DO IINT2 = I2HGH(IBATCH-1)+1, I2HGH(IBATCH) 1477 IDLSTL = INTMED2(1,IINT2) 1478 IDLSTR = INTMED2(3,IINT2) 1479 ISYML = ILSTSYM(LISTL,IDLSTL) 1480 ISYMR = ILSTSYM(LISTR,IDLSTR) 1481 ISYRES = MULD2H(ISYML,ISYMR) 1482 1483 IINTL = ICCSET2(INTMEDL,LISTL,IDLSTL, 1484 & 'R0 ',0,NINTL,MAXSIM,NOAPPEND) 1485 1486 DTIME = SECOND() 1487 CALL CCFINT2(IDEL,ISYDEL,ISYML,WORK(KXINT), 1488 & WORK(KBFRHF),WORK(K3OINT), 1489 & WORK(KMGDL(IINT2)),WORK(KZDEN(IINT2)), 1490 & WORK(KFIM(IINT2)), WORK(KGIM(IINT2)), 1491 & WORK(KRHO1(IINT2)), 1492 & LUBFI,FNBFI,IADRBFI(1,IINT2),ISTARTBFI, 1493 & LUPQR0,FNPQR0,IADRPQ0(1,IINTL), 1494 & LUPQRR,FNPQRR,IADRPQR(1,IINT2), 1495 & LUBDZ0,FNBDZ0,IADRZ0(1,IINTL), 1496 & WORK(KLAMPA(IINT2)),WORK(KLAMHA(IINT2)),ISYMR, 1497 & WORK(KLAMP0), WORK(KLAMH0), 1498 & WORK(KEND3), LWRK3, 1499 & TIMFCK, TIMBF, TIMGZ, TIMIO, TIMFG ) 1500 TIMIM2 = TIMIM2 + SECOND() - DTIME 1501 1502 END DO 1503 1504C ---------------------------------------------------------- 1505C for CC2 calculate the 21CD term contributions. 1506C ---------------------------------------------------------- 1507 IF (CC2) THEN 1508C ........................................................ 1509C DSRHF -- (**|kdel) integrals, transformed with XLAMH0: 1510C ........................................................ 1511 DTIME = SECOND() 1512 CALL CCTRBT(WORK(KXINT),WORK(KDSRHF),WORK(KLAMH0),ISYM0, 1513 & WORK(KEND3),LWRK3,ISYDEL) 1514 TIMTRBT = TIMTRBT + SECOND() - DTIME 1515 TIMIM0 = TIMIM0 + SECOND() - DTIME 1516 1517 DO IINT2 = I2HGH(IBATCH-1)+1, I2HGH(IBATCH) 1518 IDLSTL = INTMED2(1,IINT2) 1519 IDLSTR = INTMED2(3,IINT2) 1520 ISYML = ILSTSYM(LISTL,IDLSTL) 1521 ISYMR = ILSTSYM(LISTR,IDLSTR) 1522 ISYRES = MULD2H(ISYML,ISYMR) 1523 1524 IINTL = ICCSET2(INTMEDL,LISTL,IDLSTL, 1525 & 'R0 ',0,NINTL,MAXSIM,NOAPPEND) 1526 1527C ............................................. 1528C calculate the first term depending on DSRHF0: 1529C ............................................. 1530 LLSAME = .FALSE. 1531 FACT = ONE 1532 DTIME = SECOND() 1533 CALL CCG_21CD( WORK(KRHO1(IINT2)),ISYRES, 1534 & WORK(KCTR2(IINT2)),ISYML, 1535 & WORK(KLAMH0), WORK(KLAMP0), 1536 & WORK(KLAMHA(IINT2)),WORK(KLAMPA(IINT2)),ISYMR, 1537 & WORK(KLAMH0), WORK(KLAMP0), ISYM0, 1538 & WORK(KDSRHF), ISYDEL, IDEL, ISYDEL, 1539 & LLSAME, FACT, WORK(KEND3), LWRK3 ) 1540 1541C .............................................. 1542C calculate the second term depending on DSRHFA: 1543C .............................................. 1544 DTIME = SECOND() 1545 CALL CCTRBT(WORK(KXINT),WORK(KDSRHFA), 1546 & WORK(KLAMHA(IINT2)),ISYMR, 1547 & WORK(KEND3),LWRK3,ISYDEL) 1548 TIMTRBT = TIMTRBT + SECOND() - DTIME 1549 1550 1551 LLSAME = .TRUE. 1552 FACT = HALF 1553 ISYMXA = MULD2H(ISYDEL,ISYMR) 1554 DTIME = SECOND() 1555 CALL CCG_21CD( WORK(KRHO1(IINT2)),ISYRES, 1556 & WORK(KCTR2(IINT2)),ISYML, 1557 & WORK(KLAMH0), WORK(KLAMP0), 1558 & WORK(KLAMH0), WORK(KLAMP0), ISYM0, 1559 & WORK(KLAMH0), WORK(KLAMP0), ISYM0, 1560 & WORK(KDSRHFA), ISYMXA, IDEL, ISYDEL, 1561 & LLSAME, FACT, WORK(KEND3), LWRK3 ) 1562 TIMIM2 = TIMIM2 + SECOND() - DTIME 1563 1564 END DO 1565 END IF 1566 1567 END DO ! IDEL2 1568*---------------------------------------------------------------------* 1569* end of the loop over integral distributions: 1570* if batched I/O algorithm used, save result on disc: 1571*---------------------------------------------------------------------* 1572 IF (NBATCH.GT.1) THEN 1573 DTIME = SECOND() 1574 CALL CCFSAVE(IBATCH, IRHGH, I2HGH, INTMEDR, INTMED2, 1575 & KRHO2, LUBF, BFFIL, LENBF, 1576 & KFOCK, LUFK, FILFCK,LENFK, 1577 & KRIM, LURIM, FILRIM,LENR, 1578 & KFIM, LUFIM, FIMFIL,LENFIM, 1579 & KGIM, LUGIM, GIMFIL,LENGIM, 1580 & KRHO1, LURHO, RHOFIL,LENRHO, 1581 & NINTR, NINT2, WORK, LWORK ) 1582 TIMIO = TIMIO + SECOND() - DTIME 1583 END IF 1584 1585 1586 END DO ! IBATCH 1587 END DO ! ILLL 1588 END DO ! ISYMD1 1589*=====================================================================* 1590* End of Loop over AO-integrals 1591*=====================================================================* 1592 1593*---------------------------------------------------------------------* 1594* if in-core algorithm used, save results now on disc: 1595*---------------------------------------------------------------------* 1596 IF (NBATCH.EQ.1) THEN 1597 DTIME = SECOND() 1598 CALL CCFSAVE(1, IRHGH, I2HGH, INTMEDR, INTMED2, 1599 & KRHO2, LUBF, BFFIL, LENBF, 1600 & KFOCK, LUFK, FILFCK,LENFK, 1601 & KRIM, LURIM, FILRIM,LENR, 1602 & KFIM, LUFIM, FIMFIL,LENFIM, 1603 & KGIM, LUGIM, GIMFIL,LENGIM, 1604 & KRHO1, LURHO, RHOFIL,LENRHO, 1605 & NINTR, NINT2, WORK, LWORK ) 1606 TIMIO = TIMIO + SECOND() - DTIME 1607 END IF 1608 1609*---------------------------------------------------------------------* 1610* close & delete some files which are no longer needed: 1611*---------------------------------------------------------------------* 1612 DTIME = SECOND() 1613 1614 IF (CCS.OR.CC2) THEN 1615 CALL WCLOSE2(LUZDN, FILZDN,'DELETE') 1616 ELSE IF (.NOT.(CCS.OR.CC2)) THEN 1617 CALL WCLOSE2(LUMGD, FILMGD,'DELETE') 1618 CALL WCLOSE2(LUPQR0,FNPQR0,'DELETE') 1619 CALL WCLOSE2(LUPQRR,FNPQRR,'DELETE') 1620 CALL WCLOSE2(LUBDZ0,FNBDZ0,'DELETE') 1621 CALL WCLOSE2(LUO3 , FILO3, 'DELETE') 1622 END IF 1623 1624 TIMPRE = TIMPRE + SECOND() - DTIME 1625 1626*---------------------------------------------------------------------* 1627* recover workspace: 1628*---------------------------------------------------------------------* 1629 KEND1 = KEND0 1630 LWRK1 = LWRK0 1631 1632 1633 IF (LOCDBG) THEN 1634 WRITE (LUPRI,*) MSGDBG,'Loop over AO-integrals completed ', 1635 & ' & AO intermediates saved on file.' 1636 WRITE (LUPRI,*) MSGDBG,'recover work space: KEND1,LWRK1=', 1637 & KEND1,LWRK1 1638 WRITE (LUPRI,*) MSGDBG,'norm of XLAMH0:', 1639 & DDOT(NLAMDT,WORK(KLAMH0),1,WORK(KLAMH0),1) 1640 CALL FLSHFO(LUPRI) 1641 END IF 1642 1643*---------------------------------------------------------------------* 1644* correct response BF intermediate for CCSD(R12) and higher * 1645*---------------------------------------------------------------------* 1646 IF (CCSDR12) THEN 1647 DO IINTR = 1, NINTR 1648 IDLSTR = INTMEDR(1,IINTR) 1649 ISYM = ILSTSYM(LISTR,IDLSTR) 1650 !allocate scratch memory 1651 KSCR1 = KEND1 1652 KEND2 = KSCR1 + NT2AO(ISYM) 1653C KEND2 = KSCR1 + 2*NT2ORT(ISYAMP) 1654 LWRK2 = LWORK - KEND2 1655 IF (LWRK2 .LE. 0) THEN 1656 CALL QUIT('Insufficient work space in CC_FMAT. (7)') 1657 END IF 1658 1659 CALL DZERO(WORK(KSCR1),NT2AO(ISYM)) 1660C CALL DZERO(WORK(KSCR1),2*NT2ORT(ISYAMP)) 1661 1662 !calculate contribution: 1663 IOPT = 1 1664C IOPT = 0 1665 IAMP = 2 1666 CALL CCRHS_BP(WORK(KSCR1),ISYM,IOPT,IAMP,DUMMY,IDUMMY, 1667 & IDUMMY,LISTR,IDLSTR,KETSCL,WORK(KEND2),LWRK2) 1668 1669 !read in response BF intermediate 1670 KSCR2 = KEND2 1671 KEND2 = KSCR2 + NT2AOIJ(ISYM) 1672 LWRK2 = LWORK - KEND2 1673 IF (LWRK2 .LE. 0) THEN 1674 CALL QUIT('Insufficient work space in CC_FMAT. (7)') 1675 END IF 1676 DTIME = SECOND() 1677 CALL CC_RVEC(LUBF,BFFIL,LENBF,NT2AOIJ(ISYM),IINTR, 1678 & WORK(KSCR2)) 1679 TIMIO = TIMIO + SECOND() - DTIME 1680 1681 !transform beta index to occupied and add on BF intermediate: 1682 OSQSAV = OMEGSQ 1683 OORSAV = OMEGOR 1684 OMEGSQ = .FALSE. 1685 OMEGOR = .FALSE. 1686C OMEGOR = .TRUE. 1687 ICON = 4 1688 CALL CC_T2MO(DUMMY,DUMMY,1,WORK(KSCR1),DUMMY,WORK(KSCR2), 1689 & WORK(KLAMP0),WORK(KLAMP0),1,WORK(KEND2),LWRK2, 1690 & ISYM,ICON) 1691 OMEGSQ = OSQSAV 1692 OMEGOR = OORSAV 1693 1694 !write out response BF intermediate 1695 DTIME = SECOND() 1696 CALL CC_WVEC(LUBF,BFFIL,LENBF,NT2AOIJ(ISYM),IINTR, 1697 & WORK(KSCR2)) 1698 TIMIO = TIMIO + SECOND() - DTIME 1699 1700 END DO 1701 END IF 1702 1703*---------------------------------------------------------------------* 1704* open files for CBAR, DBAR, X & Y intermediates: 1705*---------------------------------------------------------------------* 1706 LUCBAR = -1 1707 LUDBAR = -1 1708 LUPQMO = -1 1709 LUXIM = -1 1710 LUYIM = -1 1711 LUMIM = -1 1712 IF (.NOT.(CCS.OR.CC2)) THEN 1713 CALL WOPEN2(LUCBAR,CBAFIL, 64,0) 1714 CALL WOPEN2(LUDBAR,DBAFIL, 64,0) 1715 CALL WOPEN2(LUPQMO,FNPQMO, 64,0) 1716 END IF 1717 1718 IF (.NOT.CCS) THEN 1719 CALL WOPEN2(LUXIM, FILXIM, 64,0) 1720 CALL WOPEN2(LUYIM, FILYIM, 64,0) 1721 END IF 1722 1723 IF (CCSDR12) CALL WOPEN2(LUMIM, FILMIM, 64,0) 1724 1725*---------------------------------------------------------------------* 1726* save a copy of the zeroth-order AO Fock matrix: 1727*---------------------------------------------------------------------* 1728 DTIME = SECOND() 1729 1730 CALL DCOPY(N2BAST,WORK(KFOCK0),1,WORK(KFOCK0AO),1) 1731 1732 TIMFCK = TIMFCK + SECOND() - DTIME 1733 TIMIM0 = TIMIM0 + SECOND() - DTIME 1734 1735*---------------------------------------------------------------------* 1736* do some printout: 1737*---------------------------------------------------------------------* 1738 IF (IPRINT.GT.2) THEN 1739 1740 WRITE (LUPRI,'(1X,A,F10.2,A)') 1741 & '| time for zero order intermediat.:',TIMIM0 ,' secs.|' 1742 WRITE (LUPRI,'(1X,A,I3,A,F10.2,A)') 1743 & '| time for',NINTL,' sets of left interm.:',TIMIML ,' secs.|' 1744 WRITE (LUPRI,'(1X,A,I3,A,F10.2,A)') 1745 & '| time for',NINTR,' sets of right inter.:',TIMIMR ,' secs.|' 1746 WRITE (LUPRI,'(1X,A,I3,A,F10.2,A)') 1747 & '| time for',NINT2,' sets of 2. ord. int.:',TIMIM2 ,' secs.|' 1748 WRITE (LUPRI,'(1X,A,I3,A)') 1749 & '| intermediates calculated in ',NBATCH,' batches |' 1750 1751 WRITE (LUPRI,'(1X,A1,50("-"),A1)') '+','+' 1752 1753 IF (IOPTRES.EQ.5) THEN 1754 WRITE (LUPRI,'(1X,A)') 1755 & '| L vector | R vector | # products | |' 1756 WRITE (LUPRI,'(1X,3(A,A3),A)') '| ',LISTL(1:3), ' No. | ', 1757 & LISTR(1:3),' No. | with ', FILFMA(1:3), 1758 & ' | time/secs |' 1759 ELSE 1760 WRITE (LUPRI,'(1X,A)') 1761 & '| L vector | R vector | result | |' 1762 WRITE (LUPRI,'(1X,A2,A3,2(A,A3),A)') '| ',LISTL(1:3), 1763 & ' No. | ', LISTR(1:3),' No. | ', FILFMA(1:3), 1764 & ' No. | time/secs |' 1765 END IF 1766 1767 WRITE (LUPRI,'(1X,A1,50("-"),A1)') '+','+' 1768 END IF 1769 1770*=====================================================================* 1771* calculate F matrix transformations: 1772*=====================================================================* 1773 IADRTH = 1 1774 DO ITRAN = 1, NFTRAN 1775 1776 IDLSTL = IFTRAN(1,ITRAN) 1777 IDLSTR = IFTRAN(2,ITRAN) 1778 1779 ISYCTR = ILSTSYM(LISTL,IDLSTL) 1780 ISYAMP = ILSTSYM(LISTR,IDLSTR) 1781 ISYRES = MULD2H(ISYAMP,ISYCTR) 1782 1783 IINTL = ICCSET2(INTMEDL,LISTL,IDLSTL, 1784 & 'R0 ',0, NINTL,MAXSIM,NOAPPEND) 1785 IINTR = ICCSET1(INTMEDR,LISTR,IDLSTR,NINTR,MAXSIM,NOAPPEND) 1786 IINT2 = ICCSET2(INTMED2,LISTL,IDLSTL, 1787 & LISTR,IDLSTR,NINT2,MAXSIM,NOAPPEND) 1788 1789 TIMTRN = SECOND() 1790 1791 IF (LOCDBG) THEN 1792 WRITE (LUPRI,*) MSGDBG,'F matrix transformation for ITRAN,', 1793 & ITRAN 1794 WRITE (LUPRI,*) MSGDBG,'IADRTH:',IADRTH 1795 WRITE (LUPRI,*) MSGDBG,'LISTL,IDLSTL:',LISTL(1:3),IDLSTL 1796 WRITE (LUPRI,*) MSGDBG,'LISTR,IDLSTR:',LISTR(1:3),IDLSTR 1797 WRITE (LUPRI,*) MSGDBG,'ISYCTR,ISYAMP,ISYRES:',ISYCTR,ISYAMP, 1798 & ISYRES 1799 WRITE (LUPRI,*) MSGDBG,'IINTL,IINTR,IINT2:',IINTL,IINTR,IINT2 1800 END IF 1801 1802*---------------------------------------------------------------------* 1803* read the single excitation part of the response vectors and 1804* lagrangian multipliers and calculate the response Lambda matrices: 1805*---------------------------------------------------------------------* 1806 KTHETA1 = KEND1 1807 KEND2 = KTHETA1 + NT1AM(ISYRES) 1808 KTHETA2 = KDUM 1809 1810 CALL DZERO(WORK(KTHETA1),NT1AM(ISYRES)) 1811 1812 DTIME = SECOND() 1813 1814 KZETA1 = KEND2 1815 KT1AMPA = KZETA1 + NT1AM(ISYCTR) 1816 KLAMDPA = KT1AMPA + NT1AM(ISYAMP) 1817 KLAMDHA = KLAMDPA + NGLMDT(ISYAMP) 1818 KEND2 = KLAMDHA + NGLMDT(ISYAMP) 1819 LWRK2 = LWORK - KEND2 1820 1821 IF (LWRK2 .LE. 0) THEN 1822 CALL QUIT('Insufficient work space in CC_FMAT. (8)') 1823 END IF 1824 1825 1826 IOPT = 1 1827 1828 CALL CC_RDRSP(LISTL,IDLSTL,ISYCTR,IOPT,MODEL, 1829 & WORK(KZETA1),WORK(KDUM)) 1830 1831 CALL CC_RDRSP(LISTR,IDLSTR,ISYAMP,IOPT,MODEL, 1832 & WORK(KT1AMPA),WORK(KDUM)) 1833 1834 CALL CCLR_LAMTRA(WORK(KLAMP0),WORK(KLAMDPA), WORK(KLAMH0), 1835 & WORK(KLAMDHA),WORK(KT1AMPA),ISYAMP) 1836 1837 TIMPRE = TIMPRE + SECOND() - DTIME 1838 1839*---------------------------------------------------------------------* 1840* do some memory allocation and read (ia|jb) integrals: 1841*---------------------------------------------------------------------* 1842 IF (.NOT.CCS) THEN 1843 KXBARA = KEND2 1844 KYBARA = KXBARA + NMATIJ(ISYAMP) 1845 KA2IM = KYBARA + NMATAB(ISYAMP) 1846 KEND2 = KA2IM + NT1AM(ISYRES) 1847 1848 KLIAJB = KEND2 1849 KXIAJB = KLIAJB + NT2SQ(ISYM0) 1850 KEND3 = KXIAJB + NT2AM(ISYM0) 1851 ELSE IF (LISTL(1:2).EQ.'L0') THEN 1852 KXIAJB = KEND2 1853 KEND3 = KXIAJB + NT2AM(ISYM0) 1854 ELSE 1855 KEND3 = KEND2 1856 END IF 1857 1858 LWRK3 = LWORK - KEND3 1859 1860 IF (LWRK3 .LE. 0) THEN 1861 CALL QUIT('Insufficient work space in CC_FMAT.') 1862 END IF 1863 1864C ----------------------- 1865C read (ia|jb) integrals: 1866C ----------------------- 1867 IF ((.NOT.CCS) .OR. LISTL(1:2).EQ.'L0') THEN 1868 CALL CCG_RDIAJB(WORK(KXIAJB),NT2AM(ISYM0)) 1869 END IF 1870 1871C -------------------------------------------------------- 1872C for CCSD resort (ia|jb) integrals to full square matrix: 1873C -------------------------------------------------------- 1874 IF (.NOT.CCS) THEN 1875 CALL CC_T2SQ(WORK(KXIAJB),WORK(KLIAJB),ISYM0) 1876 END IF 1877 1878*---------------------------------------------------------------------* 1879* calculate the projection on <HF|: 1880*---------------------------------------------------------------------* 1881 IF (LISTL(1:2).EQ.'L0') THEN 1882 1883C ----------------------------------------------------------- 1884C calculate packed L(ia,jb) integrals and evaluate the 1885C projection contribution <HF|[[H,T],t_mu1]|CC> 1886C ----------------------------------------------------------- 1887 IOPTTCME = 1 1888 CALL CCSD_TCMEPK(WORK(KXIAJB),ONE,ISYM0,IOPTTCME) 1889 1890 IOPT = 0 1891 CALL CCG_LXD(WORK(KTHETA1),ISYRES,WORK(KT1AMPA),ISYAMP, 1892 & WORK(KXIAJB),ISYM0,IOPT) 1893 CALL DSCAL(NT1AM(ISYRES),TWO,WORK(KTHETA1),1) 1894 1895 IF (LOCDBG) THEN 1896 WRITE (LUPRI,*) MSGDBG, 1897 & 'norm of THETA1 after <HF| contrib.:', 1898 & DDOT(NT1AM(ISYRES),WORK(KTHETA1),1,WORK(KTHETA1),1) 1899C WRITE (LUPRI,*) 1900C & MSGDBG,'result vector after <HF| contribution:' 1901C CALL CC_PRP(WORK(KTHETA1),WORK(KTHETA2),ISYRES,1,0) 1902 END IF 1903 END IF 1904 1905*---------------------------------------------------------------------* 1906* for CCSD calculate second 21G contribution 1907*---------------------------------------------------------------------* 1908 IF (CCSD .OR. CCSDT) THEN 1909 1910C ---------------------------------------------------------- 1911C memory allocation, reusing space for packed (ia|jb) array: 1912C ---------------------------------------------------------- 1913 KXIDJL = KLIAJB + NT2SQ(ISYM0) 1914 KXJLID = KXIDJL + NTRAOC(ISYAMP) 1915 KEND3 = KXJLID + NTRAOC(ISYAMP) 1916 LWRK3 = LWORK - KEND3 1917 1918 IF (LWRK3 .LE. 0) THEN 1919 CALL QUIT('Insufficient work space in CC_FMAT. (21G)') 1920 END IF 1921 1922C ---------------------------------------------------------- 1923C calculate g_iljd = (ia|jd) * T_al and compute 21G contrib. 1924C ---------------------------------------------------------- 1925 CALL CCG_TRANS4(WORK(KXIDJL),ISYAMP,WORK(KLIAJB),ISYM0, 1926 & WORK(KT1AMPA),ISYAMP) 1927 1928 IOPT = 5 1929 CALL CCG_SORT1(WORK(KXIDJL),WORK(KXJLID),ISYAMP,IOPT) 1930 1931 CALL CC_21GMO(WORK(KTHETA1),ISYRES,WORK(KXJLID),ISYAMP, 1932 & LUPQMO,FNPQMO,IADRPQMO(1,IINTL),ISYCTR, 1933 & WORK(KEND3),LWRK3) 1934 1935 IF (LOCDBG) THEN 1936 WRITE (LUPRI,*) 'Theta1 after 21GMO contribution:' 1937 CALL CC_PRP(WORK(KTHETA1),WORK(KTHETA2),ISYRES,1,0) 1938 END IF 1939 END IF 1940 1941*---------------------------------------------------------------------* 1942* for CC2 & CCSD precalculate the YBARA and A2 intermediates 1943* for CC2 precalculate also the XBARA intermediate 1944*---------------------------------------------------------------------* 1945 IF (.NOT.CCS) THEN 1946 1947 KT2AMPA = KLIAJB + NT2SQ(ISYM0) 1948 KEND3 = KT2AMPA + MAX(NT2AM(ISYAMP),NT2AM(ISYM0)) 1949 LWRK3 = LWORK - KEND3 1950 1951 IF (LWRK3 .LE. 0) THEN 1952 CALL QUIT('Insufficient work space in CC_FMAT. '// 1953 & '(YBARA/A2IM)') 1954 END IF 1955 1956C ---------------------------------------------- 1957C calculate Liajb and read response amplitudes: 1958C ---------------------------------------------- 1959 CALL CCRHS_T2TR(WORK(KLIAJB),WORK(KEND3),LWRK3,ISYM0) 1960 1961 IOPT = 2 1962 CALL CC_RDRSP(LISTR,IDLSTR,ISYAMP,IOPT,MODEL, 1963 & WORK(KDUM),WORK(KT2AMPA)) 1964 CALL CCLR_DIASCL(WORK(KT2AMPA),TWO,ISYAMP) 1965 1966C ---------------------------------------------- 1967C XBARA_kj = sum_dem T_djem L_mekd 1968C ---------------------------------------------- 1969 IF (CC2) THEN 1970 CALL CC_XBAR(WORK(KXBARA),WORK(KLIAJB),ISYM0, 1971 & WORK(KT2AMPA),ISYAMP,WORK(KEND3),LWRK3) 1972 END IF 1973 1974C ---------------------------------------------- 1975C YBARA_ba = sum_dlk T_dlbk L_ldka 1976C ---------------------------------------------- 1977 CALL CC_YBAR(WORK(KYBARA),WORK(KLIAJB),ISYM0, 1978 & WORK(KT2AMPA),ISYAMP,WORK(KEND3),LWRK3) 1979 1980C ---------------------------------------------- 1981C A2IM = sum_bj (2T_aibj - T_ajbi) Zeta_bj 1982C ---------------------------------------------- 1983 IOPTTCME = 1 1984 CALL CCSD_TCMEPK(WORK(KT2AMPA),ONE,ISYAMP,IOPTTCME) 1985 CALL CCG_LXD(WORK(KA2IM),ISYRES,WORK(KZETA1),ISYCTR, 1986 & WORK(KT2AMPA),ISYAMP,0) 1987 1988 END IF 1989 1990*---------------------------------------------------------------------* 1991* allocate & initialize doubles result vector: 1992*---------------------------------------------------------------------* 1993 IF (.NOT.CCS) THEN 1994 KTHETA2 = KEND2 1995 KEND2 = KTHETA2 + NT2AM(ISYRES) 1996 CALL DZERO(WORK(KTHETA2),NT2AM(ISYRES)) 1997 END IF 1998 1999*---------------------------------------------------------------------* 2000* allocate work space some `small' intermediates: 2001*---------------------------------------------------------------------* 2002 IF (CCS.OR.CC2) THEN 2003 KFOCKA = KEND2 2004 KGZETA = KFOCKA + N2BST(ISYAMP) 2005 KEND2 = KGZETA + N2BST(ISYRES) 2006 END IF 2007 2008 IF (.NOT.CCS) THEN 2009 KRIMA = KEND2 2010 KXINT = KRIMA + NEMAT1(ISYAMP) 2011 KYINT = KXINT + NMATIJ(ISYRES) 2012 KFINT = KYINT + NMATAB(ISYRES) 2013 KEND2 = KFINT + NT1AO(ISYRES) 2014 END IF 2015 2016 LWRK2 = LWORK - KEND2 2017 IF (LWRK2 .LE. NT1AM(ISYRES)) THEN 2018 CALL QUIT('Insufficient work space in CC_FMAT. (8)') 2019 END IF 2020 2021*---------------------------------------------------------------------* 2022* for CCS/CC2 read the response Fock matrix and the Fock matrix like 2023* CCS/CC2 GZeta intermediate and calculate the CC_11A contribution 2024*---------------------------------------------------------------------* 2025 IF (CCS.OR.CC2) THEN 2026 IF (LOCDBG) THEN 2027 WRITE (LUPRI,*) 'Norm^2 of THETA1 before GZETA contrib.: ', 2028 & DDOT(NT1AM(ISYRES),WORK(KTHETA1),1,WORK(KTHETA1),1) 2029 END IF 2030 2031 CALL CC_RVEC(LUFK,FILFCK,LENFK,N2BST(ISYAMP),IINTR, 2032 & WORK(KFOCKA)) 2033 CALL CC_FCKMO(WORK(KFOCKA),WORK(KLAMP0),WORK(KLAMH0), 2034 & WORK(KEND2),LWRK2,ISYAMP,ISYM0,ISYM0) 2035 2036 CALL CC_RVEC(LUGIM,GIMFIL,LENGIM,N2BST(ISYRES),IINT2, 2037 & WORK(KGZETA)) 2038 CALL CC_11A(WORK(KTHETA1),WORK(KGZETA),ISYRES,WORK(KLAMH0), 2039 & WORK(KLAMP0),WORK(KEND2),LWRK2) 2040 2041 IF (LOCDBG) THEN 2042 WRITE (LUPRI,*) 'CCS GZETA intermediate:' 2043 CALL CC_PRFCKAO(WORK(KGZETA),ISYRES) 2044 WRITE (LUPRI,*) 'Norm^2 of THETA1 after GZETA contrib.: ', 2045 & DDOT(NT1AM(ISYRES),WORK(KTHETA1),1,WORK(KTHETA1),1) 2046 END IF 2047 END IF 2048 2049*---------------------------------------------------------------------* 2050* for CC2 and CCSD read precalculated X and Y intermediates from file, 2051* and add the precalculated contributions from Rho1 to result vector: 2052*---------------------------------------------------------------------* 2053 IF (.NOT.CCS) THEN 2054 CALL CC_RVEC(LUXIM,FILXIM,LENX,NMATIJ(ISYRES),IINT2, 2055 & WORK(KXINT)) 2056 2057 CALL CC_RVEC(LUYIM,FILYIM,LENY,NMATAB(ISYRES),IINT2, 2058 & WORK(KYINT)) 2059 2060 CALL CC_RVEC(LURHO,RHOFIL,LENRHO,NT1AM(ISYRES),IINT2, 2061 & WORK(KEND2)) 2062 2063 CALL DAXPY(NT1AM(ISYRES),ONE,WORK(KEND2),1,WORK(KTHETA1),1) 2064 2065 IF (LOCDBG) THEN 2066 WRITE (LUPRI,*) 'restored X and Y intermediates:' 2067 CALL CC_PREI(WORK(KYINT),WORK(KXINT),ISYRES,1) 2068 WRITE (LUPRI,*) LUXIM,FILXIM,LENX,NMATIJ(ISYRES),IINT2 2069 WRITE (LUPRI,*) LUYIM,FILYIM,LENY,NMATAB(ISYRES),IINT2 2070 2071 WRITE (LUPRI,*) 'Rho1 read from file:' 2072 CALL CC_PRP(WORK(KEND2),WORK(KEND2),ISYRES,1,0) 2073 WRITE (LUPRI,*) 'result vector after Rho1 was added:' 2074 CALL CC_PRP(WORK(KTHETA1),WORK(KTHETA2),ISYRES,1,0) 2075 2076 CALL FLSHFO(LUPRI) 2077 END IF 2078 2079 END IF 2080 2081*---------------------------------------------------------------------* 2082* for CCSD read also the F and R intermediates; 2083* the F intermediate is transformed to MO and added to result vector 2084*---------------------------------------------------------------------* 2085 IF (.NOT.(CCS.OR.CC2)) THEN 2086 CALL CC_RVEC(LUFIM,FIMFIL,LENFIM,NT1AO(ISYRES),IINT2, 2087 & WORK(KFINT)) 2088 CALL CC_T1AM(WORK(KTHETA1),ISYRES,WORK(KFINT),ISYRES, 2089 & WORK(KLAMH0),ISYM0,ONE) 2090 2091 CALL CC_RVEC(LURIM,FILRIM,LENR,NEMAT1(ISYAMP),IINTR, 2092 & WORK(KRIMA)) 2093 END IF 2094 2095*---------------------------------------------------------------------* 2096* calculate the remaining contributions: 2097*---------------------------------------------------------------------* 2098 2099 IF (LOCDBG) THEN 2100 WRITE (LUPRI,*) MSGDBG,'norm of THETA1 before CCFMAT2:', 2101 & DDOT(NT1AM(ISYRES),WORK(KTHETA1),1,WORK(KTHETA1),1) 2102 IF (.NOT.CCS) 2103 & WRITE (LUPRI,*) MSGDBG,'norm of THETA2 before CCFMAT2:', 2104 & DDOT(NT2AM(ISYRES),WORK(KTHETA2),1,WORK(KTHETA2),1) 2105 WRITE (LUPRI,*) 'result vector before entering CCFMAT2:' 2106 CALL CC_PRP(WORK(KTHETA1),WORK(KTHETA2),ISYRES,1,0) 2107 END IF 2108 2109 CALL CCFMAT2(WORK(KTHETA1),WORK(KTHETA2),ISYRES, 2110 & LISTL,IDLSTL, WORK(KZETA1), ISYCTR, 2111 & LISTR,IDLSTR, WORK(KT1AMPA),ISYAMP, 2112 & WORK(KLAMDPA),WORK(KLAMDHA), 2113 & WORK(KLAMP0), WORK(KLAMH0), 2114 & WORK(KXINT), WORK(KYINT), 2115 & WORK(KXBARA), WORK(KYBARA), 2116 & WORK(KRIMA), WORK(KA2IM), 2117 & WORK(KFOCKA), WORK(KFOCK0AO), WORK(KFCKC0), 2118 & LUBFI,FNBFI,IADRBFI(1,IINT2), 2119 & LUBF, BFFIL,LENBF,IINTR, 2120 & LUC, CTFIL, LUCBAR, CBAFIL, 2121 & LUD, DTFIL, LUDBAR, DBAFIL, 2122 & IOFFCD(IINTR), WORK(KEND2),LWRK2) 2123 2124*---------------------------------------------------------------------* 2125* allocate memory for R12-part of result vector (stored as symmetry 2126* packed triangular matrix): 2127*---------------------------------------------------------------------* 2128 IF (CCR12) THEN 2129 KTHETAR12 = KEND2 2130 KEND2 = KTHETAR12 + NTR12AM(ISYRES) 2131 LWRK2 = LWORK - KEND2 2132 IF (LWRK2 .LT. 0) THEN 2133 CALL QUIT('Insufficient memory for R12-FMATRIX!') 2134 END IF 2135 CALL DZERO(WORK(KTHETAR12),NTR12AM(ISYRES)) 2136 2137 END IF 2138C 2139*---------------------------------------------------------------------* 2140* calculate R12 contributions: 2141* 2142* C. Neiss, Autumn 2004 2143*---------------------------------------------------------------------* 2144C 2145 IF (CCR12) THEN 2146 IF (LOCDBG) THEN 2147 CALL AROUND('Call R12 section for FMATRIX') 2148 END IF 2149C 2150 CALL CC_R12FMAT(WORK(KTHETA1),WORK(KTHETA2), 2151 & WORK(KTHETAR12), ISYRES, 2152 & LISTL,IDLSTL, WORK(KZETA1), ISYCTR, 2153 & LISTR,IDLSTR, ISYAMP, 2154 & WORK(KLAMDPA),WORK(KLAMDHA), 2155 & WORK(KLAMP0), WORK(KLAMH0), 2156 & WORK(KXINT),LUMIM,FILMIM,LENM,IINT2, 2157 & WORK(KEND2),LWRK2) 2158 IF (LOCDBG) THEN 2159 CALL AROUND('Returned from R12 section for FMATRIX') 2160 END IF 2161C 2162 ENDIF 2163C 2164C------------------------------------------- 2165C SLV98,OC, CCMM JK, NYQMMM KS 10 2166C Calculate T^{gB} solvent contribution. 2167C Unsquared T2 is required! 2168C------------------------------------------- 2169C 2170 IF (CCSLV.AND.LSTBTR) THEN 2171C 2172 KT2AMPA = KEND2 2173 KEND3 = KT2AMPA + NT2AM(ISYAMP) 2174 LWRK3 = LWORK - KEND3 2175 IF (LWRK3 .LT. 0) THEN 2176 CALL QUIT('Insuff. work in CC_FMAT. in CCSLV part') 2177 END IF 2178C 2179C 2180 IF (LISTL(1:2).EQ.'L0') THEN 2181 IOPT = 2 2182 CALL CC_RDRSP(LISTR,IDLSTR,ISYAMP,IOPT,MODEL, 2183 * WORK(KDUM),WORK(KT2AMPA)) 2184 2185 IF (.NOT. CCMM) CALL CCSL_LTRB(WORK(KTHETA1),WORK(KTHETA2), 2186 * WORK(KT1AMPA),WORK(KT2AMPA),ISYAMP,'F', 2187 * WORK(KEND3),LWRK3) 2188 IF (CCMM) THEN 2189 IF (.NOT. NYQMMM) THEN 2190 CALL CCMM_LTRB(WORK(KTHETA1),WORK(KTHETA2), 2191 * WORK(KT1AMPA),WORK(KT2AMPA),ISYAMP,'F', 2192 * WORK(KEND3),LWRK3) 2193 ELSE IF (NYQMMM .AND. (.NOT. HFFLD)) THEN 2194 CALL CCMM_TRANSFORMER(WORK(KTHETA1),WORK(KTHETA2), 2195 * WORK(KT1AMPA),WORK(KT2AMPA),MODEL,ISYAMP, 2196 * 'F',WORK(KEND3),LWRK3) 2197 END IF 2198 END IF 2199 ELSE IF (LISTL(1:2).NE.'L0') THEN 2200 IF (.NOT. CCMM) CALL CCSL_PBTR(WORK(KTHETA1),WORK(KTHETA2), 2201 * ISYRES, 2202 * LISTL,IDLSTL, WORK(KZETA1), 2203 * ISYCTR, 2204 * LISTR,IDLSTR, WORK(KT1AMPA), 2205 * ISYAMP,MODEL, 2206 * WORK(KEND2),LWRK2) 2207 IF (CCMM) THEN 2208 IF (.NOT. NYQMMM) THEN 2209 CALL CCMM_PBTR(WORK(KTHETA1),WORK(KTHETA2), 2210 * ISYRES, 2211 * LISTL,IDLSTL, WORK(KZETA1), 2212 * ISYCTR, 2213 * LISTR,IDLSTR, WORK(KT1AMPA), 2214 * ISYAMP,MODEL, 2215 * WORK(KEND2),LWRK2) 2216 ELSE IF (NYQMMM .AND. (.NOT. HFFLD ) ) THEN 2217 RSPTYP='F' 2218 CALL CCMM_QRTRANSFORMER(WORK(KTHETA1),WORK(KTHETA2), 2219 * ISYRES, 2220 * LISTL,IDLSTL,ISYCTR, 2221 * LISTR,IDLSTR,ISYAMP, 2222 * MODEL,RSPTYP,WORK(KEND2),LWRK2) 2223 END IF 2224 END IF 2225 END IF 2226 ENDIF 2227C 2228 IF (USE_PELIB().AND.LSTBTR) THEN 2229C 2230 KT2AMPA = KEND2 2231 KEND3 = KT2AMPA + NT2AM(ISYAMP) 2232 LWRK3 = LWORK - KEND3 2233 IF (LWRK3 .LT. 0) THEN 2234 CALL QUIT('Insuff. work in CC_FMAT. in PECC part') 2235 END IF 2236C 2237C 2238 IF (LISTL(1:2).EQ.'L0') THEN 2239 IOPT = 2 2240 CALL CC_RDRSP(LISTR,IDLSTR,ISYAMP,IOPT,MODEL, 2241 & WORK(KDUM),WORK(KT2AMPA)) 2242 IF (.NOT. HFFLD) THEN 2243 CALL PELIB_IFC_TRANSFORMER(WORK(KTHETA1),WORK(KTHETA2), 2244 & WORK(KT1AMPA),WORK(KT2AMPA),MODEL,ISYAMP, 2245 & 'F',WORK(KEND3),LWRK3) 2246 END IF 2247C 2248 ELSE IF (LISTL(1:2).NE.'L0') THEN 2249 IF (.NOT. HFFLD) THEN 2250 RSPTYP='F' 2251 CALL PELIB_IFC_QRTRANSFORMER(WORK(KTHETA1),WORK(KTHETA2), 2252 & ISYRES,LISTL,IDLSTL,ISYCTR, 2253 & LISTR,IDLSTR,ISYAMP,MODEL,RSPTYP, 2254 & WORK(KEND2),LWRK2) 2255 END IF 2256 END IF 2257 END IF 2258C---------------------------------------------------------- 2259C 2260 2261 IF (CCSDT) THEN 2262 2263 ! allocate memory for effective rhs vector 2264 IF ( IOPTRES.GE.0 .AND. IOPTRES.LE.4 ) THEN 2265 KTHETA1EFF = KEND2 2266 KTHETA2EFF = KTHETA1EFF + NT1AM(ISYRES) 2267 KEND2 = KTHETA2EFF + NT2AM(ISYRES) 2268 LWRK2 = LWORK - KEND2 2269 IF (LWRK2 .LE. NT1AM(ISYRES)) THEN 2270 CALL QUIT('Insufficient work space in CC_FMAT. (9)') 2271 END IF 2272 END IF 2273 2274 IF ( (.NOT. NODDY_FMAT)) THEN 2275 2276C ---------------------------------------- 2277C Calculate <L_2|[[H,T^B_3],\tau_nu_1|HF> 2278C ---------------------------------------- 2279 IF (IOPTRES.LT.5) THEN 2280 2281 CALL CC3_FT3B(LISTL,IDLSTL,LISTR,IDLSTR,IOPTRES, 2282 * WORK(KTHETA1),ISYRES,WORK(KEND2),LWRK2, 2283 * FNCKJD,FNDKBC,FNDELD,FNTOC) 2284 2285 END IF 2286 2287 LUCKJD = -1 2288 LUDKBC = -1 2289 LUTOC = -1 2290 LU3VI = -1 2291 LU3VI2 = -1 2292 LUDKBC3 = -1 2293 LU3FOP = -1 2294 LU3FOP2 = -1 2295 LU3FOPX = -1 2296 LU3FOP2X = -1 2297 2298 CALL WOPEN2(LUCKJD,FNCKJD,64,0) 2299 CALL WOPEN2(LUDKBC,FNDKBC,64,0) 2300 CALL WOPEN2(LUTOC,FNTOC,64,0) 2301 CALL WOPEN2(LU3VI,FN3VI,64,0) 2302 CALL WOPEN2(LU3VI2,FN3VI2,64,0) 2303 CALL WOPEN2(LUDKBC3,FNDKBC3,64,0) 2304 CALL WOPEN2(LU3FOP,FN3FOP,64,0) 2305 CALL WOPEN2(LU3FOP2,FN3FOP2,64,0) 2306 CALL WOPEN2(LU3FOPX,FN3FOPX,64,0) 2307 CALL WOPEN2(LU3FOP2X,FN3FOP2X,64,0) 2308 2309C 2310C --------------------------------------------- 2311C Calculate <L_3|[[H^B,T^0_2],\tau_nu_1]|HF> 2312C + <L_3|[[H^0,T^B_2],\tau_nu_1]|HF> 2313C and 2314C <L_3|[H^B,\tau_nu_2]|HF> 2315C --------------------------------------------- 2316 IF (IOPTRES.EQ.5) THEN 2317 ! in CC3_BFMAT Omegaeff is used as scratch array. 2318 KTHETA1EFF = KEND2 2319 KTHETA2EFF = KTHETA1EFF + NT1AM(ISYRES) 2320 KEND3 = KTHETA2EFF + NT2AM(ISYRES) 2321 LWRK3 = LWORK - KEND3 2322 2323 IF (LWRK3 .LT. 0) THEN 2324 CALL QUIT('Insufficient work space in CC_FMAT. (9b)') 2325 END IF 2326 2327 ELSE 2328 KEND3 = KEND2 2329 LWRK3 = LWRK2 2330 END IF 2331 2332 CALL DZERO(WORK(KTHETA1EFF),NT1AM(ISYRES)) 2333 CALL DZERO(WORK(KTHETA2EFF),NT2AM(ISYRES)) 2334 2335 CALL CC3_BFMAT(LISTL,IDLSTL,LISTR,IDLSTR, 2336 * WORK(KTHETA1), WORK(KTHETA2), 2337 * WORK(KTHETA1EFF),WORK(KTHETA2EFF), 2338 * ISYRES, 2339 * WORK(KEND3),LWRK3, 2340 * LUTOC,FNTOC, 2341 * LU3VI,FN3VI,LUDKBC3,FNDKBC3, 2342 * LU3FOPX,FN3FOPX,LU3FOP2X,FN3FOP2X, 2343 * LU3VI2,FN3VI2,LU3FOP,FN3FOP,LU3FOP2,FN3FOP2) 2344 2345C 2346C ------------------------------------------ 2347C Calculate <F_3|[[H,T^0_2],\tau_nu_1]|HF> 2348C and <F_3|[H,\tau_nu_2]|HF> 2349C ------------------------------------------ 2350 IF (IOPTRES.LT.5) THEN 2351 2352 CALL DZERO(WORK(KTHETA1EFF),NT1AM(ISYRES)) 2353 CALL DZERO(WORK(KTHETA2EFF),NT2AM(ISYRES)) 2354 2355 CALL CC3_FMATSD(LISTL,IDLSTL,LISTR,IDLSTR, 2356 * WORK(KTHETA1),WORK(KTHETA2),WORK(KTHETA1EFF), 2357 * WORK(KTHETA2EFF),ISYRES,WORK(KEND2),LWRK2, 2358 * LUDKBC,FNDKBC,LUCKJD,FNCKJD,LUTOC,FNTOC, 2359 * LU3VI,FN3VI,LU3VI2,FN3VI2, 2360 * LU3FOP,FN3FOP,LU3FOP2,FN3FOP2) 2361 ENDIF 2362 2363 CALL WCLOSE2(LUCKJD,FNCKJD,'KEEP') 2364 CALL WCLOSE2(LUDKBC,FNDKBC,'KEEP') 2365 CALL WCLOSE2(LUTOC,FNTOC,'KEEP') 2366 CALL WCLOSE2(LU3VI,FN3VI,'KEEP') 2367 CALL WCLOSE2(LU3VI2,FN3VI2,'KEEP') 2368 CALL WCLOSE2(LUDKBC3,FNDKBC3,'KEEP') 2369 CALL WCLOSE2(LU3FOP,FN3FOP,'KEEP') 2370 CALL WCLOSE2(LU3FOP2,FN3FOP2,'KEEP') 2371 CALL WCLOSE2(LU3FOPX,FN3FOPX,'KEEP') 2372 CALL WCLOSE2(LU3FOP2X,FN3FOP2X,'KEEP') 2373 2374 ELSE 2375 2376 IF ( IOPTRES.GE.0 .AND. IOPTRES.LE.4 ) THEN 2377 CALL DZERO(WORK(KTHETA1EFF),NT1AM(ISYRES)) 2378 CALL DZERO(WORK(KTHETA2EFF),NT2AM(ISYRES)) 2379 END IF 2380 2381 CALL CCSDT_FMAT_NODDY(LISTL,IDLSTL,LISTR,IDLSTR, 2382 & IOPTRES,CCSDT_F_ALTER, 2383 & WORK(KTHETA1), WORK(KTHETA2), 2384 & WORK(KTHETA1EFF),WORK(KTHETA2EFF), 2385 & IFDOTS,FCONS,FILFMA,ITRAN, 2386 & NFTRAN,MXVEC,WORK(KEND2),LWRK2) 2387 2388 ENDIF 2389 2390 IF (IPRINT .GT. 110) THEN 2391 CALL AROUND('Omega1 and Omega2 after F(cc3) : ') 2392 CALL CC_PRP(WORK(KTHETA1),WORK(KTHETA2),1,1,1) 2393 ENDIF 2394 2395 END IF 2396 2397 2398 2399*---------------------------------------------------------------------* 2400* write result vector to output: 2401*---------------------------------------------------------------------* 2402 DTIME = SECOND() 2403 2404 KTHETA0 = -999999 2405 2406 IF (IOPTRES .EQ. 0 .OR. IOPTRES .EQ. 1) THEN 2407 2408* write to a common direct access file, 2409* store start address in IFTRAN(3,ITRAN) 2410 2411 IFTRAN(3,ITRAN) = IADRTH 2412 2413 CALL PUTWA2(LUFMAT,FILFMA,WORK(KTHETA1),IADRTH,NT1AM(ISYRES)) 2414 IADRTH = IADRTH + NT1AM(ISYRES) 2415 2416 IF (.NOT.CCS) THEN 2417 CALL PUTWA2(LUFMAT,FILFMA,WORK(KTHETA2),IADRTH,NT2AM(ISYRES)) 2418 IADRTH = IADRTH + NT2AM(ISYRES) 2419 END IF 2420 2421 IF (CCR12) THEN 2422 IF (LOCDBG) THEN 2423 WRITE(LUPRI,*) 'THETAR12 in CC_FMATNEW:' 2424 WRITE(LUPRI,*) 'Will write to file ', FILFMA, LUFMAT 2425 CALL OUTPAK(WORK(KTHETAR12),NMATKI(ISYRES),1,LUPRI) 2426 END IF 2427 CALL PUTWA2(LUFMAT,FILFMA,WORK(KTHETAR12),IADRTH, 2428 & NTR12AM(ISYRES)) 2429 IADRTH = IADRTH + NTR12AM(ISYRES) 2430 END IF 2431 2432C some debug output: 2433 IF (LOCDBG) THEN 2434 WRITE (LUPRI,*) 'F matrix transformation nb. ',ITRAN, 2435 & ' saved on file.' 2436 WRITE (LUPRI,*) 'ADRESS, LENGTH:', 2437 & IFTRAN(3,ITRAN),IADRTH-IFTRAN(3,ITRAN) 2438 XNORM = DDOT(NT1AM(ISYRES),WORK(KTHETA1),1,WORK(KTHETA1),1) 2439 IF (.NOT.CCS) XNORM = XNORM + 2440 & DDOT(NT2AM(ISYRES),WORK(KTHETA2),1,WORK(KTHETA2),1) 2441 IF (CCR12) XNORM = XNORM + 2442 & DDOT(NTR12AM(ISYRES),WORK(KTHETAR12),1, 2443 & WORK(KTHETAR12),1) 2444 WRITE (LUPRI,*) 'Norm:', XNORM 2445 2446 Call AROUND('F matrix transformation written to file:') 2447 NDBLE = 1 2448 IF (CCS) NDBLE = 0 2449 Call CC_PRP(WORK(KTHETA1),WORK(KTHETA2),ISYRES,1,NDBLE) 2450 IF (CCR12) THEN 2451 write(lupri,*) 'THETA_R12 in CC_FMATNEW:' 2452 CALL OUTPAK(WORK(KTHETAR12),NMATKI(ISYRES),1,LUPRI) 2453 END IF 2454 END IF 2455C 2456 2457 ELSE IF ( IOPTRES .EQ. 3 .OR. IOPTRES .EQ. 4 ) THEN 2458 ! write to a sequential file by a call to CC_WRRSP/CC_WARSP, 2459 ! use FILFMA as LIST type and IFTRAN(3,ITRAN) as index 2460 IF (IOPTRES.EQ.3) THEN 2461 CALL CC_WRRSP(FILFMA,IFTRAN(3,ITRAN),ISYRES,IOPTW,MODELW, 2462 & WORK(KTHETA0),WORK(KTHETA1),WORK(KTHETA2), 2463 & WORK(KEND2),LWRK2) 2464 IF (CCR12) THEN 2465 CALL CC_WRRSP(FILFMA,IFTRAN(3,ITRAN),ISYRES,IOPTWR12, 2466 & MODELW,DUMMY,DUMMY,WORK(KTHETAR12), 2467 & WORK(KEND2),LWRK2) 2468 END IF 2469 IF (CCSDT) THEN 2470 CALL CC_WRRSP(FILFMA,IFTRAN(3,ITRAN),ISYRES,IOPTWE,MODELW, 2471 & WORK(KTHETA0),WORK(KTHETA1EFF), 2472 & WORK(KTHETA2EFF),WORK(KEND2),LWRK2) 2473 END IF 2474 ELSE IF (IOPTRES.EQ.4) THEN 2475 CALL CC_WARSP(FILFMA,IFTRAN(3,ITRAN),ISYRES,IOPTW,MODELW, 2476 & WORK(KTHETA0),WORK(KTHETA1),WORK(KTHETA2), 2477 & WORK(KEND2),LWRK2) 2478 IF (CCR12) THEN 2479 CALL CC_WARSP(FILFMA,IFTRAN(3,ITRAN),ISYRES,IOPTWR12, 2480 & MODELW,DUMMY,DUMMY,WORK(KTHETAR12), 2481 & WORK(KEND2),LWRK2) 2482 END IF 2483 IF (CCSDT) THEN 2484 CALL CC_WARSP(FILFMA,IFTRAN(3,ITRAN),ISYRES,IOPTWE,MODELW, 2485 & WORK(KTHETA0),WORK(KTHETA1EFF), 2486 & WORK(KTHETA2EFF),WORK(KEND2),LWRK2) 2487 END IF 2488 END IF 2489 2490 IF (LOCDBG) THEN 2491 WRITE (LUPRI,*) 'Write ',LISTL(1:3),' * B * ',LISTR(1:3), 2492 & ' transformation',' as ',FILFMA(1:3),' type vector to file.' 2493 WRITE (LUPRI,*) 'index of inp. left vector:',IFTRAN(1,ITRAN) 2494 WRITE (LUPRI,*) 'index of inp. right vector:',IFTRAN(2,ITRAN) 2495 WRITE (LUPRI,*) 'index of result vector:',IFTRAN(3,ITRAN) 2496 XNORM = DDOT(NT1AM(ISYRES),WORK(KTHETA1),1,WORK(KTHETA1),1) 2497 IF (.NOT.CCS) XNORM = XNORM + 2498 & DDOT(NT2AM(ISYRES),WORK(KTHETA2),1,WORK(KTHETA2),1) 2499 IF (CCR12) THEN 2500 XNORM = XNORM + DDOT(NTR12AM(ISYRES),WORK(KTHETAR12),1, 2501 & WORK(KTHETAR12),1) 2502 END IF 2503 WRITE (LUPRI,*) 'norm^2 of result vector:',XNORM 2504 CALL CC_PRP(WORK(KTHETA1),WORK(KTHETA2),ISYRES,1,1) 2505 IF (CCR12) CALL CC_PRPR12(WORK(KTHETAR12),ISYRES,1,.TRUE.) 2506 IF (CCSDT) THEN 2507 WRITE(LUPRI,*) 'CCSDT effective vector:' 2508 XNORM = DDOT(NT1AM(ISYRES),WORK(KTHETA1EFF),1, 2509 & WORK(KTHETA1EFF),1) + DDOT(NT2AM(ISYRES), 2510 & WORK(KTHETA2EFF),1,WORK(KTHETA2EFF),1) 2511 WRITE (LUPRI,*) 'norm^2:',XNORM 2512 CALL CC_PRP(WORK(KTHETA1EFF),WORK(KTHETA2EFF),ISYRES,1,1) 2513 END IF 2514 END IF 2515 ELSE IF (IOPTRES.EQ.5) THEN 2516 IF (LOCDBG.AND.CCSDT) THEN 2517 WRITE(LUPRI,*) 'FCONS TRIPLES CONTRIBUTION:' 2518 IVEC = 1 2519 DO WHILE (IFDOTS(IVEC,ITRAN).NE.0 .AND. IVEC.LE.MXVEC) 2520 WRITE(LUPRI,*)'FCONS:',IVEC,ITRAN,FCONS(IVEC,ITRAN),IOPTW 2521 IVEC = IVEC + 1 2522 END DO 2523 END IF 2524 2525 IF (CCR12) THEN 2526 CALL CCDOTRSP(IFDOTS,FCONS,IOPTWR12,FILFMA,ITRAN,NFTRAN, 2527 & MXVEC,DUMMY,WORK(KTHETAR12),ISYRES, 2528 & WORK(KEND2),LWRK2) 2529 END IF 2530 2531 IF (LOCDBG) THEN 2532 write(lupri,*)'FCONS R12 doubles contributions:' 2533 IVEC = 1 2534 DO WHILE (IFDOTS(IVEC,ITRAN).NE.0 .AND. IVEC.LE.MXVEC) 2535 WRITE(LUPRI,*)'FCONS:',IVEC,ITRAN,FCONS(IVEC,ITRAN),IOPTW 2536 IVEC = IVEC + 1 2537 END DO 2538 END IF 2539 2540 IF (.NOT.CCS) CALL CCLR_DIASCL(WORK(KTHETA2),TWO,ISYRES) 2541 CALL CCDOTRSP(IFDOTS,FCONS,IOPTW,FILFMA,ITRAN,NFTRAN,MXVEC, 2542 & WORK(KTHETA1),WORK(KTHETA2),ISYRES, 2543 & WORK(KEND2),LWRK2) 2544 IF (LOCDBG) THEN 2545 IVEC = 1 2546 DO WHILE (IFDOTS(IVEC,ITRAN).NE.0 .AND. IVEC.LE.MXVEC) 2547 WRITE(LUPRI,*)'FCONS:',IVEC,ITRAN,FCONS(IVEC,ITRAN),IOPTW 2548 IVEC = IVEC + 1 2549 END DO 2550 END IF 2551 ELSE 2552 CALL QUIT('Illegal value for IOPTRES in CC_FMAT.') 2553 END IF 2554 2555 TIMIO = TIMIO + SECOND() - DTIME 2556 2557 TIMTRN = SECOND() - TIMTRN 2558 2559 IF (IPRINT.GT.2) THEN 2560 2561 IF (IOPTRES.EQ.5) THEN 2562 IVEC = 1 2563 DO WHILE (IFDOTS(IVEC,ITRAN).NE.0 .AND. IVEC.LE.MXVEC) 2564 IVEC = IVEC + 1 2565 END DO 2566 WRITE (LUPRI,'(1X,2(A,I5),A,I6,A,F10.2,A)')'| ',IDLSTL, 2567 & ' | ',IDLSTR,' | ',IVEC-1,' | ',TIMTRN,' |' 2568 ELSE 2569 WRITE (LUPRI,'(1X,2(A,I5),A,I6,A,F10.2,A)') '| ',IDLSTL, 2570 & ' | ',IDLSTR,' | ',IFTRAN(3,ITRAN),' | ', 2571 & TIMTRN,' |' 2572 END IF 2573 2574 END IF 2575 2576*---------------------------------------------------------------------* 2577* End of loop over F matrix transformations 2578*---------------------------------------------------------------------* 2579 END DO 2580 2581*---------------------------------------------------------------------* 2582* close & remove scratch files: 2583*---------------------------------------------------------------------* 2584 DTIME = SECOND() 2585 2586 IF (.NOT. (CCS.OR.CC2)) THEN 2587 CALL WCLOSE2(LUBF, BFFIL, 'DELETE') 2588 CALL WCLOSE2(LUC, CTFIL, 'DELETE') 2589 CALL WCLOSE2(LUD, DTFIL, 'DELETE') 2590 CALL WCLOSE2(LUCBAR, CBAFIL, 'DELETE') 2591 CALL WCLOSE2(LUDBAR, DBAFIL, 'DELETE') 2592 CALL WCLOSE2(LURIM, FILRIM, 'DELETE') 2593 CALL WCLOSE2(LUFIM, FIMFIL, 'DELETE') 2594 CALL WCLOSE2(LUPQMO, FNPQMO, 'DELETE') 2595 CALL WCLOSE2(LUBFI, FNBFI, 'DELETE') 2596 ELSE 2597 CALL WCLOSE2(LUGIM, GIMFIL, 'DELETE') 2598 CALL WCLOSE2(LUFK, FILFCK, 'DELETE') 2599 END IF 2600 2601 IF (.NOT.CCS) THEN 2602 CALL WCLOSE2(LUXIM,FILXIM, 'DELETE') 2603 CALL WCLOSE2(LUYIM,FILYIM, 'DELETE') 2604 CALL WCLOSE2(LURHO,RHOFIL, 'DELETE') 2605 END IF 2606 IF (CCSD .OR. CCSDT) CALL WCLOSE2(LUBFD, FNBFD, 'DELETE') 2607 2608 IF (CCSDR12) CALL WCLOSE2(LUMIM,FILMIM,'DELETE') 2609 2610 2611 TIMIO = TIMIO + SECOND() - DTIME 2612 2613*---------------------------------------------------------------------* 2614* if IOPTRES=1 and enough work space available, read result 2615* vectors back into memory: 2616*---------------------------------------------------------------------* 2617 DTIME = SECOND() 2618 2619* check size of work space: 2620 IF (IOPTRES .EQ. 1) THEN 2621 LENALL = IADRTH-1 2622 IF (LENALL .GT. LWORK) IOPTRES = 0 2623 END IF 2624 2625* read the result vectors back into memory: 2626 IF (IOPTRES .EQ. 1) THEN 2627 2628 CALL GETWA2(LUFMAT,FILFMA,WORK(1),1,LENALL) 2629 2630 IF (LOCDBG) THEN 2631 DO ITRAN = 1, NFTRAN 2632 IF (ITRAN.LT.NFTRAN) THEN 2633 LEN = IFTRAN(3,ITRAN+1)-IFTRAN(3,ITRAN) 2634 ELSE 2635 LEN = IADRTH-IFTRAN(3,NFTRAN) 2636 END IF 2637 KTHETA1 = IFTRAN(3,ITRAN) 2638 XNORM = DDOT(LEN, WORK(KTHETA1),1, WORK(KTHETA1),1) 2639 WRITE (LUPRI,*) 'Read F matrix transformation nb. ',NFTRAN 2640 WRITE (LUPRI,*) 'Adress, length, NORM:',IFTRAN(3,NFTRAN), 2641 & LEN,XNORM 2642 END DO 2643 CALL FLSHFO(LUPRI) 2644 END IF 2645 END IF 2646 2647 TIMIO = TIMIO + SECOND() - DTIME 2648*---------------------------------------------------------------------* 2649* close F matrix file, print timings & return 2650*---------------------------------------------------------------------* 2651 DTIME = SECOND() 2652 2653 IF (IOPTRES.EQ.0 ) THEN 2654 CALL WCLOSE2(LUFMAT, FILFMA, 'KEEP') 2655 ELSE IF (IOPTRES.EQ.1) THEN 2656 CALL WCLOSE2(LUFMAT, FILFMA, 'DELETE') 2657 ELSE IF (IOPTRES.EQ.3 .OR. IOPTRES.EQ.4 .OR. IOPTRES.EQ.5) THEN 2658 CONTINUE 2659 ELSE 2660 CALL QUIT('Illegal value of IOPTRES in CC_FMAT.') 2661 END IF 2662 2663 TIMIO = TIMIO + SECOND() - DTIME 2664 TIMALL = SECOND() - TIMALL 2665 2666 IF (IPRINT.GT.2) THEN 2667 WRITE (LUPRI,'(1X,A1,50("-"),A1)') '+','+' 2668 WRITE (LUPRI,'(1X,A,I4,A,F10.2,A)') 2669 & '| total time for',NFTRAN,' F transforms.:',TIMALL,' secs.|' 2670 IF (TIMALL.GE.1.0D0) THEN 2671 CONVRT = 100.0D0/TIMALL 2672 WRITE (LUPRI,'(1X,A1,50("-"),A1)')'+','+' 2673 WRITE 2674 & (LUPRI,'(1X,"| % of time used in ",A,":",F10.2," |")') 2675 & 'start up org.', TIMPRE*CONVRT, 2676 & 'Fock interm. ', TIMFCK*CONVRT, 2677 & 'ERI ', TIMINT*CONVRT, 2678 & 'CCRDAO ', TIMRDAO*CONVRT, 2679 & '(**|k del) ', TIMTRBT*CONVRT, 2680 & 'BF term ', TIMBF*CONVRT, 2681 & 'C term ', TIMC*CONVRT, 2682 & 'GZ terms ', TIMGZ*CONVRT, 2683 & 'FG terms ', TIMFG*CONVRT, 2684 & 'D term ', TIMD*CONVRT, 2685 & 'I/O ', TIMIO*CONVRT 2686 END IF 2687 WRITE (LUPRI,'(1X,A1,50("="),A1,//)') '+','+' 2688 END IF 2689*---------------------------------------------------------------------* 2690* that's it; return: 2691*---------------------------------------------------------------------* 2692 2693 CALL QEXIT('CC_FMATNEW') 2694 2695 RETURN 2696 END 2697*=====================================================================* 2698* END OF SUBROUTINE CC_FMAT 2699*=====================================================================* 2700*=====================================================================* 2701C /* Deck ccfpreint */ 2702 SUBROUTINE CCFPREINT(INTMED2,NINT2,XLAMDP,XLAMDH, 2703 & FILXIM,LENX, 2704 & FILYIM,LENY, 2705 & FILMIM,LENM, 2706 & FILMGD,LENMGD, 2707 & FILZDN,LENZDN, 2708 & FILPQ, IADRPQ, 2709 & FILPQMO,IADRPQMO, 2710 & FNBDZ0, IADRZ0, 2711 & TIMXYM,TIMBF,TIMIO,TIMPQ, 2712 & IOPTRSP,WORK,LWORK) 2713*---------------------------------------------------------------------* 2714* 2715* Purpose: precalculate some intermediates for the F matrix 2716* transformation which depend on left and right vectors: 2717* 2718* X, Y, M, Chi, Yps, P, Q, and the BZeta density 2719* 2720* the results are written to direct acces files 2721* 2722* IOPTRSP = 1 computes zeroth-order version of BZeta density 2723* and P and Q intermediates in MO and backtransformed 2724* X, Y are only used local 2725* 2726* IOPTRSP = 2 computes F matrix version of BZeta density, 2727* and only backtransformed P and Q intermediates 2728* X, Y are written to file 2729* 2730* Written by Christof Haettig November 1998 2731*---------------------------------------------------------------------* 2732#if defined (IMPLICIT_NONE) 2733 IMPLICIT NONE 2734#else 2735# include "implicit.h" 2736#endif 2737#include "priunit.h" 2738#include "ccorb.h" 2739#include "maxorb.h" 2740#include "ccsdsym.h" 2741#include "ccsdinp.h" 2742#include "cciccset.h" 2743#include "ccfield.h" 2744#include "second.h" 2745 2746 INTEGER ISYM0 2747 PARAMETER (ISYM0 = 1) 2748 2749 INTEGER LUXIM, LUYIM, LUBDZ0, LUPQ, LUPQMO, LUMGD, LUZDN, LUMIM 2750 2751 CHARACTER*(*) FILYIM,FILXIM,FILMGD,FILPQ,FILPQMO,FNBDZ0,FILZDN, 2752 & FILMIM 2753 INTEGER LWORK,IDLSTL,IDLSTR,NINT2,LENX,LENY,LENMGD,LENM,IOPTRSP 2754 INTEGER IADRPQ(MXCORB_CC,*), IADRZ0(MXCORB_CC,*), INTMED2(4,*) 2755 INTEGER IADRPQMO(MXCORB_CC,*), LENZDN 2756 2757#if defined (SYS_CRAY) 2758 REAL WORK(LWORK), XLAMDP(*), XLAMDH(*) 2759 REAL TIMXYM, TIMBF, TIMIO, TIMPQ, DTIME 2760 REAL TWO, ONE, ZERO, DDOT, DUMMY, XNORM 2761#else 2762 DOUBLE PRECISION WORK(LWORK), XLAMDP(*), XLAMDH(*) 2763 DOUBLE PRECISION TIMXYM, TIMBF, TIMIO, TIMPQ, DTIME 2764 DOUBLE PRECISION TWO, ONE, ZERO, DUMMY 2765#endif 2766 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 2767 2768 CHARACTER*(10) MODEL 2769 CHARACTER*(3) LISTL, LISTR, LISTLOLD, LISTROLD 2770 LOGICAL NEWLVEC, NEWRVEC 2771 INTEGER IOPT, IOPTY, ISYML, ISYMR, ISYINT 2772 INTEGER KXINT, KYINT, KMINT, KMGDL, KCHI, KEND1, LWRK1, ISYM 2773 INTEGER ISYMK, ISYMC, ISYMI, KOFF1, KOFF2, KOFF3, NRHFK, NVIRC 2774 INTEGER KZETA1, KZETA2, KT1AMP, KT2AMP, IINT2, ISTARTPQ 2775 INTEGER MT1AM, MT2AM, MT2SQ, M3ORHF, ILASTL, ILASTR 2776 INTEGER ILSTSYM, MGLMRH, ISTARTZ0, ISTARTPQMO, IDUMMY 2777 INTEGER MGLMDT, KLAMPA, KLAMHA, KZDENS 2778 2779 CALL QENTER('CCFPREINT') 2780 2781*---------------------------------------------------------------------* 2782* set some symmetries & dimensions and do some initializations: 2783*---------------------------------------------------------------------* 2784 IF (CCS .AND. IOPTRSP.EQ.1) THEN 2785 CALL QEXIT('CCFPREINT') 2786 RETURN 2787 END IF 2788 2789 LENX = 0 2790 LENY = 0 2791 LENMGD = 0 2792 LENZDN = 0 2793 MT1AM = 0 2794 MT2AM = 0 2795 MT2SQ = 0 2796 MGLMRH = 0 2797 MGLMDT = 0 2798 M3ORHF = 0 2799 DO ISYM = 1, NSYM 2800 LENX = MAX(LENX ,NMATIJ(ISYM)) 2801 LENY = MAX(LENY ,NMATAB(ISYM)) 2802 LENMGD = MAX(LENMGD,NT2AOIJ(ISYM)) 2803 LENZDN = MAX(LENZDN,N2BST(ISYM)) 2804 MT1AM = MAX(MT1AM, NT1AM(ISYM)) 2805 MT2AM = MAX(MT2AM, NT2AM(ISYM)) 2806 MT2SQ = MAX(MT2SQ, NT2SQ(ISYM)) 2807 MGLMRH = MAX(MGLMRH,NGLMRH(ISYM)) 2808 MGLMDT = MAX(MGLMDT,NGLMDT(ISYM)) 2809 M3ORHF = MAX(M3ORHF,N3ORHF(ISYM)) 2810 END DO 2811 LENM = M3ORHF 2812 2813 ISTARTZ0 = 1 2814 ISTARTPQ = 1 2815 ISTARTPQMO = 1 2816 2817*---------------------------------------------------------------------* 2818* open direct access files: 2819*---------------------------------------------------------------------* 2820 2821 LUPQ = -1 2822 LUZDN = -1 2823 LUXIM = -1 2824 LUYIM = -1 2825 LUMGD = -1 2826 LUBDZ0 = -1 2827 LUPQMO = -1 2828 LUMIM = -1 2829C 2830 IF (CCSD .OR. CCSDT) THEN 2831 CALL WOPEN2(LUPQ, FILPQ, 64,0) 2832 END IF 2833 2834 IF (IOPTRSP.EQ.2 .AND. (CCS.OR.CC2)) THEN 2835 CALL WOPEN2(LUZDN,FILZDN,64,0) 2836 END IF 2837 2838 IF (IOPTRSP.EQ.2 .AND. (.NOT.CCS)) THEN 2839 CALL WOPEN2(LUXIM,FILXIM,64,0) 2840 CALL WOPEN2(LUYIM,FILYIM,64,0) 2841 IF (.NOT.CC2) CALL WOPEN2(LUMGD,FILMGD,64,0) 2842 END IF 2843 2844 IF (IOPTRSP.EQ.1 .AND. (.NOT.CCS) .AND. (.NOT.CC2)) THEN 2845 CALL WOPEN2(LUBDZ0,FNBDZ0, 64,0) 2846 CALL WOPEN2(LUPQMO,FILPQMO,64,0) 2847 END IF 2848 2849 IF (IOPTRSP.EQ.2 .AND. (.NOT.CCS) .AND. (.NOT.CC2) 2850 & .AND. CCR12) THEN 2851 CALL WOPEN2(LUMIM,FILMIM,64,0) 2852 END IF 2853 2854*---------------------------------------------------------------------* 2855* allocate work space: 2856*---------------------------------------------------------------------* 2857 KZETA1 = 1 2858 KT1AMP = KZETA1 + MT1AM 2859 KEND1 = KT1AMP + MT1AM 2860 2861 IF (CCS.OR.CC2) THEN 2862 KLAMHA = KEND1 2863 KLAMPA = KLAMHA + MGLMDT 2864 KZDENS = KLAMPA + MGLMDT 2865 KEND1 = KZDENS + LENZDN 2866 END IF 2867 2868 IF (.NOT.CCS) THEN 2869 KZETA2 = KEND1 2870 KT2AMP = KZETA2 + MT2SQ 2871 KXINT = KT2AMP + MT2AM 2872 KYINT = KXINT + LENX 2873 KEND1 = KYINT + LENY 2874 END IF 2875 2876 IF (CCSD .OR. CCSDT) THEN 2877 KCHI = KEND1 2878 KMINT = KCHI + MAX(LENX,MGLMRH) 2879 KMGDL = KMINT + M3ORHF 2880 KEND1 = KMGDL + LENMGD 2881 END IF 2882 2883 LWRK1 = LWORK - KEND1 2884 IF (LWRK1 .LT. 0) THEN 2885 CALL QUIT('Insufficient memory in CCFPRE2INT.') 2886 END IF 2887 2888*---------------------------------------------------------------------* 2889* loop over pairs of left and right vectors: 2890*---------------------------------------------------------------------* 2891 LISTLOLD = ' ' 2892 LISTROLD = ' ' 2893 ILASTL = 0 2894 ILASTR = 0 2895 DO IINT2 = 1, NINT2 2896 LISTL = VTABLE(INTMED2(2,IINT2)) 2897 IDLSTL = INTMED2(1,IINT2) 2898 LISTR = VTABLE(INTMED2(4,IINT2)) 2899 IDLSTR = INTMED2(3,IINT2) 2900 ISYML = ILSTSYM(LISTL,IDLSTL) 2901 ISYMR = ILSTSYM(LISTR,IDLSTR) 2902 ISYINT = MULD2H(ISYML,ISYMR) 2903 2904 NEWLVEC = (IINT2.EQ.1) .OR. (LISTL.NE.LISTLOLD) 2905 & .OR. (IDLSTL.NE.ILASTL) 2906 2907 NEWRVEC = (IINT2.EQ.1) .OR. (LISTR.NE.LISTROLD) 2908 & .OR. (IDLSTR.NE.ILASTR) 2909 2910*---------------------------------------------------------------------* 2911* if needed read new lagrangian multiplier and/or amplitude vector: 2912* Note that if new lagrangian multipliers are read, the amplitudes 2913* are overwritten and must then also be read from file. 2914*---------------------------------------------------------------------* 2915 DTIME = SECOND() 2916 2917 IF (NEWLVEC) THEN 2918 IOPT = 1 2919 IF ( (.NOT.(CCS.OR.CC2)) .OR. (CC2.AND.NONHF) ) IOPT = 3 2920 CALL CC_RDRSP(LISTL,IDLSTL,ISYML,IOPT,MODEL, 2921 & WORK(KZETA1),WORK(KT2AMP)) 2922 IF (IOPT.EQ.3) CALL CC_T2SQ(WORK(KT2AMP),WORK(KZETA2),ISYML) 2923 2924 LISTLOLD = LISTL 2925 ILASTL = IDLSTL 2926 END IF 2927 2928 IF (NEWRVEC .OR. NEWLVEC) THEN 2929 IOPT = 1 2930 IF ( (.NOT.(CCS.OR.CC2)) .OR. (CC2.AND.NONHF) ) IOPT = 3 2931 CALL CC_RDRSP(LISTR,IDLSTR,ISYMR,IOPT,MODEL, 2932 & WORK(KT1AMP),WORK(KT2AMP)) 2933 IF ( (IOPT.EQ.3) .AND. (.NOT. LISTR(1:2).EQ.'R0') ) THEN 2934 CALL CCLR_DIASCL(WORK(KT2AMP),TWO,ISYMR) 2935 END IF 2936 2937 LISTROLD = LISTR 2938 ILASTR = IDLSTR 2939 END IF 2940 2941 TIMIO = TIMIO + SECOND() - DTIME 2942 2943*---------------------------------------------------------------------* 2944* for CCS & CC2 calculate the double AO backtransformed Zeta1 2945* vector (the `Zeta' density): 2946*---------------------------------------------------------------------* 2947 IF ((CCS.OR.CC2) .AND. IOPTRSP.EQ.2) THEN 2948 2949 CALL CCLR_LAMTRA(XLAMDP,WORK(KLAMPA),XLAMDH,WORK(KLAMHA), 2950 & WORK(KT1AMP),ISYMR) 2951 2952 CALL DZERO(WORK(KZDENS),LENZDN) 2953 2954 IOPT = 2 2955 CALL CC_ZETADEN(WORK(KZDENS),WORK(KZETA1),ISYML, 2956 & XLAMDP,XLAMDH,ISYM0, 2957 & WORK(KLAMPA),WORK(KLAMHA),ISYMR, 2958 & IOPT,WORK(KEND1),LWRK1) 2959 2960 CALL CC_WVEC(LUZDN,FILZDN,LENZDN,N2BST(ISYINT),IINT2, 2961 & WORK(KZDENS)) 2962 2963 END IF 2964*---------------------------------------------------------------------* 2965* calculate X, Y and (MO) Chi intermediate: 2966* save X and Y on file, Chi is only needed for BFZeta density 2967*---------------------------------------------------------------------* 2968 DTIME = SECOND() 2969 2970 IF ( (.NOT.(CCS.OR.CC2)) .OR. (CC2.AND.NONHF) ) THEN 2971 CALL CC_XI(WORK(KXINT),WORK(KZETA2),ISYML,WORK(KT2AMP),ISYMR, 2972 & WORK(KEND1),LWRK1) 2973 2974 CALL CC_YI(WORK(KYINT),WORK(KZETA2),ISYML,WORK(KT2AMP),ISYMR, 2975 & WORK(KEND1),LWRK1) 2976 END IF 2977 2978 IF (CCSD .OR. CCSDT) THEN 2979 2980 IF (IOPTRSP.EQ.1) THEN 2981 2982 CALL CCLT_CHI(WORK(KZETA1),WORK(KXINT),ISYML, 2983 & XLAMDP,ISYM0,WORK(KCHI),1) 2984 2985 ELSE IF (IOPTRSP.EQ.2) THEN 2986 2987 DO ISYMK = 1, NSYM 2988 ISYMC = MULD2H(ISYMK,ISYMR) 2989 ISYMI = MULD2H(ISYMC,ISYML) 2990 2991 KOFF1 = KT1AMP + IT1AM(ISYMC,ISYMK) 2992 KOFF2 = KZETA1 + IT1AM(ISYMC,ISYMI) 2993 KOFF3 = KCHI + IMATIJ(ISYMK,ISYMI) 2994 2995 NRHFK = MAX(NRHF(ISYMK),1) 2996 NVIRC = MAX(NVIR(ISYMC),1) 2997 2998 CALL DGEMM('T','N',NRHF(ISYMK),NRHF(ISYMI),NVIR(ISYMC), 2999 * -ONE,WORK(KOFF1),NVIRC,WORK(KOFF2),NVIRC, 3000 * ZERO,WORK(KOFF3),NRHFK) 3001 3002 END DO 3003 CALL DAXPY(NMATIJ(ISYINT),-ONE,WORK(KXINT),1,WORK(KCHI),1) 3004 3005 END IF 3006 3007 END IF 3008 3009 TIMXYM = TIMXYM + SECOND() - DTIME 3010 3011 IF (IOPTRSP.EQ.2 .AND. (.NOT.CCS)) THEN 3012 DTIME = SECOND() 3013 CALL CC_WVEC(LUXIM,FILXIM,LENX,NMATIJ(ISYINT),IINT2,WORK(KXINT)) 3014 CALL CC_WVEC(LUYIM,FILYIM,LENY,NMATAB(ISYINT),IINT2,WORK(KYINT)) 3015 TIMIO = TIMIO + SECOND() - DTIME 3016 END IF 3017 3018*---------------------------------------------------------------------* 3019* for CCSD calculate also the M intermediates 3020*---------------------------------------------------------------------* 3021 IF (CCSD .OR. CCSDT) THEN 3022 DTIME = SECOND() 3023 CALL CC_MI(WORK(KMINT),WORK(KZETA2),ISYML,WORK(KT2AMP),ISYMR, 3024 & WORK(KEND1),LWRK1) 3025 TIMXYM = TIMXYM + SECOND() - DTIME 3026 3027 IF (IOPTRSP.EQ.2 .AND. CCR12) THEN 3028 DTIME = SECOND() 3029 CALL CC_WVEC(LUMIM,FILMIM,LENM,N3ORHF(ISYINT),IINT2, 3030 & WORK(KMINT)) 3031 TIMIO = TIMIO + SECOND() - DTIME 3032 END IF 3033 END IF 3034 3035*---------------------------------------------------------------------* 3036* calculate the P and Q intermediates and write them to file: 3037*---------------------------------------------------------------------* 3038 IF (CCSD .OR. CCSDT) THEN 3039 DTIME = SECOND() 3040 3041 IF (IOPTRSP.EQ.1) THEN 3042 3043 IOPTY = 1 3044 CALL CC_PQIMO(WORK(KZETA2),ISYML,WORK(KT2AMP),ISYMR, 3045 & WORK(KYINT),IOPTY,FILPQMO,LUPQMO, 3046 & IADRPQMO,ISTARTPQMO,IINT2, 3047 & WORK(KEND1),LWRK1) 3048 3049 CALL CC_PQIAO(FILPQMO,LUPQMO,IADRPQMO,ISYINT, 3050 & FILPQ,LUPQ,IADRPQ,ISTARTPQ,IINT2, 3051 & XLAMDH,ISYM0,WORK(KEND1),LWRK1) 3052 3053 ELSE IF (IOPTRSP.EQ.2) THEN 3054 3055 IOPTY = 1 3056 CALL CC_PQI2(WORK(KZETA2),ISYML,WORK(KT2AMP),ISYMR, 3057 & WORK(KYINT),IOPTY,FILPQ,LUPQ, 3058 & IADRPQ,ISTARTPQ,IINT2,WORK(KEND1),LWRK1) 3059 3060 END IF 3061 3062 TIMPQ = TIMPQ + SECOND() - DTIME 3063 END IF 3064 3065*---------------------------------------------------------------------* 3066* calculate effective density for left BFZeta intermediate: 3067*---------------------------------------------------------------------* 3068 IF (CCSD .OR. CCSDT) THEN 3069 3070 IF (IOPTRSP.EQ.1) THEN 3071 3072 IOPT = 8 3073 CALL CC_BFDEN(WORK(KZETA2),ISYML,WORK(KMINT),ISYINT, 3074 & XLAMDP, ISYM0, XLAMDP,ISYM0, 3075 & WORK(KCHI),ISYINT, DUMMY, IDUMMY, 3076 & FNBDZ0,LUBDZ0,IADRZ0,ISTARTZ0, 3077 & IINT2,IOPT,.FALSE.,WORK(KEND1),LWRK1) 3078 3079 ELSE IF (IOPTRSP.EQ.2) THEN 3080 3081 DTIME = SECOND() 3082 CALL CC_BFDENF(WORK(KZETA2), ISYML, WORK(KMINT),ISYINT, 3083 & XLAMDP,ISYM0, WORK(KCHI), ISYINT, 3084 & WORK(KT1AMP), ISYMR, WORK(KMGDL), 3085 & WORK(KEND1),LWRK1) 3086 3087 CALL CC_WVEC(LUMGD,FILMGD,LENMGD,NT2AOIJ(ISYINT), 3088 & IINT2,WORK(KMGDL)) 3089 3090 END IF 3091 TIMBF = TIMBF + SECOND() - DTIME 3092 END IF 3093 3094*---------------------------------------------------------------------* 3095* and loop 3096*---------------------------------------------------------------------* 3097 END DO 3098 3099*---------------------------------------------------------------------* 3100* close direct access files and return: 3101*---------------------------------------------------------------------* 3102 IF (CCSD .OR. CCSDT) THEN 3103 CALL WCLOSE2(LUPQ, FILPQ, 'KEEP') 3104 END IF 3105 3106 IF (IOPTRSP.EQ.2 .AND. (CCS.OR.CC2)) THEN 3107 CALL WCLOSE2(LUZDN,FILZDN,'KEEP') 3108 END IF 3109 3110 IF (IOPTRSP.EQ.2 .AND. (.NOT.CCS)) THEN 3111 CALL WCLOSE2(LUXIM,FILXIM,'KEEP') 3112 CALL WCLOSE2(LUYIM,FILYIM,'KEEP') 3113 IF (.NOT.CC2) CALL WCLOSE2(LUMGD,FILMGD,'KEEP') 3114 END IF 3115 3116 IF (IOPTRSP.EQ.1 .AND. (.NOT.CCS) .AND. (.NOT.CC2)) THEN 3117 CALL WCLOSE2(LUBDZ0,FNBDZ0, 'KEEP') 3118 CALL WCLOSE2(LUPQMO,FILPQMO,'KEEP') 3119 END IF 3120 3121 IF (IOPTRSP.EQ.2 .AND. (.NOT.CCS) .AND. (.NOT.CC2) 3122 & .AND. CCR12) THEN 3123 CALL WCLOSE2(LUMIM,FILMIM,'KEEP') 3124 END IF 3125 3126 CALL QEXIT('CCFPREINT') 3127 3128 RETURN 3129 END 3130*=====================================================================* 3131*---------------------------------------------------------------------* 3132c/* Deck CCFOPEN*/ 3133*=====================================================================* 3134 SUBROUTINE CCFOPEN( LUBF, LUFK, LUR, LUFIM, LUGIM, LURHO, 3135 & BFFIL, FKFIL, RFIL, FIMFIL, GIMFIL, RHOFIL, 3136 & LENBF, LENFK, LENR, LENFIM, LENGIM, LENRHO, 3137 & NINT1, NINT2, WORK, LWORK ) 3138*---------------------------------------------------------------------* 3139* Purpose: open files for intermediates in F matrix transformation 3140* which are used under the loop over AO integrals and 3141* need to be initialized at the beginning. 3142* 3143* Written by Christof Haettig, November 1998. 3144*=====================================================================* 3145#if defined (IMPLICIT_NONE) 3146 IMPLICIT NONE 3147#else 3148# include "implicit.h" 3149#endif 3150#include "priunit.h" 3151#include "ccsdinp.h" 3152#include "ccorb.h" 3153#include "ccsdsym.h" 3154 3155* local parameters: 3156 LOGICAL LOCDBG 3157 PARAMETER (LOCDBG = .FALSE.) 3158 3159 CHARACTER*(*) BFFIL, FKFIL, RFIL, FIMFIL, GIMFIL, RHOFIL 3160 INTEGER LUBF, LUFK, LUR, LUFIM, LUGIM, LURHO 3161 INTEGER LENBF, LENFK, LENR, LENFIM, LENGIM, LENRHO 3162 INTEGER NINT1, NINT2, LWORK 3163 INTEGER IINT1, IINT2, ISYM 3164 3165#if defined (SYS_CRAY) 3166 REAL WORK(LWORK) 3167#else 3168 DOUBLE PRECISION WORK(LWORK) 3169#endif 3170 3171 CALL QENTER('CCFOPEN') 3172 3173*---------------------------------------------------------------------* 3174* open files for local intermediates: 3175*---------------------------------------------------------------------* 3176 LUBF = -1 3177 LUR = -1 3178 LUFIM = -1 3179 LUGIM = -1 3180 LUFK = -1 3181 LURHO = -1 3182 IF (.NOT. (CCS.OR.CC2)) THEN 3183 CALL WOPEN2(LUBF, BFFIL, 64,0) 3184 CALL WOPEN2(LUR, RFIL, 64,0) 3185 CALL WOPEN2(LUFIM,FIMFIL,64,0) 3186 ELSE IF (CCS.OR.CC2) THEN 3187 CALL WOPEN2(LUGIM, GIMFIL,64, 0) 3188 CALL WOPEN2(LUFK, FKFIL, 64, 0) 3189 END IF 3190 3191 IF (.NOT.CCS) THEN 3192 CALL WOPEN2(LURHO,RHOFIL,64,0) 3193 END IF 3194 3195*---------------------------------------------------------------------* 3196* calculate a common vector length for BF intermediates and 3197* initialize them with zeros: 3198*---------------------------------------------------------------------* 3199 IF (.NOT. (CCS.OR.CC2)) THEN 3200 LENBF = 0 3201 DO ISYM = 1, NSYM 3202 LENBF = MAX(LENBF,NT2AOIJ(ISYM)) 3203 END DO 3204 3205 IF (LWORK .LT. LENBF) THEN 3206 CALL QUIT('OUT OF MEMORY IN CCFOPEN.') 3207 END IF 3208 3209 CALL DZERO(WORK,LENBF) 3210 3211 DO IINT1 = 1, NINT1 3212 CALL CC_WVEC(LUBF,BFFIL,LENBF,LENBF,IINT1,WORK) 3213 END DO 3214 3215 END IF 3216 3217*---------------------------------------------------------------------* 3218* calculate a common vector length for R intermediates and 3219* initialize them with zeros: 3220*---------------------------------------------------------------------* 3221 IF (.NOT. (CCS.OR.CC2)) THEN 3222 LENR = 0 3223 DO ISYM = 1, NSYM 3224 LENR = MAX(LENR,NEMAT1(ISYM)) 3225 END DO 3226 3227 IF (LWORK .LT. LENR) THEN 3228 CALL QUIT('OUT OF MEMORY IN CCFOPEN.') 3229 END IF 3230 3231 CALL DZERO(WORK,LENR) 3232 3233 DO IINT1 = 1, NINT1 3234 CALL CC_WVEC(LUR,RFIL,LENR,LENR,IINT1,WORK) 3235 END DO 3236 3237 END IF 3238 3239*---------------------------------------------------------------------* 3240* calculate a common vector length for the response Fock matrices and 3241* initialize them with zeros: 3242*---------------------------------------------------------------------* 3243 IF (CCS.OR.CC2) THEN 3244 LENFK = 0 3245 DO ISYM = 1, NSYM 3246 LENFK = MAX(LENFK,N2BST(ISYM)) 3247 END DO 3248 LENGIM = LENFK 3249 3250 IF (LWORK .LT. LENFK) THEN 3251 CALL QUIT('OUT OF MEMORY IN CCFOPEN.') 3252 END IF 3253 3254 CALL DZERO(WORK,LENFK) 3255 3256 DO IINT1 = 1, NINT1 3257 CALL CC_WVEC(LUFK,FKFIL,LENFK,LENFK,IINT1,WORK) 3258 END DO 3259 3260 DO IINT2 = 1, NINT2 3261 CALL CC_WVEC(LUGIM,GIMFIL,LENGIM,LENGIM,IINT2,WORK) 3262 END DO 3263 END IF 3264 3265*---------------------------------------------------------------------* 3266* calculate a common vector length for the RHO1 vectors and 3267* initialize them with zeros: 3268*---------------------------------------------------------------------* 3269 IF (.NOT.CCS) THEN 3270 LENRHO = 0 3271 DO ISYM = 1, NSYM 3272 LENRHO = MAX(LENRHO,NT1AM(ISYM)) 3273 END DO 3274 3275 IF (LWORK .LT. LENRHO) THEN 3276 CALL QUIT('OUT OF MEMORY IN CCFOPEN.') 3277 END IF 3278 3279 CALL DZERO(WORK,LENRHO) 3280 3281 DO IINT2 = 1, NINT2 3282 CALL CC_WVEC(LURHO,RHOFIL,LENRHO,LENRHO,IINT2,WORK) 3283 END DO 3284 END IF 3285 3286*---------------------------------------------------------------------* 3287* calculate a common vector length for the F intermediates and 3288* initialize them with zeros: 3289*---------------------------------------------------------------------* 3290 IF (.NOT. (CCS.OR.CC2)) THEN 3291 LENFIM = 0 3292 DO ISYM = 1, NSYM 3293 LENFIM = MAX(LENFIM,NT1AO(ISYM)) 3294 END DO 3295 3296 IF (LWORK .LT. LENFIM) THEN 3297 CALL QUIT('OUT OF MEMORY IN CCFOPEN.') 3298 END IF 3299 3300 CALL DZERO(WORK,LENFIM) 3301 3302 DO IINT2 = 1, NINT2 3303 CALL CC_WVEC(LUFIM,FIMFIL,LENFIM,LENFIM,IINT2,WORK) 3304 END DO 3305 3306 END IF 3307 3308 CALL QEXIT('CCFOPEN') 3309 3310 RETURN 3311 END 3312*=====================================================================* 3313* END OF SUBROUTINE CCFOPEN * 3314*=====================================================================* 3315*---------------------------------------------------------------------* 3316c/* Deck CCFSAVE*/ 3317*=====================================================================* 3318 SUBROUTINE CCFSAVE(IBATCH, I1HGH, I2HGH, INTMED1, INTMED2, 3319 & KRHO2, LUBF, BFFIL, LENBF, 3320 & KFOCK, LUFK, FKFIL, LENFK, 3321 & KRIM, LUR, RFIL, LENR, 3322 & KFIM, LUFIM,FIMFIL,LENFIM, 3323 & KGIM, LUGIM,GIMFIL,LENGIM, 3324 & KRHO1, LURHO,RHOFIL,LENRHO, 3325 & NINT1, NINT2,WORK, LWORK ) 3326*---------------------------------------------------------------------* 3327* Purpose: save intermediates for F matrix transformation on file 3328* 3329* Written by Christof Haettig, November 1998. 3330*=====================================================================* 3331#if defined (IMPLICIT_NONE) 3332 IMPLICIT NONE 3333#else 3334# include "implicit.h" 3335#endif 3336#include "priunit.h" 3337#include "ccsdinp.h" 3338#include "ccsdsym.h" 3339#include "ccorb.h" 3340#include "cciccset.h" 3341 3342* local parameters: 3343 LOGICAL LOCDBG 3344 PARAMETER (LOCDBG = .FALSE.) 3345 3346 INTEGER LUBF, LUFK, LUR, LUFIM, LUGIM, LURHO 3347 INTEGER LENBF,LENFK,LENR,LENFIM,LENGIM, LENRHO 3348 INTEGER NINT1,NINT2,LWORK,IBATCH 3349 CHARACTER*(*) BFFIL, FIMFIL, GIMFIL, FKFIL, RFIL, RHOFIL 3350 INTEGER I1HGH(0:IBATCH), I2HGH(0:IBATCH) 3351 INTEGER INTMED1(2,NINT1), INTMED2(4,NINT2) 3352 INTEGER KRHO2(NINT1), KFOCK(NINT1), KRIM(NINT1) 3353 INTEGER KFIM(NINT2), KRHO1(NINT2), KGIM(NINT2) 3354 3355 CHARACTER*(3) LIST, LISTL, LISTR 3356 INTEGER IDLST, IDLSTL, IDLSTR, ISYM, ISYML, ISYMR, ISYRES, LEN 3357 INTEGER IINT1, IINT2 3358 3359#if defined (SYS_CRAY) 3360 REAL WORK(LWORK) 3361 REAL XNORM 3362 REAL DDOT 3363#else 3364 DOUBLE PRECISION WORK(LWORK) 3365 DOUBLE PRECISION XNORM 3366 DOUBLE PRECISION DDOT 3367#endif 3368 3369* external function: 3370 INTEGER ILSTSYM 3371 3372 CALL QENTER('CCFSAVE') 3373 3374*---------------------------------------------------------------------* 3375* Fock, BF, and R intermediates: 3376*---------------------------------------------------------------------* 3377 DO IINT1 = I1HGH(IBATCH-1)+1, I1HGH(IBATCH) 3378 LIST = VTABLE(INTMED1(2,IINT1)) 3379 IDLST = INTMED1(1,IINT1) 3380 ISYM = ILSTSYM(LIST,IDLST) 3381 3382* BF intermediate: 3383 IF (.NOT. (CCS .OR. CC2)) THEN 3384 LEN = NT2AOIJ(ISYM) 3385 CALL CC_WVEC(LUBF,BFFIL, LENBF,LEN,IINT1,WORK(KRHO2(IINT1))) 3386 IF (LOCDBG) THEN 3387 XNORM = DDOT(LEN,WORK(KRHO2(IINT1)),1,WORK(KRHO2(IINT1)),1) 3388 WRITE (LUPRI,*) 'Norm of saved BF nb. ',IINT1,' is ', 3389 & XNORM 3390 END IF 3391 END IF 3392 3393* R intermediate: 3394 IF (.NOT. (CCS .OR. CC2)) THEN 3395 LEN = NEMAT1(ISYM) 3396 CALL CC_WVEC(LUR,RFIL,LENR,LEN,IINT1,WORK(KRIM(IINT1))) 3397 IF (LOCDBG) THEN 3398 XNORM = DDOT(LEN,WORK(KRIM(IINT1)),1,WORK(KRIM(IINT1)),1) 3399 WRITE (LUPRI,*) 'Norm of saved RIM nb. ',IINT1,' is ', 3400 & XNORM 3401 END IF 3402 END IF 3403 3404* Fock intermediate: 3405 IF (CCS .OR. CC2) THEN 3406 LEN = N2BST(ISYM) 3407 CALL CC_WVEC (LUFK,FKFIL,LENFK,LEN,IINT1,WORK(KFOCK(IINT1))) 3408 IF (LOCDBG) THEN 3409 XNORM = DDOT(LEN,WORK(KFOCK(IINT1)),1,WORK(KFOCK(IINT1)),1) 3410 WRITE (LUPRI,*) 'Norm of saved FOCK mat. nb. ',IINT1,' is ', 3411 & XNORM 3412 END IF 3413 END IF 3414 3415 END DO 3416 3417*---------------------------------------------------------------------* 3418* RHO1, F and the CCS G intermediate: 3419*---------------------------------------------------------------------* 3420 DO IINT2 = I2HGH(IBATCH-1)+1, I2HGH(IBATCH) 3421 LISTL = VTABLE(INTMED2(2,IINT2)) 3422 LISTR = VTABLE(INTMED2(4,IINT2)) 3423 IDLSTL = INTMED2(1,IINT2) 3424 IDLSTR = INTMED2(3,IINT2) 3425 ISYML = ILSTSYM(LISTL,IDLSTL) 3426 ISYMR = ILSTSYM(LISTR,IDLSTR) 3427 ISYRES = MULD2H(ISYML,ISYMR) 3428 3429 3430* Rho1: 3431 IF (.NOT.CCS) THEN 3432 LEN = NT1AM(ISYRES) 3433 CALL CC_WVEC(LURHO,RHOFIL,LENRHO,LEN,IINT2,WORK(KRHO1(IINT2))) 3434 IF (LOCDBG) THEN 3435 XNORM = DDOT(LEN,WORK(KRHO1(IINT2)),1,WORK(KRHO1(IINT2)),1) 3436 WRITE (LUPRI,*) 'Norm of saved Rho1 nb. ',IINT2,' is ', 3437 & XNORM 3438 END IF 3439 END IF 3440 3441* F intermediate: 3442 IF (.NOT.(CCS.OR.CC2)) THEN 3443 LEN = NT1AO(ISYRES) 3444 CALL CC_WVEC(LUFIM,FIMFIL,LENFIM,LEN,IINT2,WORK(KFIM(IINT2))) 3445 IF (LOCDBG) THEN 3446 XNORM = DDOT(LEN,WORK(KFIM(IINT2)),1,WORK(KFIM(IINT2)),1) 3447 WRITE (LUPRI,*) 'Norm of saved FIM nb. ',IINT2,' is ', 3448 & XNORM 3449 END IF 3450 END IF 3451 3452* the CCS/CC2 G intermediate: 3453 IF (CCS.OR.CC2) THEN 3454 LEN = N2BST(ISYRES) 3455 CALL CC_WVEC(LUGIM,GIMFIL,LENGIM,LEN,IINT2,WORK(KGIM(IINT2))) 3456 IF (LOCDBG) THEN 3457 XNORM = DDOT(LEN,WORK(KGIM(IINT2)),1,WORK(KGIM(IINT2)),1) 3458 WRITE (LUPRI,*) 'Norm of saved GIM nb. ',IINT2,' is ', 3459 & XNORM 3460 END IF 3461 END IF 3462 3463 END DO 3464 3465 3466 CALL QEXIT('CCFSAVE') 3467 3468 RETURN 3469 END 3470*=====================================================================* 3471* END OF SUBROUTINE CCFSAVE 3472*=====================================================================* 3473*---------------------------------------------------------------------* 3474c/* Deck CCFPRE2 */ 3475*=====================================================================* 3476 SUBROUTINE CCFPRE2(INTMED2,ISTART,IEND, 3477 & KFIM, LUFIM,FIMFIL,LENFIM, 3478 & KGIM, LUGIM,GIMFIL,LENGIM, 3479 & KRHO1,LURHO,RHOFIL,LENRHO, 3480 & KMGDL,LUMGD,MGDFIL,LENMGD, 3481 & KZDEN,LUZDN,FILZDN,LENZDN, 3482 & KCTR2,KLAMPA,KLAMHA,XLAMDP,XLAMDH, 3483 & WORK,LWORK,KENDIN,KENDOUT ) 3484*---------------------------------------------------------------------* 3485* Purpose: prepare for calculation of intermediates that depend 3486* on the AO integrals and left and right vectors 3487* 3488* N.B.: this routine allocates work space for CC_FMAT 3489* INPUT end of used space: KENDIN 3490* OUTPUT end of used space: KENDOUT 3491* 3492* Written by Christof Haettig, November 1998. 3493*=====================================================================* 3494#if defined (IMPLICIT_NONE) 3495 IMPLICIT NONE 3496#else 3497# include "implicit.h" 3498#endif 3499#include "priunit.h" 3500#include "ccsdinp.h" 3501#include "ccsdsym.h" 3502#include "ccorb.h" 3503#include "cciccset.h" 3504 3505* local parameters: 3506 LOGICAL LOCDBG 3507 PARAMETER (LOCDBG = .FALSE.) 3508 INTEGER KDUM 3509 PARAMETER (KDUM = +99 999 999) ! dummy address on work space 3510 3511 INTEGER LWORK, KENDIN, KENDOUT 3512 INTEGER ISTART, IEND 3513 INTEGER LUFIM,LENFIM, LURHO,LENRHO, LUMGD,LENMGD 3514 INTEGER LUZDN,LENZDN, LUGIM,LENGIM 3515 INTEGER INTMED2(4,IEND) 3516 INTEGER KLAMPA(IEND), KLAMHA(IEND), KZDEN(IEND), KCTR2(IEND) 3517 INTEGER KFIM(IEND), KRHO1(IEND), KMGDL(IEND), KGIM(IEND) 3518 3519 CHARACTER*(*) FIMFIL, RHOFIL, MGDFIL, FILZDN, GIMFIL 3520 CHARACTER*(3) LISTL, LISTR 3521 CHARACTER*(10) MODEL 3522 INTEGER KT1AMPA, KZETA2 3523 INTEGER LEN, KEND1, IOPT 3524 INTEGER IDLSTL, IDLSTR, ISYML, ISYMR, ISYRES, IINT 3525 3526#if defined (SYS_CRAY) 3527 REAL WORK(LWORK) 3528 REAL XLAMDP(NLAMDT), XLAMDH(NLAMDT), ddot 3529#else 3530 DOUBLE PRECISION WORK(LWORK) 3531 DOUBLE PRECISION XLAMDP(NLAMDT), XLAMDH(NLAMDT), ddot 3532#endif 3533 3534* external functions: 3535 INTEGER ILSTSYM 3536 3537 CALL QENTER('CCFPRE2') 3538 3539*---------------------------------------------------------------------* 3540* begin: 3541*---------------------------------------------------------------------* 3542 KENDOUT = KENDIN 3543 3544 DO IINT = ISTART, IEND 3545 LISTL = VTABLE(INTMED2(2,IINT)) 3546 LISTR = VTABLE(INTMED2(4,IINT)) 3547 IDLSTL = INTMED2(1,IINT) 3548 IDLSTR = INTMED2(3,IINT) 3549 ISYML = ILSTSYM(LISTL,IDLSTL) 3550 ISYMR = ILSTSYM(LISTR,IDLSTR) 3551 ISYRES = MULD2H(ISYML,ISYMR) 3552 3553 IF (CCS) THEN 3554 KCTR2(IINT) = KDUM 3555 KLAMPA(IINT) = KDUM 3556 KLAMHA(IINT) = KDUM 3557 KRHO1(IINT) = KDUM 3558 KMGDL(IINT) = KDUM 3559 KFIM(IINT) = KDUM 3560 KGIM(IINT) = KENDOUT 3561 KZDEN(IINT) = KGIM(IINT) + N2BST(ISYRES) 3562 KENDOUT = KZDEN(IINT) + N2BST(ISYRES) 3563 ELSE IF (CC2) THEN 3564 KMGDL(IINT) = KDUM 3565 KFIM(IINT) = KDUM 3566 KLAMPA(IINT) = KENDOUT 3567 KLAMHA(IINT) = KLAMPA(IINT) + NGLMDT(ISYMR) 3568 KGIM(IINT) = KLAMHA(IINT) + NGLMDT(ISYMR) 3569 KZDEN(IINT) = KGIM(IINT) + N2BST(ISYRES) 3570 KRHO1(IINT) = KZDEN(IINT) + N2BST(ISYRES) 3571 KCTR2(IINT) = KRHO1(IINT) + NT1AM(ISYRES) 3572 KENDOUT = KCTR2(IINT) + NT2SQ(ISYML) 3573 ELSE 3574 KCTR2(IINT) = KDUM 3575 KGIM(IINT) = KDUM 3576 KZDEN(IINT) = KDUM 3577 KLAMPA(IINT) = KENDOUT 3578 KLAMHA(IINT) = KLAMPA(IINT) + NGLMDT(ISYMR) 3579 KRHO1(IINT) = KLAMHA(IINT) + NGLMDT(ISYMR) 3580 KMGDL(IINT) = KRHO1(IINT) + NT1AM(ISYRES) 3581 KFIM(IINT) = KMGDL(IINT) + NT2AOIJ(ISYRES) 3582 KENDOUT = KFIM(IINT) + NT1AO(ISYRES) 3583 END IF 3584 3585 KT1AMPA = KENDOUT 3586 KEND1 = KT1AMPA + NT1AM(ISYMR) 3587 3588 IF (CC2) THEN 3589 KZETA2 = KEND1 3590 KEND1 = KZETA2 + NT2AM(ISYML) 3591 END IF 3592 3593 IF ( (LWORK-KEND1) .LE. 0 ) THEN 3594 CALL QUIT('Insufficient work space in CCFPRE2.') 3595 END IF 3596 3597 3598* recover Zeta density for CCS/CC2: 3599 IF (CCS.OR.CC2) THEN 3600 LEN = N2BST(ISYRES) 3601 CALL CC_RVEC(LUZDN,FILZDN,LENZDN,LEN,IINT,WORK(KZDEN(IINT))) 3602 END IF 3603 3604* recover G intermediate for CCS/CC2: 3605 IF (CCS.OR.CC2) THEN 3606 LEN = N2BST(ISYRES) 3607 CALL CC_RVEC(LUGIM,GIMFIL,LENGIM,LEN,IINT,WORK(KGIM(IINT))) 3608 END IF 3609 3610 3611* recover RHO1: 3612 IF (.NOT.CCS) THEN 3613 LEN = NT1AM(ISYRES) 3614 CALL CC_RVEC(LURHO,RHOFIL,LENRHO,LEN,IINT,WORK(KRHO1(IINT))) 3615 END IF 3616 3617* read singles response vector and calculate response Lambda: 3618 IF (.NOT.CCS) THEN 3619 IOPT = 1 ! read singles response vector 3620 CALL CC_RDRSP(LISTR,IDLSTR,ISYMR,IOPT,MODEL, 3621 & WORK(KT1AMPA),WORK(KDUM) ) 3622 3623 ! calculate response Lambda matrices: 3624 CALL CCLR_LAMTRA(XLAMDP,WORK(KLAMPA(IINT)), 3625 & XLAMDH,WORK(KLAMHA(IINT)), 3626 & WORK(KT1AMPA),ISYMR) 3627 END IF 3628 3629* read doubles Zeta response vector: 3630 IF (CC2) THEN 3631 IOPT = 2 3632 CALL CC_RDRSP(LISTL,IDLSTL,ISYML,IOPT,MODEL, 3633 & WORK(KDUM),WORK(KZETA2)) 3634 CALL CC_T2SQ(WORK(KZETA2),WORK(KCTR2(IINT)),ISYML) 3635 END IF 3636 3637* recover FIM intermediate: 3638 IF (CCSD .OR. CCSDT) THEN 3639 LEN = NT1AO(ISYRES) 3640 CALL CC_RVEC(LUFIM,FIMFIL,LENFIM,LEN,IINT,WORK(KFIM(IINT))) 3641 END IF 3642 3643* recover XMGDL density: 3644 IF (CCSD .OR. CCSDT) THEN 3645 LEN = NT2AOIJ(ISYRES) 3646 CALL CC_RVEC(LUMGD,MGDFIL,LENMGD,LEN,IINT,WORK(KMGDL(IINT))) 3647 END IF 3648 3649 END DO 3650 3651 3652 CALL QEXIT('CCFPRE2') 3653 3654 RETURN 3655 END 3656*=====================================================================* 3657* END OF SUBROUTINE CCFPRE2 3658*=====================================================================* 3659*---------------------------------------------------------------------* 3660c/* Deck CCFINT2 */ 3661*=====================================================================* 3662 SUBROUTINE CCFINT2(IDEL, ISYDEL, ISYCTR, 3663 & XINT, BFRHF, X3OINT, XMGDL, ZDEN, 3664 & FIM, GIM, RHO1, 3665 & LUBFI, FNBFI, IADRBFI, ISTARTBFI, 3666 & LUPQR0, FNPQR0, IADRPQ0, 3667 & LUPQRR, FNPQRR, IADRPQR, 3668 & LUBDZ0, FNBDZ0, IADRZ0, 3669 & XLAMPA, XLAMHA, ISYMA, 3670 & XLAMP0, XLAMH0, WORK, LWORK, 3671 & TIMFCK, TIMBF, TIMGZ, TIMIO, TIMFG ) 3672*---------------------------------------------------------------------* 3673* 3674* Purpose: calculate intermediates for F matrix transformation 3675* which require AO integrals and depend on both the 3676* response amplitude and the Lagrangian multiplier vectors 3677* 3678* BFZeta -- written to file 3679* FIM -- intermediate for 21F term 3680* GIM -- intermediate for CCS case 3681* RHO1 -- contribution to 21G term 3682* (GZIM -- intermediate for E term, added to FIM) 3683* 3684* input: XINT - AO integral distribution 3685* BFRHF - (**|kdel) integrals sorted for CC_BFIF 3686* X3OINT - (ij|kdel) integrals for CC_21H 3687* XMGDL - effective density for CC_BFIF 3688* ZDEN - effective density for GIM intermediate 3689* XLAMPA, XLAMHA - response A Lambda matrices 3690* XLAMP0, XLAMH0 - ordinary zeroth order Lambda matrices 3691* 3692* Written by Christof Haettig, November 1998. 3693* 3694*=====================================================================* 3695#if defined (IMPLICIT_NONE) 3696 IMPLICIT NONE 3697#else 3698# include "implicit.h" 3699#endif 3700#include "priunit.h" 3701#include "ccsdinp.h" 3702#include "ccsdsym.h" 3703#include "maxorb.h" 3704#include "ccorb.h" 3705#include "second.h" 3706 3707* local parameters: 3708 LOGICAL LOCDBG 3709 PARAMETER (LOCDBG = .FALSE.) 3710 INTEGER ISYM0 3711 PARAMETER (ISYM0 = 1) 3712 3713 CHARACTER*(*) FNBFI, FNPQR0, FNPQRR, FNBDZ0 3714 INTEGER LUBFI, LUPQR0, LUPQRR, LUBDZ0, ISTARTBFI 3715 INTEGER IDEL, ISYDEL, ISYMA, LWORK, ISYRES 3716 INTEGER ISYCTR, ISYMGD, ISYABK, ISYZT0, ISYZT1, LEN0, LEN1 3717 INTEGER KDSRHA, KMGD, KEND1, LWRK1, IADR, IOPT 3718 INTEGER KPINT0, KQINT0, KPINT1, KQINT1, KWINT1 3719 INTEGER IADRPQ0(MXCORB_CC),IADRPQR(MXCORB_CC) 3720 INTEGER IADRBFI(MXCORB_CC),IADRZ0(MXCORB_CC) 3721 3722 3723#if defined (SYS_CRAY) 3724 REAL WORK(LWORK), XINT(*), BFRHF(*), X3OINT(*) 3725 REAL XLAMP0(*), XLAMH0(*), XLAMPA(*), XLAMHA(*) 3726 REAL XMGDL(*), ZDEN(*), FIM(*), GIM(*), RHO1(*) 3727 REAL TIMBF, TIMGZ, TIMIO, TIMFG, TIMFCK 3728 REAL DUMMY, DTIME 3729#else 3730 DOUBLE PRECISION WORK(LWORK), XINT(*), BFRHF(*), X3OINT(*) 3731 DOUBLE PRECISION XLAMP0(*), XLAMH0(*), XLAMPA(*), XLAMHA(*) 3732 DOUBLE PRECISION XMGDL(*), ZDEN(*), FIM(*), GIM(*), RHO1(*) 3733 DOUBLE PRECISION TIMBF, TIMGZ, TIMIO, TIMFG, TIMFCK 3734 DOUBLE PRECISION DUMMY, DTIME 3735#endif 3736 3737 3738 CALL QENTER('CCFINT2') 3739 3740 3741*---------------------------------------------------------------------* 3742* begin: 3743*---------------------------------------------------------------------* 3744 ISYRES = MULD2H(ISYCTR,ISYMA) 3745 3746*---------------------------------------------------------------------* 3747* calculate the CCS/CC2 GZeta intermediate: 3748*---------------------------------------------------------------------* 3749 IF (CCS.OR.CC2) THEN 3750 DTIME = SECOND() 3751 CALL CC_AOFOCK(XINT,ZDEN,GIM,WORK,LWORK,IDEL,ISYDEL,.FALSE., 3752 & DUMMY,ISYRES) 3753 TIMFCK = TIMFCK + SECOND() - DTIME 3754 END IF 3755 3756*---------------------------------------------------------------------* 3757* calculate BZeta intermediate: 3758*---------------------------------------------------------------------* 3759 IF (CCSD .OR. CCSDT) THEN 3760 DTIME = SECOND() 3761 CALL CC_BFIF(BFRHF,ISYDEL,XMGDL,ISYRES,LUBFI,FNBFI, 3762 & IADRBFI,ISTARTBFI,IDEL,WORK,LWORK) 3763 TIMBF = TIMBF + SECOND() - DTIME 3764 END IF 3765 3766*---------------------------------------------------------------------* 3767* calculate the GZeta intermediate: 3768*---------------------------------------------------------------------* 3769 IF (CCSD .OR. CCSDT) THEN 3770 DTIME = SECOND() 3771 3772 ISYABK = MULD2H(ISYDEL,ISYMA) 3773 ISYMGD = MULD2H(ISYDEL,ISYCTR) 3774 3775 KDSRHA = 1 3776 KMGD = KDSRHA + NDSRHF(ISYABK) 3777 KEND1 = KMGD + NT2BGD(ISYMGD) 3778 LWRK1 = LWORK - KEND1 3779 3780 IF (LWRK1 .LT. 0) THEN 3781 CALL QUIT('Insufficient memory in CCFINT2. (21F & 21G)') 3782 END IF 3783 3784 CALL CCTRBT(XINT,WORK(KDSRHA),XLAMHA,ISYMA, 3785 & WORK(KEND1),LWRK1,ISYDEL) 3786 3787 IADR = IADRZ0(IDEL) 3788 CALL GETWA2(LUBDZ0,FNBDZ0,WORK(KMGD),IADR,NT2BGD(ISYMGD)) 3789 3790 IOPT = 0 3791 CALL CC_GIM(WORK(KDSRHA),ISYABK,WORK(KMGD),ISYMGD,FIM,IOPT, 3792 & WORK(KEND1),LWRK1) 3793 3794 TIMGZ = TIMGZ + SECOND() - DTIME 3795 END IF 3796 3797*---------------------------------------------------------------------* 3798* calculate 21F term and one part of the 21G term: 3799*---------------------------------------------------------------------* 3800 IF (CCSD .OR. CCSDT) THEN 3801 3802 ISYZT0 = MULD2H(ISYDEL,ISYCTR) 3803 ISYZT1 = MULD2H(ISYZT0,ISYMA) 3804 LEN0 = NT2BCD(ISYZT0) 3805 LEN1 = NT2BCD(ISYZT1) 3806 3807 KPINT0 = 1 3808 KQINT0 = KPINT0 + LEN0 3809 KPINT1 = KQINT0 + LEN0 3810 KQINT1 = KPINT1 + LEN1 3811 KWINT1 = KQINT1 + LEN1 3812 KEND1 = KWINT1 + LEN1 3813 LWRK1 = LWORK - KEND1 3814 3815 IF (LWRK1 .LT. 0) THEN 3816 CALL QUIT('Insufficient memory in CCFINT2. (21F & 21G)') 3817 END IF 3818 3819 DTIME = SECOND() 3820 3821 IADR = IADRPQ0(IDEL) 3822 CALL GETWA2(LUPQR0,FNPQR0,WORK(KPINT0),IADR,LEN0) 3823 IADR = IADRPQ0(IDEL) + LEN0 3824 CALL GETWA2(LUPQR0,FNPQR0,WORK(KQINT0),IADR,LEN0) 3825 3826 IADR = IADRPQR(IDEL) 3827 CALL GETWA2(LUPQRR,FNPQRR,WORK(KPINT1),IADR,LEN1) 3828 IADR = IADRPQR(IDEL) + LEN1 3829 CALL GETWA2(LUPQRR,FNPQRR,WORK(KQINT1),IADR,LEN1) 3830 3831 TIMIO = TIMIO + SECOND() - DTIME 3832 3833 IOPT = 2 3834 DTIME = SECOND() 3835 CALL CC_21I2(FIM,XINT,ISYDEL,DUMMY,0, 3836 * WORK(KPINT0),WORK(KQINT0),ISYZT0, 3837 * WORK(KPINT1),WORK(KQINT1),ISYZT1, 3838 * XLAMP0,XLAMH0,ISYM0,XLAMPA,ISYMA, 3839 * WORK(KEND1),LWRK1,IOPT, 3840 * .TRUE.,.FALSE.,.FALSE.) 3841 TIMFG = TIMFG + SECOND() - DTIME 3842 3843 CALL DZERO(WORK(KWINT1),LEN1) 3844 3845 DTIME = SECOND() 3846 CALL CC_21H(RHO1,ISYRES, WORK(KQINT1), 3847 * WORK(KWINT1),WORK(KPINT1),ISYZT1, 3848 * X3OINT,ISYM0, 3849 * WORK(KEND1),LWRK1,ISYDEL) 3850 TIMFG = TIMFG + SECOND() - DTIME 3851 3852 END IF 3853 3854*---------------------------------------------------------------------* 3855* that's it, return: 3856*---------------------------------------------------------------------* 3857 3858 CALL QEXIT('CCFINT2') 3859C 3860 RETURN 3861 3862 END 3863*=====================================================================* 3864* END OF SUBROUTINE CCFINT2 3865*=====================================================================* 3866*---------------------------------------------------------------------* 3867c/* Deck CCFMAT2 */ 3868*=====================================================================* 3869 SUBROUTINE CCFMAT2(THETA1, THETA2, ISYRES, 3870 & LISTL, IDLSTL, ZETA1, ISYCTR, 3871 & LISTR, IDLSTR, T1AMPA, ISYAMP, 3872 & XLAMPA, XLAMHA, XLAMP0, XLAMH0, 3873 & XINT, YINT, XBARA, YBARA, 3874 & RIMA, A2IM, FOCKA, FOCK0AO, FCKC0, 3875 & LUBFI, FNBFI, IADRBFI, 3876 & LUBF, BFFIL, LENBF, IINTR, 3877 & LUC, CTFIL, LUCBAR, CBAFIL, 3878 & LUD, DTFIL, LUDBAR, DBAFIL, 3879 & IOFFCD, WORK, LWORK ) 3880*---------------------------------------------------------------------* 3881* 3882* Purpose: calculate final contributions to F matrix transformation 3883* from precalculated intermediates 3884* 3885* Written by Christof Haettig, November 1998. 3886* 3887*=====================================================================* 3888 USE PELIB_INTERFACE, ONLY: USE_PELIB 3889#if defined (IMPLICIT_NONE) 3890 IMPLICIT NONE 3891#else 3892# include "implicit.h" 3893#endif 3894#include "priunit.h" 3895#include "dummy.h" 3896#include "ccsdinp.h" 3897#include "ccsdsym.h" 3898#include "ccsections.h" 3899#include "maxorb.h" 3900#include "ccorb.h" 3901#include "ccfield.h" 3902#include "ccslvinf.h" 3903#include "qm3.h" 3904!#include "qmmm.h" 3905 3906* local parameters: 3907 CHARACTER MSGDBG*(17) 3908 PARAMETER (MSGDBG='[debug] CCFMAT2> ') 3909 LOGICAL LOCDBG 3910 PARAMETER (LOCDBG = .FALSE.) 3911 INTEGER ISYM0, LUBF0 3912 PARAMETER (ISYM0 = 1) 3913 3914 CHARACTER*(*) LISTL, LISTR, FNBFI, BFFIL 3915 CHARACTER*(*) CTFIL, CBAFIL, DTFIL, DBAFIL 3916 INTEGER IDLSTL, IDLSTR, ISYRES, ISYCTR, ISYAMP, LWORK 3917 INTEGER LUBFI, IADRBFI(MXCORB_CC), IOFFCD 3918 INTEGER LUC, LUCBAR, LUD, LUDBAR, LUBF, LENBF, IINTR 3919 INTEGER IOPTTCME 3920C 3921C 3922#if defined (SYS_CRAY) 3923 REAL WORK(LWORK) 3924 REAL THETA1(*), THETA2(*), ZETA1(*), T1AMPA(*) 3925 REAL XLAMP0(*), XLAMH0(*), XLAMPA(*), XLAMHA(*) 3926 REAL XINT(*), YINT(*), XBARA(*), YBARA(*) 3927 REAL RIMA(*), A2IM(*), FOCKA(*), FOCK0AO(*), FCKC0(*) 3928 REAL XNORM, DDOT, ZERO, ONE, HALF, FACB, FF 3929#else 3930 DOUBLE PRECISION WORK(LWORK) 3931 DOUBLE PRECISION THETA1(*), THETA2(*), ZETA1(*), T1AMPA(*) 3932 DOUBLE PRECISION XLAMP0(*), XLAMH0(*), XLAMPA(*), XLAMHA(*) 3933 DOUBLE PRECISION XINT(*), YINT(*), XBARA(*), YBARA(*), RIMA(*) 3934 DOUBLE PRECISION A2IM(*), FOCKA(*), FOCK0AO(*), FCKC0(*) 3935 DOUBLE PRECISION XNORM, DDOT, ZERO, ONE, HALF, FACB, FF 3936#endif 3937 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, HALF = 0.5D0) 3938 3939 LOGICAL LGAMMA, LBFZETA, LO3BF, LRCON, LGCON, FCKCON 3940 CHARACTER*(10) MODEL 3941 INTEGER KGZETA, KGIMA, KEND0, LWRK0, ISYMB, ISYMJ, KFOCKA 3942 INTEGER IOPTG, ICON, KEND1, LWRK1, KBZETA, KFOCK 3943 INTEGER IOPT, KGAMMA, KBF, KZETA2, KOFF1, KOFF2, KLIAJB 3944 INTEGER KYPS, KCHI, ISYMA, ISYMD, ISYMK, NVIRA, NVIRD, KOFF3 3945 INTEGER IOPTE, IOPTB, KCDBAR, KEMAT1, KEMAT2 3946 INTEGER KONEHR, KONEHG, KONEH, KA2CON 3947 INTEGER NRHFK, NVIRC, ISYMC, ISYMI, IF, IOPTR12 3948 INTEGER KFIELD, KFIELDR, KFIELDG, ISYLMA 3949 3950 REAL*8, ALLOCATABLE :: FOCKMAT(:), FOCKTEMP(:) 3951 3952 3953 CALL QENTER('CCFMAT2') 3954 3955*---------------------------------------------------------------------* 3956* begin: 3957*---------------------------------------------------------------------* 3958 IF (CCS) THEN 3959 KONEHR = 1 3960 KONEHG = KONEHR + MAX(N2BST(ISYM0),N2BST(ISYAMP)) 3961 KEMAT1 = KONEHG + MAX(N2BST(ISYM0),N2BST(ISYAMP)) 3962 KEMAT2 = KEMAT1 + NMATAB(ISYAMP) 3963 KEND0 = KEMAT2 + NMATIJ(ISYAMP) 3964 ELSE 3965 KGZETA = 1 3966 KGIMA = KGZETA + NT1AO(ISYRES) 3967 KFOCKA = KGIMA + NT1AO(ISYAMP) 3968 KYPS = KFOCKA + NT1AM(ISYAMP) 3969 KEMAT1 = KYPS + NMATAB(ISYRES) 3970 KEMAT2 = KEMAT1 + MAX(NMATAB(ISYAMP),NEMAT1(ISYAMP)) 3971 KONEH = KEMAT2 + NMATIJ(ISYAMP) 3972 KONEHR = KONEH + N2BST(ISYM0) 3973 KONEHG = KONEHR + MAX(N2BST(ISYM0),N2BST(ISYAMP)) 3974 KA2CON = KONEHG + MAX(N2BST(ISYM0),N2BST(ISYAMP)) 3975 KEND0 = KA2CON + NT1AM(ISYRES) 3976 END IF 3977 3978 IF (CC2.AND.NONHF) THEN 3979 KFIELD = KEND0 3980 KFIELDR = KFIELD + N2BST(ISYM0) 3981 KFIELDG = KFIELDR + MAX(N2BST(ISYM0),N2BST(ISYAMP)) 3982 KEND0 = KFIELDG + MAX(N2BST(ISYM0),N2BST(ISYAMP)) 3983 END IF 3984 3985 LWRK0 = LWORK - KEND0 3986 IF (LWRK0.LT.0) THEN 3987 CALL QUIT('Insufficient memory in CCFMAT2. (0)') 3988 END IF 3989 3990*---------------------------------------------------------------------* 3991* for CCS complete here the calculation of the FOCK^A matrix and 3992* put it into ONEHG and ONEHR arrays used for construction of E1 & E2: 3993*---------------------------------------------------------------------* 3994 IF (CCS.OR.CC2) THEN 3995 3996 CALL DCOPY(N2BST(ISYM0),FOCK0AO,1,WORK(KONEHR),1) 3997 CALL DCOPY(N2BST(ISYM0),FOCK0AO,1,WORK(KONEHG),1) 3998 3999 CALL CC_FCKMO(WORK(KONEHR),XLAMPA,XLAMH0,WORK(KEND0),LWRK0, 4000 & ISYM0,ISYAMP,ISYM0) 4001 4002 CALL CC_FCKMO(WORK(KONEHG),XLAMP0,XLAMHA,WORK(KEND0),LWRK0, 4003 & ISYM0,ISYM0,ISYAMP) 4004 4005 CALL DAXPY(N2BST(ISYAMP),ONE,WORK(KONEHG),1,WORK(KONEHR),1) 4006 CALL DAXPY(N2BST(ISYAMP),ONE,FOCKA, 1,WORK(KONEHR),1) 4007 4008 CALL DCOPY(N2BST(ISYAMP),WORK(KONEHR),1,WORK(KONEHG),1) 4009 4010 END IF 4011 4012*---------------------------------------------------------------------* 4013* get one-electron hamiltonian and transform to MO: 4014* if frozen core used, add core-only zeroth-order Fock matrix 4015*---------------------------------------------------------------------* 4016 IF (.NOT.(CCS.OR.CC2)) THEN 4017 CALL CCRHS_ONEAO(WORK(KONEH),WORK(KEND0),LWRK0) 4018 DO IF = 1, NFIELD 4019 FF = EFIELD(IF) 4020 CALL CC_ONEP(WORK(KONEH),WORK(KEND0),LWRK0,FF,1,LFIELD(IF)) 4021 END DO 4022C 4023C--------------------------------------------------------------------------- 4024C Solvent contribution to one-electron integrals. 4025C T^g contribution to transformation. 4026C SLV98,OC, NYQMMM10 KS 4027C------------------------------------------------------------------------------ 4028C 4029 IF (CCSLV .AND. (.NOT. CCMM)) THEN 4030 CALL CCSL_RHSTG(WORK(KONEH),WORK(KEND0),LWRK0) 4031 ENDIF 4032C 4033 IF (CCMM) THEN 4034 IF (.NOT. NYQMMM) THEN 4035 CALL CCMM_RHSTG(WORK(KONEH),WORK(KEND0),LWRK0) 4036 ELSE IF (NYQMMM) THEN 4037 IF (HFFLD) THEN 4038 CALL CCMM_ADDGHF(WORK(KONEH),WORK(KEND0),LWRK0) 4039 ELSE 4040 CALL CCMM_ADDG(WORK(KONEH),WORK(KEND0),LWRK0) 4041 END IF 4042 END IF 4043 ENDIF 4044 IF (USE_PELIB()) THEN 4045 ALLOCATE(FOCKMAT(NNBASX),FOCKTEMP(N2BST(ISYM0))) 4046 IF (HFFLD) THEN 4047 CALL GET_FROM_FILE('FOCKMHF',NNBASX,FOCKMAT) 4048 ELSE 4049 CALL GET_FROM_FILE('FOCKMAT',NNBASX,FOCKMAT) 4050 END IF 4051 CALL DSPTSI(NBAS,FOCKMAT,FOCKTEMP) 4052 CALL DAXPY(N2BST(ISYM0),1.0d0,FOCKTEMP,1,WORK(KONEH),1) 4053 DEALLOCATE(FOCKMAT,FOCKTEMP) 4054 END IF 4055C 4056C----------------------------------------------------------------- 4057 4058 IF (FROIMP.OR.FROEXP) THEN 4059 CALL DAXPY(N2BST(ISYM0),ONE,FCKC0,1,WORK(KONEH),1) 4060 END IF 4061 4062 CALL DCOPY(N2BST(ISYM0),WORK(KONEH),1,WORK(KONEHR),1) 4063 CALL DCOPY(N2BST(ISYM0),WORK(KONEH),1,WORK(KONEHG),1) 4064 4065 CALL CC_FCKMO(WORK(KONEH),XLAMP0,XLAMH0,WORK(KEND0),LWRK0, 4066 & ISYM0,ISYM0,ISYM0) 4067 CALL CC_FCKMO(WORK(KONEHR),XLAMPA,XLAMH0,WORK(KEND0),LWRK0, 4068 & ISYM0,ISYAMP,ISYM0) 4069 CALL CC_FCKMO(WORK(KONEHG),XLAMP0,XLAMHA,WORK(KEND0),LWRK0, 4070 & ISYM0,ISYM0,ISYAMP) 4071 END IF 4072 4073*---------------------------------------------------------------------* 4074* get external fields added after the SCF and transform to MO: 4075*---------------------------------------------------------------------* 4076 IF (CC2.AND.NONHF) THEN 4077 CALL DZERO(WORK(KFIELD),N2BST(ISYM0)) 4078 DO IF = 1, NFIELD 4079 FF = EFIELD(IF) 4080 CALL CC_ONEP(WORK(KFIELD),WORK(KEND0),LWRK0,FF,1,LFIELD(IF)) 4081 END DO 4082 4083C 4084C----------------------------------------------------------------------------------------- 4085C Solvent contribution to one-electron integrals. 4086C T^g contribution to transformation. 4087C SLV98,OC, NYQMMM10 KS 4088C------------------------------------------------------------------------------ 4089C 4090 IF (CCSLV .AND. (.NOT. CCMM)) THEN 4091 CALL CCSL_RHSTG(WORK(KFIELD),WORK(KEND0),LWRK0) 4092 ENDIF 4093C 4094 IF (CCMM) THEN 4095 IF (.NOT. NYQMMM) THEN 4096 CALL CCMM_RHSTG(WORK(KFIELD),WORK(KEND0),LWRK0) 4097 ELSE IF (NYQMMM) THEN 4098 IF (HFFLD) THEN 4099 CALL CCMM_ADDGHF(WORK(KFIELD),WORK(KEND0),LWRK0) 4100 ELSE 4101 CALL CCMM_ADDG(WORK(KFIELD),WORK(KEND0),LWRK0) 4102 END IF 4103 END IF 4104 ENDIF 4105 IF (USE_PELIB()) THEN 4106 ALLOCATE(FOCKMAT(NNBASX),FOCKTEMP(N2BST(ISYM0))) 4107 IF (HFFLD) THEN 4108 CALL GET_FROM_FILE('FOCKMHF',NNBASX,FOCKMAT) 4109 ELSE 4110 CALL GET_FROM_FILE('FOCKMAT',NNBASX,FOCKMAT) 4111 END IF 4112 CALL DSPTSI(NBAS,FOCKMAT,FOCKTEMP) 4113 CALL DAXPY(N2BST(ISYM0),1.0d0,FOCKTEMP,1,WORK(KFIELD),1) 4114 DEALLOCATE(FOCKMAT,FOCKTEMP) 4115 END IF 4116C 4117C------------------------------------------------------------------------------ 4118C 4119 CALL DCOPY(N2BST(ISYM0),WORK(KFIELD),1,WORK(KFIELDR),1) 4120 CALL DCOPY(N2BST(ISYM0),WORK(KFIELD),1,WORK(KFIELDG),1) 4121 4122 CALL CC_FCKMO(WORK(KFIELD),XLAMP0,XLAMH0,WORK(KEND0),LWRK0, 4123 & ISYM0,ISYM0,ISYM0) 4124 CALL CC_FCKMO(WORK(KFIELDR),XLAMPA,XLAMH0,WORK(KEND0),LWRK0, 4125 & ISYM0,ISYAMP,ISYM0) 4126 CALL CC_FCKMO(WORK(KFIELDG),XLAMP0,XLAMHA,WORK(KEND0),LWRK0, 4127 & ISYM0,ISYM0,ISYAMP) 4128 END IF 4129 4130*---------------------------------------------------------------------* 4131* calculate BZeta contribution and GZeta intermediate: 4132*---------------------------------------------------------------------* 4133 IF (CCSD .OR. CCSDT) THEN 4134 KBZETA = KEND0 4135 KEND1 = KBZETA + NT2AO(ISYRES) 4136 LWRK1 = LWORK - KEND1 4137 4138 IF (LWRK1.LT.0) THEN 4139 CALL QUIT('Insufficient memory in CCFMAT2.') 4140 END IF 4141 4142 CALL CC_BFIFSORT(WORK(KBZETA),ISYRES,LUBFI,FNBFI,IADRBFI, 4143 & WORK(KEND1),LWRK1) 4144 4145 CALL DZERO(WORK(KGZETA),NT1AO(ISYRES)) 4146 4147 ICON = 1 4148 IOPTG = 1 4149 LGAMMA = .FALSE. 4150 LO3BF = .FALSE. 4151 LBFZETA = .TRUE. 4152 CALL CC_T2MO3(DUMMY,DUMMY,1,WORK(KBZETA), 4153 & THETA2,DUMMY,WORK(KGZETA),DUMMY, 4154 & XLAMH0,ISYM0,XLAMH0,ISYM0, 4155 & WORK(KEND1),LWRK1,ISYRES, 4156 & ICON,LGAMMA,IOPTG,LO3BF,LBFZETA) 4157 4158 IF (LOCDBG) THEN 4159 WRITE (LUPRI,*) MSGDBG,'norm of THETA1 after B+E term:', 4160 & DDOT(NT1AM(ISYRES),THETA1,1,THETA1,1) 4161 IF (.NOT.CCS) 4162 & WRITE (LUPRI,*) MSGDBG,'norm of THETA2 after B+E term:', 4163 & DDOT(NT2AM(ISYRES),THETA2,1,THETA2,1) 4164C CALL AROUND('THETA2 after B+E term:') 4165C CALL CC_PRP(THETA1,THETA2,ISYRES,0,1) 4166 END IF 4167 4168 END IF 4169 4170*---------------------------------------------------------------------* 4171* calculate the E1' contribution from the GZeta intermediate: 4172*---------------------------------------------------------------------* 4173 IF (CCSD .OR. CCSDT) THEN 4174 4175 CALL CC_T1AM(THETA1,ISYRES,WORK(KGZETA),ISYRES, 4176 & XLAMH0,ISYM0,ONE) 4177 4178 IF (LOCDBG) THEN 4179 CALL AROUND("THETA1 after E1' term:") 4180 CALL CC_PRP(THETA1,THETA2,ISYRES,1,0) 4181 END IF 4182 4183 END IF 4184 4185*---------------------------------------------------------------------* 4186* calculate D2 contribution to Theta1 results vector: 4187*---------------------------------------------------------------------* 4188 IF (.NOT.(CCS.OR.CC2)) THEN 4189 4190 CALL CC_21EFM(THETA1,WORK(KONEH),ISYM0,XINT,YINT,ISYRES, 4191 & WORK(KEND0),LWRK0) 4192 4193 END IF 4194 4195 IF (LOCDBG .AND. .NOT.CCS) THEN 4196 CALL AROUND('THETA after singles excit. D2 term:') 4197 CALL CC_PRP(THETA1,THETA2,ISYRES,1,0) 4198 END IF 4199 4200*---------------------------------------------------------------------* 4201* for CC2/NONHF and CCSD read double excitation part of the Lagrangian 4202* multipliers and expand to a full square matrix: 4203*---------------------------------------------------------------------* 4204 IF ((CC2.AND.NONHF) .OR. CCSD .OR. CCSDT) THEN 4205 KZETA2 = KEND0 4206 KEND1 = KZETA2 + NT2SQ(ISYCTR) 4207 LWRK1 = LWORK - KEND1 4208 4209 IF (LWRK1.LT.NT2AM(ISYCTR)) THEN 4210 CALL QUIT('Insufficient memory in CCFMAT2.') 4211 END IF 4212 4213 IOPT = 2 4214 CALL CC_RDRSP(LISTL,IDLSTL,ISYCTR,IOPT,MODEL,DUMMY,WORK(KEND1)) 4215 CALL CC_T2SQ(WORK(KEND1),WORK(KZETA2),ISYCTR) 4216 4217 ELSE 4218 KEND1 = KEND0 4219 LWRK1 = LWORK - KEND1 4220 END IF 4221 4222*---------------------------------------------------------------------* 4223* calculate the first E2 term contribution from the zero-order BF: 4224*---------------------------------------------------------------------* 4225 IF (CCSD .OR. CCSDT) THEN 4226 KZETA2 = KEND0 4227 KBF = KZETA2 + NT2SQ(ISYCTR) 4228 KEND1 = KBF + 2 * NT2ORT(ISYM0) 4229 LWRK1 = LWORK - KEND1 4230 4231 IF (LWRK1.LT.0) THEN 4232 CALL QUIT('Insufficient memory in CCFMAT2.') 4233 END IF 4234 4235 LUBF0 = -1 4236 CALL GPOPEN(LUBF0,'CC_BFIM','OLD',' ','UNFORMATTED',IDUMMY, 4237 & .FALSE.) 4238 READ(LUBF0) (WORK(KBF-1+I),I=1,2*NT2ORT(ISYM0)) 4239 CALL GPCLOSE(LUBF0,'KEEP') 4240 4241 4242 ICON = 3 4243 IOPTG = 0 4244 LGAMMA = .FALSE. 4245 LO3BF = .FALSE. 4246 LBFZETA = .FALSE. 4247 CALL CC_T2MO3(THETA1,WORK(KZETA2),ISYCTR,WORK(KBF), 4248 & DUMMY,DUMMY,DUMMY,DUMMY, 4249 & XLAMP0,ISYM0,XLAMPA,ISYAMP, 4250 & WORK(KEND1),LWRK1,ISYM0, 4251 & ICON,LGAMMA,IOPTG,LO3BF,LBFZETA) 4252 4253 IF (LOCDBG) THEN 4254 WRITE (LUPRI,*) MSGDBG, 4255 & 'norm(THETA1) after first E2 (LT21C) term:', 4256 & DDOT(NT1AM(ISYRES),THETA1,1,THETA1,1) 4257 IF (.NOT.CCS) 4258 & WRITE (LUPRI,*) MSGDBG, 4259 & 'norm(THETA2) after first E2 (LT21C) term:', 4260 & DDOT(NT2AM(ISYRES),THETA2,1,THETA2,1) 4261C CALL AROUND('THETA after first E2 (LT21C) term:') 4262C CALL CC_PRP(THETA1,THETA2,ISYRES,1,1) 4263 END IF 4264 4265 END IF 4266 4267*---------------------------------------------------------------------* 4268* calculate second E2 term contribution and the GAMMA and G 4269* intermediates from the response BF intermediate: 4270*---------------------------------------------------------------------* 4271 IF (CCSD .OR. CCSDT) THEN 4272 KZETA2 = KEND0 4273 KBF = KZETA2 + NT2SQ(ISYCTR) 4274 KGAMMA = KBF + NT2AOIJ(ISYAMP) 4275 KEND1 = KGAMMA + NGAMMA(ISYAMP) 4276 LWRK1 = LWORK - KEND1 4277 4278 IF (LWRK1.LT.NT2AM(ISYCTR)) THEN 4279 CALL QUIT('Insufficient memory in CCFMAT2.') 4280 END IF 4281 4282 CALL CC_RVEC(LUBF,BFFIL,LENBF,NT2AOIJ(ISYAMP),IINTR,WORK(KBF)) 4283 4284 CALL DZERO(WORK(KGAMMA),NGAMMA(ISYAMP)) 4285 CALL DZERO(WORK(KGIMA), NT1AO(ISYAMP)) 4286 4287 4288 ICON = 3 4289 IOPTG = 2 4290 LGAMMA = .TRUE. 4291 LO3BF = .TRUE. 4292 LBFZETA = .FALSE. 4293 CALL CC_T2MO3(THETA1,WORK(KZETA2),ISYCTR,WORK(KBF), 4294 & DUMMY,WORK(KGAMMA),WORK(KGIMA),DUMMY, 4295 & XLAMP0,ISYM0,XLAMP0,ISYM0, 4296 & WORK(KEND1),LWRK1,ISYAMP, 4297 & ICON,LGAMMA,IOPTG,LO3BF,LBFZETA) 4298 4299 IF (LOCDBG) THEN 4300 WRITE (LUPRI,*) MSGDBG, 4301 & 'norm(THETA1) after second E2 (LT21C) term:', 4302 & DDOT(NT1AM(ISYRES),THETA1,1,THETA1,1) 4303 WRITE (LUPRI,*) MSGDBG, 4304 & 'norm(THETA2) after second E2 (LT21C) term:', 4305 & DDOT(NT2AM(ISYRES),THETA2,1,THETA2,1) 4306 WRITE (LUPRI,*) MSGDBG,'the GIMA intermediate:' 4307 WRITE (LUPRI,'(5F12.6)') (WORK(KGIMA-1+I),I=1,NT1AO(ISYAMP)) 4308C CALL AROUND('THETA after second E2 (LT21C) term:') 4309C CALL CC_PRP(THETA1,THETA2,ISYRES,1,1) 4310 XNORM = DDOT(NGAMMA(ISYAMP),WORK(KGAMMA),1,WORK(KGAMMA),1) 4311 WRITE (LUPRI,*) 'Norm(GAMMA) = ',XNORM 4312 END IF 4313 4314 END IF 4315 4316*---------------------------------------------------------------------* 4317* for CCS/CCSD calculate here E1/E1* and E2 intermediates as nedded 4318* for the B2 contribution from the G and YBAR intermediates: 4319*---------------------------------------------------------------------* 4320 LRCON = .FALSE. 4321 LGCON = .FALSE. 4322 FCKCON = .TRUE. 4323 IOPT = 1 4324 4325 IF (CCSD .OR. CCSDT) LGCON = .TRUE. 4326 4327 CALL CC_EIM(WORK(KEMAT1),WORK(KEMAT2), 4328 & DUMMY,DUMMY,WORK(KGIMA),DUMMY, 4329 & WORK(KONEHR),WORK(KONEHG), 4330 & XLAMH0,XLAMP0,ISYM0,DUMMY,DUMMY,IDUMMY, 4331 & FCKCON,LRCON,LGCON,.FALSE.,IOPT,ISYAMP) 4332 4333 IF (CC2) THEN 4334 CALL DAXPY(NMATAB(ISYAMP),ONE,YBARA,1,WORK(KEMAT1),1) 4335 CALL DAXPY(NMATIJ(ISYAMP),ONE,XBARA,1,WORK(KEMAT2),1) 4336 END IF 4337 4338 IF (CCSD .OR. CCSDT) THEN 4339 CALL DAXPY(NMATAB(ISYAMP),ONE,YBARA,1,WORK(KEMAT1),1) 4340 END IF 4341 4342 IF (LOCDBG) THEN 4343 CALL AROUND('E-intermediates out from CC_EIM:') 4344 CALL CC_PREI(WORK(KEMAT1),WORK(KEMAT2),ISYAMP,1) 4345 END IF 4346 4347*---------------------------------------------------------------------* 4348* calculate B2 contribution to Theta1 results vector: 4349*---------------------------------------------------------------------* 4350 CALL CCLR_E1C1(THETA1,ZETA1,WORK(KEMAT1),WORK(KEMAT2), 4351 & WORK(KEND1),LWRK1,ISYCTR,ISYAMP,'T') 4352 4353 IF (LOCDBG) THEN 4354 CALL AROUND('THETA1 B2 term:') 4355 CALL CC_PRP(THETA1,THETA2,ISYRES,1,0) 4356 END IF 4357 4358*---------------------------------------------------------------------* 4359* for CCSD calculate here the E intermed. from G, R and YBAR intermed.: 4360*---------------------------------------------------------------------* 4361 IF (.NOT.(CCS.OR.CC2)) THEN 4362 4363 LRCON = .TRUE. 4364 LGCON = .TRUE. 4365 FCKCON = .TRUE. 4366 IOPT = 1 4367 CALL CC_EIM(WORK(KEMAT1),WORK(KEMAT2), 4368 & RIMA,DUMMY,WORK(KGIMA),DUMMY, 4369 & WORK(KONEHR),WORK(KONEHG), 4370 & XLAMH0,XLAMP0,ISYM0,DUMMY,DUMMY,IDUMMY, 4371 & FCKCON,LRCON,LGCON,.FALSE.,IOPT,ISYAMP) 4372 4373 CALL DAXPY(NMATAB(ISYAMP),ONE,YBARA,1,WORK(KEMAT1),1) 4374 4375 IF (LOCDBG) THEN 4376 CALL AROUND('E-intermediates out from CC_EIM:') 4377 CALL CC_PREI(WORK(KEMAT1),WORK(KEMAT2),ISYAMP,1) 4378 END IF 4379 4380 END IF 4381 4382*---------------------------------------------------------------------* 4383* calculate the A term and add to Theta2 result vector: 4384*---------------------------------------------------------------------* 4385 IF (CCSD .OR. CCSDT) THEN 4386 CALL CC_GAMMA2(WORK(KGAMMA),WORK(KEMAT2),ISYAMP) 4387 4388 IOPT = 2 4389 CALL CCRHS_A(THETA2,WORK(KZETA2),WORK(KGAMMA), 4390 & WORK(KEND1),LWRK1,ISYAMP,ISYCTR,IOPT) 4391 4392 IF (LOCDBG) THEN 4393 WRITE (LUPRI,*) MSGDBG, 4394 & 'norm(THETA1) after doubles excit A term:', 4395 & DDOT(NT1AM(ISYRES),THETA1,1,THETA1,1) 4396 IF (.NOT.CCS) 4397 & WRITE (LUPRI,*) MSGDBG, 4398 & 'norm(THETA2) after doubles excit A term:', 4399 & DDOT(NT2AM(ISYRES),THETA2,1,THETA2,1) 4400C CALL AROUND('THETA after doubles excit. A term:') 4401C CALL CC_PRP(THETA1,THETA2,ISYRES,1,1) 4402 END IF 4403 4404 END IF 4405 4406*---------------------------------------------------------------------* 4407* read (ia|jb) integrals and calculate L(ia|jb) in place: 4408* (needed for H term and for A2 term) 4409*---------------------------------------------------------------------* 4410 IF (.NOT.CCS) THEN 4411 KLIAJB = KEND0 4412 KEND1 = KLIAJB + NT2AM(ISYM0) 4413 LWRK1 = LWORK - KEND1 4414 4415 IF (LWRK1.LT.0) THEN 4416 CALL QUIT('Insufficient memory in CCFMAT2.') 4417 END IF 4418 4419 CALL CCG_RDIAJB(WORK(KLIAJB),NT2AM(ISYM0)) 4420 IOPTTCME = 1 4421 CALL CCSD_TCMEPK(WORK(KLIAJB),ONE,ISYM0,IOPTTCME) 4422 END IF 4423 4424*---------------------------------------------------------------------* 4425* calculate A2 contribution: 4426*---------------------------------------------------------------------* 4427 IF (.NOT.CCS) THEN 4428 4429 CALL CCG_LXD(WORK(KA2CON),ISYRES,A2IM,ISYRES, 4430 & WORK(KLIAJB),ISYM0,0) 4431 4432 CALL DAXPY(NT1AM(ISYRES),ONE,WORK(KA2CON),1,THETA1,1) 4433 4434 IF (LOCDBG) THEN 4435 CALL AROUND('THETA1 A2 term:') 4436 CALL CC_PRP(THETA1,THETA2,ISYRES,1,0) 4437 CALL AROUND('A2IM intermediate:') 4438 CALL CC_PRP(A2IM,THETA2,ISYRES,1,0) 4439 CALL AROUND('A2CON contribution:') 4440 CALL CC_PRP(WORK(KA2CON),THETA2,ISYRES,1,0) 4441 END IF 4442 4443 END IF 4444 4445*---------------------------------------------------------------------* 4446* calculate H term: 4447*---------------------------------------------------------------------* 4448 IF (.NOT.CCS) THEN 4449 KFOCK = KEND1 4450 KCHI = KFOCK + N2BST(ISYAMP) 4451 KEND1 = KCHI + NMATIJ(ISYRES) 4452 LWRK1 = LWORK - KEND1 4453 4454 IF (LWRK1.LT.0) THEN 4455 CALL QUIT('Insufficient memory in CCFMAT2.') 4456 END IF 4457 4458C ------------------------ 4459C calculate F^A_jb matrix: 4460C ------------------------ 4461 IOPT = 0 4462 CALL CCG_LXD(WORK(KFOCKA),ISYAMP,T1AMPA,ISYAMP, 4463 & WORK(KLIAJB),ISYM0,IOPT) 4464 4465 IF (LOCDBG) THEN 4466 CALL AROUND('FOCKA_jb matrix:') 4467 CALL CC_PRP(WORK(KFOCKA),WORK,ISYRES,1,0) 4468 END IF 4469 4470CTST 4471 CALL DZERO(WORK(KFOCK),N2BST(ISYAMP)) 4472CTST 4473 4474C ----------------------------------------- 4475C resort Fock^A_jb into a full Fock matrix: 4476C ----------------------------------------- 4477 DO ISYMB = 1, NSYM 4478 ISYMJ = MULD2H(ISYAMP,ISYMB) 4479 4480 DO J = 1, NRHF(ISYMJ) 4481 DO B = 1, NVIR(ISYMB) 4482 KOFF1 = IFCVIR(ISYMJ,ISYMB) + NORB(ISYMJ)*(B-1) + J 4483 KOFF2 = IT1AM(ISYMB,ISYMJ) + NVIR(ISYMB)*(J-1) + B 4484 WORK(KFOCK-1+KOFF1) = WORK(KFOCKA-1+KOFF2) 4485 END DO 4486 END DO 4487 4488 END DO 4489 4490C --------------------------------------- 4491C calculate the simple Fock contribution: 4492C --------------------------------------- 4493 CALL CC_L1FCK(THETA2,ZETA1,WORK(KFOCK),ISYCTR,ISYAMP, 4494 & WORK(KEND1),LWRK1) 4495 4496 IF (LOCDBG) THEN 4497 WRITE (LUPRI,*) MSGDBG, 4498 & 'norm(THETA1) after doubles excit H term:', 4499 & DDOT(NT1AM(ISYRES),THETA1,1,THETA1,1) 4500 IF (.NOT.CCS) 4501 & WRITE (LUPRI,*) MSGDBG, 4502 & 'norm(THETA2) after doubles excit H term:', 4503 & DDOT(NT2AM(ISYRES),THETA2,1,THETA2,1) 4504C CALL AROUND('THETA after doubles exci. H term (Fock part):') 4505C CALL CC_PRP(THETA1,THETA2,ISYRES,1,1) 4506 XNORM = DDOT(NT1AM(ISYCTR),ZETA1,1,ZETA1,1) 4507 WRITE (LUPRI,*) 'Norm of ZETA1:',XNORM 4508 END IF 4509 4510C ---------------------------------------------- 4511C calculate the transposed Ypsilon intermediate: 4512C (for CC2 exclude the contribution from YINT) 4513C ---------------------------------------------- 4514 IF (CC2) THEN 4515 CALL DZERO(WORK(KYPS),NMATAB(ISYRES)) 4516 ELSE 4517 CALL DCOPY(NMATAB(ISYRES),YINT,1,WORK(KYPS),1) 4518 END IF 4519 4520 DO ISYMA = 1, NSYM 4521 4522 ISYMD = MULD2H(ISYMA,ISYRES) 4523 ISYMK = MULD2H(ISYMA,ISYCTR) 4524 4525 KOFF1 = IT1AM(ISYMD,ISYMK) + 1 4526 KOFF2 = IT1AM(ISYMA,ISYMK) + 1 4527 KOFF3 = KYPS + IMATAB(ISYMD,ISYMA) 4528 4529 NVIRD = MAX(NVIR(ISYMD),1) 4530 NVIRA = MAX(NVIR(ISYMA),1) 4531 4532 CALL DGEMM('N','T',NVIR(ISYMD),NVIR(ISYMA),NRHF(ISYMK), 4533 & ONE,T1AMPA(KOFF1),NVIRD,ZETA1(KOFF2),NVIRA, 4534 & ONE,WORK(KOFF3),NVIRD) 4535 END DO 4536 4537 CALL DSCAL(NMATAB(ISYRES),HALF,WORK(KYPS),1) 4538 4539C ------------------------------------------------------ 4540C calculate the Chi intermediate: 4541C (for CCSD this contribution is included in the B term) 4542C ------------------------------------------------------ 4543 CALL DZERO(WORK(KCHI),NMATIJ(ISYRES)) 4544 4545 IF (CC2) THEN 4546 4547 DO ISYMK = 1, NSYM 4548 ISYMC = MULD2H(ISYMK,ISYAMP) 4549 ISYMI = MULD2H(ISYMC,ISYCTR) 4550 4551 KOFF1 = IT1AM(ISYMC,ISYMK) + 1 4552 KOFF2 = IT1AM(ISYMC,ISYMI) + 1 4553 KOFF3 = KCHI + IMATIJ(ISYMK,ISYMI) 4554 4555 NRHFK = MAX(NRHF(ISYMK),1) 4556 NVIRC = MAX(NVIR(ISYMC),1) 4557 4558 CALL DGEMM('T','N',NRHF(ISYMK),NRHF(ISYMI),NVIR(ISYMC), 4559 * ONE,T1AMPA(KOFF1),NVIRC,ZETA1(KOFF2),NVIRC, 4560 * ONE,WORK(KOFF3),NRHFK) 4561 4562 END DO 4563 4564 CALL DSCAL(NMATIJ(ISYRES),HALF,WORK(KCHI),1) 4565 4566 END IF 4567 4568C -------------------------------- 4569C contract with (jb|id) integrals: 4570C -------------------------------- 4571 CALL CC_22EC(THETA2,WORK(KLIAJB),WORK(KCHI),WORK(KYPS), 4572 & ISYRES,WORK(KEND1),LWRK1) 4573 4574 IF (LOCDBG) THEN 4575 WRITE (LUPRI,*) MSGDBG, 4576 & 'norm(THETA1) after doubles excit H term:', 4577 & DDOT(NT1AM(ISYRES),THETA1,1,THETA1,1) 4578 IF (.NOT.CCS) 4579 & WRITE (LUPRI,*) MSGDBG, 4580 & 'norm(THETA2) after doubles excit H term:', 4581 & DDOT(NT2AM(ISYRES),THETA2,1,THETA2,1) 4582C CALL AROUND('THETA after doubles excit. H term (2. part):') 4583C CALL CC_PRP(THETA1,THETA2,ISYRES,1,1) 4584 CALL AROUND('the Ypsilon/Chi intermediates:') 4585 CALL CC_PREI(WORK(KYPS),WORK(KCHI),ISYRES,1) 4586 CALL AROUND('the Y and X intermediates:') 4587 CALL CC_PREI(YINT,XINT,ISYRES,1) 4588C CALL AROUND('the integrals:') 4589C CALL CC_PRSQ(WORK,WORK(KLIAJB),ISYM0,0,1) 4590 END IF 4591 4592 END IF 4593 4594*---------------------------------------------------------------------* 4595* calculate the D & C term contributions: 4596*---------------------------------------------------------------------* 4597 IF (.NOT. (CCS.OR.CC2)) THEN 4598 4599 KCDBAR = KEND0 4600 KZETA2 = KCDBAR + NT2SQ(ISYAMP) 4601 KEND1 = KZETA2 + NT2AM(ISYCTR) 4602 LWRK1 = LWORK - KEND1 4603 4604 IF (LWRK1 .LT. NT2AM(ISYCTR)) THEN 4605 CALL QUIT('Insufficient memory in CCFMAT2.') 4606 END IF 4607 4608 IOPT = 2 4609 CALL CC_RDRSP(LISTL,IDLSTL,ISYCTR,IOPT,MODEL, 4610 & DUMMY,WORK(KZETA2)) 4611 4612C ----------------------------------------------- 4613C complete the calculation of the D intermediate: 4614C (including the E1 contribution to the diagonal) 4615C ----------------------------------------------- 4616 IOPT = 3 4617 IOPTB = 1 4618 IOPTE = 1 4619 FACB = ONE 4620 CALL CCRHS_DIO3(WORK(KCDBAR),DUMMY,XLAMH0, 4621 & WORK(KEND1),LWRK1,1,ISYAMP, 4622 & LUD,DTFIL,LUC,CTFIL,IINTR,IOPT, 4623 & IOPTB,LUDBAR,DBAFIL,IOFFCD,FACB, 4624 & IOPTE,WORK(KEMAT1)) 4625 4626C ----------------------------------------------- 4627C contract with the lagrangian multiplier vector: 4628C ----------------------------------------------- 4629 ioptr12 = 0 4630 CALL CC_CD('D',-1,ioptr12,THETA2,ISYRES,WORK(KZETA2),ISYCTR, 4631 & WORK(KCDBAR),ISYAMP,WORK(KEND1),LWRK1) 4632 4633 IF (LOCDBG) THEN 4634 WRITE (LUPRI,*) MSGDBG, 4635 & 'norm(THETA1) after doubles excit D term:', 4636 & DDOT(NT1AM(ISYRES),THETA1,1,THETA1,1) 4637 WRITE (LUPRI,*) MSGDBG, 4638 & 'norm(THETA2) after doubles excit D term:', 4639 & DDOT(NT2AM(ISYRES),THETA2,1,THETA2,1) 4640 WRITE (LUPRI,*) MSGDBG, 4641 & 'norm(D int.) after doubles excit D term:', 4642 & DDOT(NT2SQ(ISYAMP),WORK(KCDBAR),1,WORK(KCDBAR),1) 4643C CALL AROUND('THETA2 after D term:') 4644C CALL CC_PRP(THETA1,THETA2,ISYRES,0,1) 4645 END IF 4646 4647C ----------------------------------------------- 4648C complete the calculation of the C intermediate: 4649C (including the E1 contribution to the diagonal) 4650C ----------------------------------------------- 4651 IOPT = 3 4652 IOPTB = 1 4653 IOPTE = 1 4654 FACB = ONE 4655 CALL CCRHS_CIO3(WORK(KCDBAR),DUMMY,XLAMH0, 4656 & WORK(KEND1),LWRK1,1,ISYAMP, 4657 & LUC,CTFIL,IINTR,IOPT, 4658 & IOPTB,LUCBAR,CBAFIL,IOFFCD,FACB, 4659 & IOPTE,WORK(KEMAT1)) 4660 4661C ----------------------------------------------- 4662C contract with the lagrangian multiplier vector: 4663C ----------------------------------------------- 4664 ioptr12 = 0 4665 CALL CC_CD('C',-1,ioptr12,THETA2,ISYRES,WORK(KZETA2),ISYCTR, 4666 & WORK(KCDBAR),ISYAMP,WORK(KEND1),LWRK1) 4667 4668 IF (LOCDBG) THEN 4669 WRITE (LUPRI,*) MSGDBG, 4670 & 'norm(THETA1) after doubles excit C term:', 4671 & DDOT(NT1AM(ISYRES),THETA1,1,THETA1,1) 4672 WRITE (LUPRI,*) MSGDBG, 4673 & 'norm(THETA2) after doubles excit C term:', 4674 & DDOT(NT2AM(ISYRES),THETA2,1,THETA2,1) 4675 WRITE (LUPRI,*) MSGDBG, 4676 & 'norm(C int.) after doubles excit C term:', 4677 & DDOT(NT2SQ(ISYAMP),WORK(KCDBAR),1,WORK(KCDBAR),1) 4678C CALL AROUND('THETA2 after C term:') 4679C CALL CC_PRP(THETA1,THETA2,ISYRES,0,1) 4680 END IF 4681 4682 END IF 4683 4684*---------------------------------------------------------------------* 4685* that's it, return: 4686*---------------------------------------------------------------------* 4687 4688 CALL QEXIT('CCFMAT2') 4689 4690 RETURN 4691 4692 END 4693*=====================================================================* 4694* END OF SUBROUTINE CCFMAT2 4695*=====================================================================* 4696*=====================================================================* 4697 SUBROUTINE CC_21GMO(THETA1,ISYRES,XILJD,ISYINT, 4698 & LUPQMO,FILPQMO,IADRPQMO,ISYPQI, 4699 & WORK,LWORK) 4700*---------------------------------------------------------------------* 4701* 4702* Purpose: Calculate a 21G like contribution to the F matrix transf. 4703* from (il|jd) integrals and P & Q intermediates in MO 4704* 4705* Theta1 = Theta1 - P^d_alj g^X_iljd + 1/2 (P+Q)^d_alj g^X_jlid 4706* 4707* symmetries: ISYRES - THETA1 result vector 4708* ISYINT - (il|jd) integrals 4709* ISYPQI - P and Q intermediates on file 4710* 4711* Christof Haettig, November 1998 4712*---------------------------------------------------------------------* 4713 IMPLICIT NONE 4714#include "priunit.h" 4715#include "ccorb.h" 4716#include "ccsdsym.h" 4717 4718 CHARACTER*(*) FILPQMO 4719 INTEGER LWORK, ISYRES, ISYINT, ISYPQI, LUPQMO, IADRPQMO(*) 4720 4721#if defined (SYS_CRAY) 4722 REAL XILJD(*), THETA1(*), WORK(LWORK) 4723#else 4724 DOUBLE PRECISION XILJD(*), THETA1(*), WORK(LWORK) 4725#endif 4726 4727 INTEGER IVIRD, ISYALJ, ISYJLI, ISYMD, IADR, LEN, KOFF1 4728 INTEGER KPINT, KWINT, KQINT, KEND1, LWRK1 4729 4730 CALL QENTER('CC_21GMO') 4731 4732*---------------------------------------------------------------------* 4733* check symmetries: 4734*---------------------------------------------------------------------* 4735 IF (MULD2H(ISYINT,ISYPQI).NE.ISYRES) THEN 4736 CALL QUIT('Symmetry mismatch in CC_21GMO.') 4737 END IF 4738 4739*---------------------------------------------------------------------* 4740* contract P and Q intermediates with intgrals in a loop over the 4741* virtual index D 4742*---------------------------------------------------------------------* 4743 DO ISYMD = 1, NSYM 4744 ISYALJ = MULD2H(ISYPQI,ISYMD) 4745 ISYJLI = MULD2H(ISYINT,ISYMD) 4746 4747 LEN = NT2BCD(ISYALJ) 4748 4749 KPINT = 1 4750 KQINT = KPINT + LEN 4751 KWINT = KQINT + LEN 4752 KEND1 = KWINT + LEN 4753 LWRK1 = LWORK - KEND1 4754 4755 IF (LWRK1 .LT. 0) THEN 4756 CALL QUIT('Insufficient memory in CC_21GMO.') 4757 END IF 4758 4759 CALL DZERO(WORK(KWINT),LEN) 4760 4761 DO D = 1, NVIR(ISYMD) 4762 IVIRD = IVIR(ISYMD) + D 4763 IADR = IADRPQMO(IVIRD) 4764 CALL GETWA2(LUPQMO,FILPQMO,WORK(KPINT),IADR, LEN) 4765 CALL GETWA2(LUPQMO,FILPQMO,WORK(KQINT),IADR+LEN,LEN) 4766 4767 KOFF1 = ISJIKA(ISYJLI,ISYMD) + NMAIJK(ISYJLI)*(D-1) + 1 4768 4769 CALL CC_21H(THETA1,ISYRES,WORK(KQINT),WORK(KWINT), 4770 & WORK(KPINT),ISYALJ,XILJD(KOFF1),ISYINT, 4771 & WORK(KEND1),LWRK1,ISYMD) 4772 END DO 4773 4774 END DO 4775 4776 CALL QEXIT('CC_21GMO') 4777 4778 RETURN 4779 END 4780*=====================================================================* 4781* END OF SUBROUTINE CC_21GMO * 4782*=====================================================================* 4783*---------------------------------------------------------------------* 4784c/* Deck CC_ZETADEN */ 4785*=====================================================================* 4786 SUBROUTINE CC_ZETADEN(ZDEN, ZETA1, ISYCTR, 4787 & XLAMP1, XLAMH1, ISYXL1, 4788 & XLAMP2, XLAMH2, ISYXL2, 4789 & IOPT, WORK, LWORK) 4790*---------------------------------------------------------------------* 4791* 4792* Purpose: calculate `Zeta'-density for F-matrix: a double AO 4793* backtransformed Zeta1 vector 4794* 4795* IOPT = 1 - Jacobian: XLAMP1,XLAMH1 same as XLAMP2/XLAMH2 4796* IOPT = 2 - F matrix: XLAMP1,XLAMH1 different from XLAMP2/XLAMH2 4797* 4798* Written by Christof Haettig, November 1998 4799* 4800*=====================================================================* 4801#if defined (IMPLICIT_NONE) 4802 IMPLICIT NONE 4803#else 4804# include "implicit.h" 4805#endif 4806#include "priunit.h" 4807#include "ccsdsym.h" 4808#include "ccorb.h" 4809 4810 LOGICAL LOCDBG 4811 PARAMETER (LOCDBG = .FALSE.) 4812 4813 INTEGER IOPT, ISYCTR, ISYXL1, ISYXL2, LWORK 4814 4815#if defined (SYS_CRAY) 4816 REAL ZDEN(*), ZETA1(*), WORK(LWORK 4817 REAL XLAMP1(*), XLAMH1(*), XLAMP2(*), XLAMH2(*) 4818 REAL ONE, ZERO 4819#else 4820 DOUBLE PRECISION ZDEN(*), ZETA1(*), WORK(LWORK) 4821 DOUBLE PRECISION XLAMP1(*), XLAMH1(*), XLAMP2(*), XLAMH2(*) 4822 DOUBLE PRECISION ONE, ZERO 4823#endif 4824 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0) 4825 4826 INTEGER ISYM12, ISALBE, ISYMAL, ISYMBE, ISYMB, ISYMK 4827 INTEGER KOFF1, KOFF2, KOFF3, NTOTBE, NTOTAL, NTOTB 4828 4829 CALL QENTER('CC_ZETADEN') 4830 4831*---------------------------------------------------------------------* 4832* set symmetries: 4833*---------------------------------------------------------------------* 4834 4835 ISYM12 = MULD2H(ISYXL1,ISYXL2) 4836 ISALBE = MULD2H(ISYCTR,ISYM12) 4837 4838*---------------------------------------------------------------------* 4839* LambdaP2 x LambdaH1 x Zeta 4840*---------------------------------------------------------------------* 4841 4842 DO ISYMAL = 1,NSYM 4843C 4844 ISYMBE = MULD2H(ISYMAL,ISALBE) 4845 ISYMB = MULD2H(ISYMAL,ISYXL2) 4846 ISYMK = MULD2H(ISYMBE,ISYXL1) 4847C 4848 IF (LWORK .LT. NBAS(ISYMAL)*NRHF(ISYMK) ) THEN 4849 CALL QUIT('Insufficient memory in CC_ZETADEN.') 4850 ENDIF 4851C 4852 KOFF1 = IGLMVI(ISYMAL,ISYMB) + 1 4853 KOFF2 = IT1AM(ISYMB,ISYMK) + 1 4854 NTOTAL = MAX(NBAS(ISYMAL),1) 4855 NTOTB = MAX(NVIR(ISYMB),1) 4856C 4857 CALL DGEMM('N','N',NBAS(ISYMAL),NRHF(ISYMK),NVIR(ISYMB), 4858 * ONE,XLAMP2(KOFF1),NTOTAL,ZETA1(KOFF2),NTOTB, 4859 * ZERO,WORK,NTOTAL) 4860C 4861 KOFF2 = IGLMRH(ISYMBE,ISYMK) + 1 4862 KOFF3 = IAODIS(ISYMAL,ISYMBE) + 1 4863 NTOTAL = MAX(NBAS(ISYMAL),1) 4864 NTOTBE = MAX(NBAS(ISYMBE),1) 4865C 4866 CALL DGEMM('N','T',NBAS(ISYMAL),NBAS(ISYMBE),NRHF(ISYMK), 4867 * ONE,WORK,NTOTAL,XLAMH1(KOFF2),NTOTBE, 4868 * ONE,ZDEN(KOFF3),NTOTAL) 4869 4870 END DO 4871 4872*---------------------------------------------------------------------* 4873* LambdaP1 x LambdaH2 x Zeta 4874*---------------------------------------------------------------------* 4875 4876 IF (IOPT .EQ. 2) THEN 4877C 4878 DO ISYMAL = 1,NSYM 4879C 4880 ISYMBE = MULD2H(ISYMAL,ISALBE) 4881 ISYMB = MULD2H(ISYMAL,ISYXL1) 4882 ISYMK = MULD2H(ISYMBE,ISYXL2) 4883C 4884 IF (LWORK .LT. NBAS(ISYMAL)*NRHF(ISYMK) ) THEN 4885 CALL QUIT('Insufficient memory in CC_ZETADEN.') 4886 ENDIF 4887C 4888 KOFF1 = IGLMVI(ISYMAL,ISYMB) + 1 4889 KOFF2 = IT1AM(ISYMB,ISYMK) + 1 4890 NTOTAL = MAX(NBAS(ISYMAL),1) 4891 NTOTB = MAX(NVIR(ISYMB),1) 4892C 4893 CALL DGEMM('N','N',NBAS(ISYMAL),NRHF(ISYMK),NVIR(ISYMB), 4894 * ONE,XLAMP1(KOFF1),NTOTAL,ZETA1(KOFF2),NTOTB, 4895 * ZERO,WORK,NTOTAL) 4896C 4897 KOFF2 = IGLMRH(ISYMBE,ISYMK) + 1 4898 KOFF3 = IAODIS(ISYMAL,ISYMBE) + 1 4899 NTOTAL = MAX(NBAS(ISYMAL),1) 4900 NTOTBE = MAX(NBAS(ISYMBE),1) 4901C 4902 CALL DGEMM('N','T',NBAS(ISYMAL),NBAS(ISYMBE),NRHF(ISYMK), 4903 * ONE,WORK,NTOTAL,XLAMH2(KOFF2),NTOTBE, 4904 * ONE,ZDEN(KOFF3),NTOTAL) 4905C 4906 END DO 4907C 4908 ENDIF 4909C 4910 CALL QEXIT('CC_ZETADEN') 4911C 4912 RETURN 4913 END 4914*=====================================================================* 4915C /* Deck cclt_chi */ 4916 SUBROUTINE CCLT_CHI(CTR1,XI,ISYCTR,XLAMDP,ISYLAM,CHI,IOPTX) 4917C 4918C Purpose: To calculate the Chi intermediate: 4919C 4920C Chi(alpha i) = sum_c XLAMDP(alpha c) CTR1(c i) 4921C - sum_k XLAMDP(alpha k) XI(k i) 4922C 4923C ISYCTR : symmetry of CTR1, XI 4924C ISYLAM : symmetry of XLAMDP 4925C 4926C if IOPTX <> 1 skip contribution from XI intermediate 4927C 4928C Christof Haettig, September 1998 4929C based on Asger Halkiers CCLT_CT1AO routine 4930C 4931C 4932#include "implicit.h" 4933#include "priunit.h" 4934#include "ccorb.h" 4935#include "ccsdsym.h" 4936C 4937 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 4938 DIMENSION CTR1(*),XLAMDP(*),CHI(*),XI(*) 4939C 4940 CALL QENTER('CCLT_CHI') 4941C 4942C--------------------------------------------- 4943C Half-transformation to AO-basis of CTR1. 4944C--------------------------------------------- 4945C 4946 ISYMAO = MULD2H(ISYCTR,ISYLAM) 4947 ISYMCJ = ISYCTR 4948C 4949 CALL DZERO(CHI,NGLMDT(ISYMAO)) 4950C 4951 DO ISYMAL = 1,NSYM 4952C 4953 ISYMC = MULD2H(ISYMAL,ISYLAM) 4954 ISYMJ = MULD2H(ISYMC,ISYMCJ) 4955C 4956 KOFF1 = IGLMVI(ISYMAL,ISYMC) + 1 4957 KOFF2 = IT1AM(ISYMC,ISYMJ) + 1 4958 KOFF3 = IGLMRH(ISYMAL,ISYMJ) + 1 4959C 4960 NTOTBA = MAX(NBAS(ISYMAL),1) 4961 NTOTVI = MAX(NVIR(ISYMC),1) 4962C 4963 CALL DGEMM('N','N',NBAS(ISYMAL),NRHF(ISYMJ),NVIR(ISYMC), 4964 * ONE,XLAMDP(KOFF1),NTOTBA,CTR1(KOFF2),NTOTVI, 4965 * ONE,CHI(KOFF3),NTOTBA) 4966C 4967 IF (IOPTX.EQ.1) THEN 4968C 4969 KOFF4 = IMATIJ(ISYMC,ISYMJ) + 1 4970 KOFF5 = IGLMRH(ISYMAL,ISYMC) + 1 4971C 4972 NTOTRH = MAX(NRHF(ISYMC),1) 4973C 4974 CALL DGEMM('N','N',NBAS(ISYMAL),NRHF(ISYMJ),NRHF(ISYMC), 4975 * -ONE,XLAMDP(KOFF5),NTOTBA,XI(KOFF4),NTOTRH, 4976 * ONE,CHI(KOFF3),NTOTBA) 4977C 4978 END IF 4979C 4980 END DO 4981C 4982 CALL QEXIT('CCLT_CHI') 4983C 4984 RETURN 4985 END 4986*=====================================================================* 4987*=====================================================================* 4988C /* Deck cc_bfif */ 4989 SUBROUTINE CC_BFIF(BFRHF,ISYRHF,XMGD,ISYMGD, 4990 & LUBFI,FNBFI,IADRBF,IADR,IDEL,WORK,LWORK) 4991*---------------------------------------------------------------------* 4992* 4993* Purpose: contract effective density with (**|k delta) integrals 4994* to the BF intermediate in the F matrix transformation 4995* (special version of the CC_BFI routine for F matrix) 4996* 4997* BFRHF : (**|k delta) integral distribution (symmetry ISYRHF) 4998* XMGD : effective density for BF term (symmetry ISYMGD) 4999* 5000* Written by Christof Haettig November 1998 5001*---------------------------------------------------------------------* 5002#if defined (IMPLICIT_NONE) 5003 IMPLICIT NONE 5004#else 5005# include "implicit.h" 5006#endif 5007#include "priunit.h" 5008#include "ccorb.h" 5009#include "maxorb.h" 5010#include "ccsdsym.h" 5011 5012 CHARACTER*(*) FNBFI 5013 INTEGER LWORK, ISYMGD, ISYRHF, IADR, IDEL, LUBFI, IADRBF(*) 5014 5015#if defined (SYS_CRAY) 5016 REAL BFRHF(*), XMGD(*), WORK(LWORK) 5017 REAL ZERO, ONE, HALF 5018#else 5019 DOUBLE PRECISION BFRHF(*), XMGD(*), WORK(LWORK) 5020 DOUBLE PRECISION ZERO, ONE, HALF 5021#endif 5022 PARAMETER(ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0) 5023C 5024 INTEGER NBFRHF(8), IBFRHF(8,8), ISYM, ICOUNT, ISYMAK, ISYBET 5025 INTEGER ISYMIJ, ISYMIJB, KOFF1, KOFF2, KOFF3, NBASB, NTOTAK 5026C 5027 CALL QENTER('CC_BFIF') 5028C 5029C -------------------------------------- 5030C precalculate symmetry array for BFRHF: 5031C -------------------------------------- 5032 DO ISYM = 1, NSYM 5033 ICOUNT = 0 5034 DO ISYMAK = 1, NSYM 5035 ISYBET = MULD2H(ISYMAK,ISYM) 5036 IBFRHF(ISYMAK,ISYBET) = ICOUNT 5037 ICOUNT = ICOUNT + NT1AO(ISYMAK)*NBAS(ISYBET) 5038 END DO 5039 NBFRHF(ISYM) = ICOUNT 5040 END DO 5041C 5042 ISYMIJB = MULD2H(ISYRHF,ISYMGD) 5043C 5044 IF (LWORK .LT. ND2IJG(ISYMIJB)) THEN 5045 WRITE (LUPRI,*) 'LWORK =',LWORK 5046 WRITE (LUPRI,*) 'need ',ND2IJG(ISYMIJB) 5047 CALL QUIT('Insufficient memory in CC_BFIF.') 5048 END IF 5049C 5050 DO ISYMIJ = 1, NSYM 5051C 5052 ISYMAK = MULD2H(ISYMGD,ISYMIJ) 5053 ISYBET = MULD2H(ISYMAK,ISYRHF) 5054C 5055 KOFF1 = IT2AOIJ(ISYMAK,ISYMIJ) + 1 5056 KOFF3 = IBFRHF(ISYMAK,ISYBET) + 1 5057 KOFF2 = ID2IJG(ISYMIJ,ISYBET) + 1 5058 NTOTAK = MAX(NT1AO(ISYMAK),1) 5059 NBASB = MAX(NBAS(ISYBET),1) 5060C 5061 CALL DGEMM('T','N',NBAS(ISYBET),NMATIJ(ISYMIJ),NT1AO(ISYMAK), 5062 & ONE, BFRHF(KOFF3),NTOTAK, XMGD(KOFF1),NTOTAK, 5063 & ZERO,WORK(KOFF2),NBASB) 5064C 5065 END DO 5066C 5067 IADRBF(IDEL) = IADR 5068C 5069 CALL PUTWA2(LUBFI,FNBFI,WORK,IADR,ND2IJG(ISYMIJB)) 5070C 5071 IADR = IADR + ND2IJG(ISYMIJB) 5072C 5073 CALL QEXIT('CC_BFIF') 5074C 5075 RETURN 5076 END 5077*=====================================================================* 5078*=====================================================================* 5079C /* Deck cc_bfifsort */ 5080 SUBROUTINE CC_BFIFSORT(BFI,ISYRES,LUBFI,FNBFI,IADRBF,WORK,LWORK) 5081*---------------------------------------------------------------------* 5082* 5083* Purpose: read and sort the BFZeta intermediate in F mat. transf. 5084* 5085* Written by Christof Haettig November 1998 5086*---------------------------------------------------------------------* 5087#if defined (IMPLICIT_NONE) 5088 IMPLICIT NONE 5089#else 5090# include "implicit.h" 5091#endif 5092#include "priunit.h" 5093#include "ccorb.h" 5094#include "maxorb.h" 5095#include "ccsdsym.h" 5096 5097 LOGICAL LOCDBG 5098 PARAMETER (LOCDBG = .FALSE.) 5099 5100 CHARACTER*(*) FNBFI 5101 INTEGER LWORK, ISYRES, IADR, LUBFI, IADRBF(*) 5102 5103#if defined (SYS_CRAY) 5104 REAL BFI(*), WORK(LWORK) 5105 REAL ONE, TWO 5106#else 5107 DOUBLE PRECISION BFI(*), WORK(LWORK) 5108 DOUBLE PRECISION ONE, TWO 5109#endif 5110 PARAMETER(TWO = 2.0D0, ONE = 1.0D0) 5111C 5112 INTEGER ISYDEL,ISYMIJB,ISYBET,ISYMI,ISYMJ,ISYMIJ,ISYBEJ,ISYDEI 5113 INTEGER NDI, NBJ, NIJ, NDIBJ, KOFF1, IDEL, NDB, INDEX 5114C 5115 INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 5116C 5117 CALL QENTER('CC_BFIFSORT') 5118C 5119 CALL DZERO(BFI,NT2AO(ISYRES)) 5120C 5121 DO ISYDEL = 1, NSYM 5122 5123 ISYMIJB = MULD2H(ISYDEL,ISYRES) 5124 5125 IF (LWORK .LT. ND2IJG(ISYMIJB) ) THEN 5126 CALL QUIT('Insufficient memory in CC_BFIFSORT.') 5127 END IF 5128 5129 DO D = 1, NBAS(ISYDEL) 5130 5131 IDEL = D + IBAS(ISYDEL) 5132 IADR = IADRBF(IDEL) 5133 5134 CALL GETWA2(LUBFI,FNBFI,WORK,IADR,ND2IJG(ISYMIJB)) 5135 5136 DO ISYMJ = 1, NSYM 5137 DO ISYMI = 1, NSYM 5138 5139 ISYMIJ = MULD2H(ISYMI,ISYMJ) 5140 ISYBET = MULD2H(ISYMIJB,ISYMIJ) 5141 ISYBEJ = MULD2H(ISYBET,ISYMJ) 5142 ISYDEI = MULD2H(ISYDEL,ISYMI) 5143 5144 DO J = 1, NRHF(ISYMJ) 5145 DO I = 1, NRHF(ISYMI) 5146 DO B = 1, NBAS(ISYBET) 5147 5148 NIJ = IMATIJ(ISYMI,ISYMJ) + NRHF(ISYMI)*(J-1) + I 5149 KOFF1 = ID2IJG(ISYMIJ,ISYBET) + NBAS(ISYBET)*(NIJ-1) + B 5150 5151 NBJ = IT1AO(ISYBET,ISYMJ) + NBAS(ISYBET)*(J-1) + B 5152 NDI = IT1AO(ISYDEL,ISYMI) + NBAS(ISYDEL)*(I-1) + D 5153 5154 IF (ISYDEI .EQ. ISYBEJ) THEN 5155 NDIBJ = IT2AO(ISYDEI,ISYBEJ) + INDEX(NDI,NBJ) 5156 IF (NDI.EQ.NBJ) WORK(KOFF1) = TWO * WORK(KOFF1) 5157 ELSE IF (ISYDEI .LT. ISYBEJ) THEN 5158 NDIBJ = IT2AO(ISYDEI,ISYBEJ)+NT1AO(ISYDEI)*(NBJ-1)+NDI 5159 ELSE IF (ISYDEI .GT. ISYBEJ) THEN 5160 NDIBJ = IT2AO(ISYDEI,ISYBEJ)+NT1AO(ISYBEJ)*(NDI-1)+NBJ 5161 END IF 5162 5163 BFI(NDIBJ) = BFI(NDIBJ) + WORK(KOFF1) 5164 5165 END DO 5166 END DO 5167 END DO 5168 5169 END DO 5170 END DO 5171 5172 END DO 5173 END DO 5174C 5175 IF (LOCDBG) THEN 5176 CALL AROUND("CC_BFIFSORT: The BF intermediate:") 5177 5178 WRITE (LUPRI,*) 'START OF BFI:',(BFI(I),I=1,10) 5179 5180 DO ISYMJ = 1, NSYM 5181 DO ISYMI = 1, NSYM 5182 DO ISYBET = 1, NSYM 5183 5184 ISYMIJ = MULD2H(ISYMI,ISYMJ) 5185 ISYBEJ = MULD2H(ISYMJ,ISYBET) 5186 ISYDEI = MULD2H(ISYBEJ,ISYRES) 5187 ISYDEL = MULD2H(ISYDEI,ISYMI) 5188 5189 IF (LWORK .LT. NBAS(ISYBET)*NBAS(ISYDEL) ) THEN 5190 CALL QUIT('Insufficient memory in CC_BFIFSORT.') 5191 END IF 5192 5193 DO J = 1, NRHF(ISYMJ) 5194 DO I = 1, NRHF(ISYMI) 5195 5196 DO B = 1, NBAS(ISYBET) 5197 DO D = 1, NBAS(ISYDEL) 5198 5199 NBJ = IT1AO(ISYBET,ISYMJ) + NBAS(ISYBET)*(J-1) + B 5200 NDI = IT1AO(ISYDEL,ISYMI) + NBAS(ISYDEL)*(I-1) + D 5201 NDB = NBAS(ISYDEL)*(B-1) + D 5202 5203 IF (ISYDEI .EQ. ISYBEJ) THEN 5204 NDIBJ = IT2AO(ISYDEI,ISYBEJ) + INDEX(NDI,NBJ) 5205 ELSEIF (ISYDEI .LT. ISYBEJ) THEN 5206 NDIBJ = IT2AO(ISYDEI,ISYBEJ)+NT1AO(ISYDEI)*(NBJ-1)+NDI 5207 ELSEIF (ISYDEI .GT. ISYBEJ) THEN 5208 NDIBJ = IT2AO(ISYDEI,ISYBEJ)+NT1AO(ISYBEJ)*(NDI-1)+NBJ 5209 ENDIF 5210 5211 WORK(NDB) = BFI(NDIBJ) 5212 5213 END DO 5214 END DO 5215 5216 WRITE (LUPRI,*) 'ISYMJ,ISYMI,ISYBET,ISYDEL,J,I:', 5217 & ISYMJ,ISYMI,ISYBET,ISYDEL,J,I 5218 5219 CALL OUTPUT(WORK,1,NBAS(ISYDEL),1,NBAS(ISYBET), 5220 & NBAS(ISYDEL),NBAS(ISYBET),1,LUPRI) 5221 5222 END DO 5223 END DO 5224 5225 END DO 5226 END DO 5227 END DO 5228 5229 END IF 5230C 5231 CALL QEXIT('CC_BFIFSORT') 5232C 5233 RETURN 5234 END 5235*=====================================================================* 5236*=====================================================================* 5237 SUBROUTINE CC_GIM(DSRHF,ISYRHF,XMGD,ISYMGD,GIM,IOPT,WORK,LWORK) 5238*---------------------------------------------------------------------* 5239* 5240* Purpose: Calculate G intermediate from effective BF density and 5241* (**|k del) integrals 5242* 5243* IOPT = 0 : compute left-hand-side GZeta intermediate 5244* 5245* IOPT = 1 : compute right-hand-side G intermediate, using 5246* 2 Cou - Exc combination of XMGD 5247* Note: this overwrites XMGD! 5248* 5249* symmetries: ISYRHF -- DSRHF 5250* ISYMGD -- XMGD 5251* 5252* Christof Haettig, November 1998 5253*---------------------------------------------------------------------* 5254#include "implicit.h" 5255#include "priunit.h" 5256#include "ccorb.h" 5257#include "ccsdsym.h" 5258 PARAMETER(ZERO=0.0D0,ONE=1.0D0) 5259 DIMENSION DSRHF(*),XMGD(*),GIM(*),WORK(LWORK) 5260C 5261 CALL QENTER('CC_GIM') 5262C 5263C ---------------------------------------------------------- 5264C if right-hand-side G requested replace XMGD by 2Cou - Exc: 5265C ---------------------------------------------------------- 5266 IF (IOPT.EQ.1) THEN 5267 CALL CC_MGDTCME(XMGD,ISYMGD,WORK,LWORK) 5268 END IF 5269C 5270 DO ISYMK = 1,NSYM 5271C 5272 ISYMAG = MULD2H(ISYMK,ISYRHF) 5273 ISYMGI = MULD2H(ISYMK,ISYMGD) 5274C 5275 IF (LWORK .LT. N2BST(ISYMAG)) THEN 5276 CALL QUIT('Insufficient space in CC_GIM') 5277 ENDIF 5278C 5279 DO K = 1,NRHF(ISYMK) 5280C 5281 KOFF1 = IDSRHF(ISYMAG,ISYMK) + NNBST(ISYMAG)*(K - 1) + 1 5282 CALL CCSD_SYMSQ(DSRHF(KOFF1),ISYMAG,WORK) 5283C 5284 DO ISYMI = 1,NSYM 5285C 5286 ISYMG = MULD2H(ISYMI,ISYMGI) 5287 ISYMA = MULD2H(ISYMG,ISYMAG) 5288C 5289 KOFF4 = IAODIS(ISYMA,ISYMG) + 1 5290 KOFF5 = IT2BGD(ISYMGI,ISYMK) + NT1AO(ISYMGI)*(K - 1) 5291 * + IT1AO(ISYMG,ISYMI) + 1 5292 KOFF6 = IT1AO(ISYMA,ISYMI) + 1 5293C 5294 NBASG = MAX(NBAS(ISYMG),1) 5295 NBASA = MAX(NBAS(ISYMA),1) 5296C 5297 CALL DGEMM('N','N',NBAS(ISYMA),NRHF(ISYMI),NBAS(ISYMG), 5298 * ONE,WORK(KOFF4),NBASA,XMGD(KOFF5),NBASG, 5299 * ONE,GIM(KOFF6), NBASA) 5300C 5301 END DO 5302 END DO 5303 END DO 5304 5305 CALL QEXIT('CC_GIM') 5306C 5307 RETURN 5308 END 5309*=====================================================================* 5310* END OF SUBROUTINE CC_GIM * 5311*=====================================================================* 5312*=====================================================================* 5313 SUBROUTINE CC_MGDTCME(XMGD,ISYMGD,WORK,LWORK) 5314*---------------------------------------------------------------------* 5315* 5316* Purpose: Calculate 2 Cou - Exc combination of two-index AO back- 5317* transformed amplitude batch XMGD(gam i;j)^del 5318* 5319* Christof Haettig, November 1998 5320*---------------------------------------------------------------------* 5321#include "implicit.h" 5322#include "priunit.h" 5323#include "ccorb.h" 5324#include "ccsdsym.h" 5325C 5326 PARAMETER(TWO=2.0D0,ONE=1.0D0) 5327 DIMENSION XMGD(*),WORK(LWORK) 5328C 5329 CALL QENTER('CC_MGDTCME') 5330C 5331 DO ISYMJ = 1,NSYM 5332C 5333 ISYMGI = MULD2H(ISYMJ,ISYMGD) 5334C 5335 DO J = 1,NRHF(ISYMJ) 5336C 5337 DO ISYMI = 1,ISYMJ 5338C 5339 ISYMG = MULD2H(ISYMI,ISYMGI) 5340 ISYMGJ = MULD2H(ISYMG,ISYMJ) 5341C 5342 KSCR1 = 1 5343 KSCR2 = KSCR1 + NBAS(ISYMG) 5344 KEND1 = KSCR2 + NBAS(ISYMG) 5345 LWRK1 = LWORK - KEND1 5346C 5347 IF (LWRK1 .LT. 0) THEN 5348 CALL QUIT('Insufficient space in CC_MGDTCME') 5349 ENDIF 5350C 5351 IF (ISYMI .EQ. ISYMJ) THEN 5352 NRHFI = J - 1 5353 ELSE 5354 NRHFI = NRHF(ISYMI) 5355 END IF 5356C 5357 DO I = 1,NRHFI 5358C 5359 NGIJ = IT2BGD(ISYMGI,ISYMJ) + NT1AO(ISYMGI)*(J-1) 5360 * + IT1AO(ISYMG,ISYMI) + NBAS(ISYMG)*(I-1) + 1 5361 NGJI = IT2BGD(ISYMGJ,ISYMI) + NT1AO(ISYMGJ)*(I-1) 5362 * + IT1AO(ISYMG,ISYMJ) + NBAS(ISYMG)*(J-1) + 1 5363C 5364 CALL DCOPY(NBAS(ISYMG),XMGD(NGIJ),1,WORK(KSCR1),1) 5365 CALL DCOPY(NBAS(ISYMG),XMGD(NGJI),1,WORK(KSCR2),1) 5366 CALL DSCAL(NBAS(ISYMG),TWO,XMGD(NGIJ),1) 5367 CALL DSCAL(NBAS(ISYMG),TWO,XMGD(NGJI),1) 5368 CALL DAXPY(NBAS(ISYMG),-ONE,WORK(KSCR2),1,XMGD(NGIJ),1) 5369 CALL DAXPY(NBAS(ISYMG),-ONE,WORK(KSCR1),1,XMGD(NGJI),1) 5370C 5371 END DO 5372 END DO 5373 END DO 5374 END DO 5375C 5376 CALL QEXIT('CC_MGDTCME') 5377C 5378 RETURN 5379 END 5380*=====================================================================* 5381* END OF SUBROUTINE CC_MGDTCME * 5382*=====================================================================* 5383*=====================================================================* 5384 SUBROUTINE CCSDT_F2B_SETUP(IFTRAN,IFDOTS,FCON,NFTRAN,MXVEC,LADD, 5385 & IB1TRAN,IB1DOTS,B1CON,NB1TRAN,MXB1VEC, 5386 & IB2TRAN,IB2DOTS,B2CON,NB2TRAN,MXB2VEC, 5387 & LISTL,LISTB,LISTC,MXB1TRAN,MXB2TRAN) 5388*---------------------------------------------------------------------* 5389* 5390* Purpose: transform from F matrix like loop structure to a B matrix 5391* like loop structure needed for CCSDT_FBMAT_CONT 5392* 5393* Christof Haettig, May 2002 5394*---------------------------------------------------------------------* 5395 IMPLICIT NONE 5396#include "priunit.h" 5397#include "ccorb.h" 5398#include "ccsdsym.h" 5399 5400 CHARACTER*(25) MSGDBG 5401 PARAMETER (MSGDBG = '[debug] CCSDT_F2B_SETUP> ') 5402 LOGICAL LOCDBG 5403 PARAMETER (LOCDBG = .FALSE.) 5404 5405*-------- 5406* input: 5407*-------- 5408 CHARACTER LISTL*(3), LISTC*(3), LISTB*(3) 5409 LOGICAL LADD 5410 INTEGER NFTRAN, MXVEC, MXB1VEC, MXB1TRAN, MXB2VEC, MXB2TRAN 5411 INTEGER IFTRAN(3,NFTRAN), IFDOTS(MXVEC,NFTRAN) 5412#if defined (SYS_CRAY) 5413 REAL FCON(MXVEC,NFTRAN) 5414#else 5415 DOUBLE PRECISION FCON(MXVEC,NFTRAN) 5416#endif 5417 5418*-------- 5419* output: 5420*-------- 5421 INTEGER NB1TRAN, NB2TRAN 5422 INTEGER IB1TRAN(3,MXB1TRAN), IB1DOTS(MXB1VEC,MXB1TRAN) 5423 INTEGER IB2TRAN(3,MXB2TRAN), IB2DOTS(MXB2VEC,MXB2TRAN) 5424#if defined (SYS_CRAY) 5425 REAL B1CON(MXB1VEC,MXB1TRAN), B2CON(MXB2VEC,MXB2TRAN) 5426#else 5427 DOUBLE PRECISION B1CON(MXB1VEC,MXB1TRAN), B2CON(MXB2VEC,MXB2TRAN) 5428#endif 5429 5430* local: 5431 INTEGER MXVAB1, MXVAB2, ITRAN, IZETA, ITAMPB, ITAMPC, IVEC, 5432 * JVEC, JTRAN 5433#if defined (SYS_CRAY) 5434 REAL CONT1, CONT2 5435#else 5436 DOUBLE PRECISION CONT1, CONT2 5437#endif 5438 5439 CALL QENTER('CCFCSU') 5440 5441 IF (.NOT. LADD) THEN 5442 NB1TRAN = 0 5443 NB2TRAN = 0 5444 END IF 5445 MXVAB1 = 0 5446 MXVAB2 = 0 5447 5448 DO ITRAN = 1, NFTRAN 5449 5450 IZETA = IFTRAN(1,ITRAN) 5451 IF (IZETA.EQ.0 .AND. LISTL(1:3).EQ.'L0 ') THEN 5452 IZETA = 1 5453 ELSE IF (IZETA.EQ.0) THEN 5454 CALL QUIT('Illegal index for left vector...') 5455 END IF 5456 5457 ITAMPB = IFTRAN(2,ITRAN) 5458 5459 IVEC = 1 5460 DO WHILE(IFDOTS(IVEC,ITRAN).NE.0 .AND. IVEC.LE.MXVEC) 5461 ITAMPC = IFDOTS(IVEC,ITRAN) 5462 5463 IF (LOCDBG) THEN 5464 WRITE(LUPRI,*) MSGDBG,'IZETA,ITAMPB,ITAMPC:', 5465 & IZETA,ITAMPB,ITAMPC 5466 END IF 5467 5468 ! R3=ITAMPB, R1=ITAMPC 5469 CALL CC_SETF12(IB1TRAN,IB1DOTS,MXB1TRAN,MXB1VEC, 5470 & ITAMPB,ITAMPC,IZETA,JTRAN,JVEC) 5471 NB1TRAN = MAX(NB1TRAN,JTRAN) 5472 MXVAB1 = MAX(MXVAB1,JVEC) 5473 CONT1 = B1CON(JVEC,JTRAN) 5474 5475 IF (LISTB(1:3).EQ.LISTC(1:3)) THEN 5476 ! R3=ITAMPC, R1=ITAMPB on the same list 5477 CALL CC_SETF12(IB1TRAN,IB1DOTS,MXB1TRAN,MXB1VEC, 5478 & ITAMPC,ITAMPB,IZETA,JTRAN,JVEC) 5479 NB1TRAN = MAX(NB1TRAN,JTRAN) 5480 MXVAB1 = MAX(MXVAB1,JVEC) 5481 CONT2 = B1CON(JVEC,JTRAN) 5482 ELSE 5483 ! R3=ITAMPC, R1=ITAMPB on a second list 5484 CALL CC_SETF12(IB2TRAN,IB2DOTS,MXB2TRAN,MXB2VEC, 5485 & ITAMPC,ITAMPB,IZETA,JTRAN,JVEC) 5486 NB2TRAN = MAX(NB2TRAN,JTRAN) 5487 MXVAB2 = MAX(MXVAB2,JVEC) 5488 CONT2 = B2CON(JVEC,JTRAN) 5489 END IF 5490 5491 IF (LADD) THEN 5492 FCON(IVEC,ITRAN) = FCON(IVEC,ITRAN) + CONT1 + CONT2 5493 5494 IF (LOCDBG) THEN 5495 WRITE(LUPRI,*)'ITRAN,IVEC:',ITRAN,IVEC 5496 WRITE(LUPRI,*)'FCON(before):',FCON(IVEC,ITRAN)-CONT1-CONT2 5497 WRITE(LUPRI,'(A,3I5)') 'IZETA,ITAMPB,ITAMPC:', 5498 & IZETA,ITAMPB,ITAMPC 5499 WRITE(LUPRI,*) '<L|[[H,RC_1],RB_3]|HF> = ',CONT1 5500 WRITE(LUPRI,*) '<L|[[H,RB_1],RC_3]|HF> = ',CONT2 5501 WRITE(LUPRI,*) 'FCON(after):',FCON(IVEC,ITRAN) 5502 END IF 5503 END IF 5504 5505 IVEC = IVEC + 1 5506 END DO 5507 5508 END DO 5509C 5510 IF (LOCDBG) THEN 5511 WRITE (LUPRI,*) 5512 WRITE (LUPRI,*) 'List of <L2|[[H,RC_1],RB_3]|HF>:' 5513 DO ITRAN = 1, NB1TRAN 5514 WRITE(LUPRI,'(A,2I5,5X,(12I5,20X))') MSGDBG, 5515 & (IB1TRAN(I,ITRAN),I=1,2),(IB1DOTS(I,ITRAN),I=1,MXVAB1) 5516 END DO 5517 WRITE (LUPRI,*) 5518 IF (LISTB(1:3).NE.LISTC(1:3)) THEN 5519 WRITE (LUPRI,*) 'List of <L2|[[H,RB_1],RC_3]|HF>:' 5520 DO ITRAN = 1, NB2TRAN 5521 WRITE(LUPRI,'(A,2I5,5X,(12I5,20X))') MSGDBG, 5522 & (IB2TRAN(I,ITRAN),I=1,2),(IB2DOTS(I,ITRAN),I=1,MXVAB2) 5523 END DO 5524 WRITE (LUPRI,*) 5525 END IF 5526 END IF 5527C 5528 CALL QEXIT('CCFCSU') 5529C 5530 RETURN 5531 END 5532*=====================================================================* 5533* END OF SUBROUTINE CCSDT_F2B_SETUP * 5534*=====================================================================* 5535 5536*---------------------------------------------------------------------* 5537c/* Deck CC_R12FMAT */ 5538*=====================================================================* 5539 SUBROUTINE CC_R12FMAT(THETA1, THETA2, THETAR12, 5540 & ISYRES, LISTL, IDLSTL, 5541 & ZETA1, ISYCTR, LISTR, IDLSTR, 5542 & ISYAMP, LAMDPA, LAMDHA, LAMP0, LAMH0, 5543 & XINT,LUMIM,FILMIM,LENM,IINT2,WORK,LWRK) 5544*---------------------------------------------------------------------* 5545* 5546* Purpose: calculate R12 contributions for F-matrix 5547* 5548* C. Neiss, C. Haettig autumn 2004 5549* Added finite fields, Chr. Neiss, June 2005 5550* Added CCSD(R12), Chr. Neiss, Nov. 2005 5551* 5552* THETA1 singles part of result vector 5553* THETA2 doubles part of result vector 5554* THETAR12 R12-doubles part of result vector 5555* ISYRES symmetry number of result vector 5556* LISTL type of L-vector; see CC_RDRSP.F 5557* IDLSTL index of L-vector; see CC_RDRSP.F 5558* ZETA1 singles part of Lagrangian multipliers 5559* ISYCTR symmetry of Lagrangian multipliers 5560* LISTR type of right response vector 5561* IDLSTR index of right response vector 5562* ISYAMP symmetry of right response vector 5563* LAMDPA Lambda^p-Matrix contracted with singles part of 5564* response vector 5565* LAMDHA Lambda^h-Matrix contracted with singles part of 5566* response vector 5567* LAMP0 Lambda^p-Matrix 5568* LAMH0 Lambda^h-Matrix 5569* 5570*=====================================================================* 5571 implicit none 5572#include "priunit.h" 5573#include "ccsdinp.h" 5574#include "ccsdsym.h" 5575#include "ccorb.h" 5576#include "dummy.h" 5577#include "r12int.h" 5578#include "ccr12int.h" 5579#include "ccfield.h" 5580 5581 LOGICAL LOCDBG 5582 PARAMETER (LOCDBG = .FALSE.) 5583 5584 LOGICAL LV, LVIJKL, LVAJKL 5585 5586 INTEGER LWRK, ISYRES, IDLSTL, ISYCTR, IDLSTR, ISYAMP 5587 INTEGER KEIM, ISYME, LWRK1, LWRK2, KEND1, KEND2, KEND3, LWRK3, 5588 & KR12SCR, NR12SCR, KVAJKL, 5589 & KZETA12, KCTR2, IOPT 5590 INTEGER ISYM1,ISYM2,LUNIT 5591 INTEGER LUMIM,LENM,IINT2,KMINT 5592 INTEGER KVABKL, KVINT 5593 INTEGER KT12AMP, KXINTTRI, KXINTSQ, KSCR, KFIELDAO, IFLD, IAN 5594 CHARACTER LISTL*3, LISTR*3, CDUMMY*8, MODEL*10 5595 CHARACTER*(*) FILMIM 5596 5597#if defined (SYS_CRAY) 5598 REAL WORK(LWRK), THETA1(*), THETA2(*), THETAR12(*), ZETA1(*), 5599 & LAMDPA(*), LAMDHA(*), LAMP0(*), LAMH0(*), XINT(*) 5600 REAL ZERO, ONE, TWO, HALF, DDOT 5601 REAL TIM1, TIM2, TIM3 5602#else 5603 DOUBLE PRECISION WORK(LWRK), THETA1(*), THETA2(*), THETAR12(*), 5604 & ZETA1(*), XINT(*), 5605 & LAMDPA(*), LAMDHA(*), LAMP0(*), LAMH0(*) 5606 DOUBLE PRECISION ZERO, ONE, TWO, HALF, DDOT 5607 DOUBLE PRECISION TIM1, TIM2, TIM3 5608#endif 5609 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0, HALF = 0.5D0) 5610 5611 5612 CALL QENTER('CC_R12FMAT') 5613 5614 IF (IANR12 .NE. 1) THEN 5615 CALL QUIT('Only available for CC-R12/ANSATZ 1') 5616 END IF 5617 IF (LOCDBG) THEN 5618 WRITE(LUPRI,*) "Theta1 and Theta2 before R12 section:" 5619 CALL CC_PRP(THETA1,THETA2,ISYRES,1,1) 5620 END IF 5621 5622 IF (CC2) THEN 5623C -------------------------------------------------------------- 5624C calculate G'-Term Singles contribution: 5625C do first E-intermediate calculation, then contract with 5626C singles Lagrangian multipliers and add to conventional term 5627C -------------------------------------------------------------- 5628C 5629C E-Intermediate: 5630C 5631 ISYME = ISYAMP 5632 IF (ISYRES .NE. MULD2H(ISYME,ISYCTR)) THEN 5633 WRITE(LUPRI,*) 'isyres, isyamp, isyctr: ',isyres, isyamp, isyctr 5634 CALL QUIT('Symmetry mismatch in cc_r12fmat!') 5635 END IF 5636 5637 KEIM = 1 5638 KEND1 = KEIM + NMATIJ(ISYME) 5639 LWRK1 = LWRK - KEND1 5640 IF (LWRK1 .LT. 0) THEN 5641 CALL QUIT('Insufficient work space in cc_r12fmat') 5642 END IF 5643 5644 CALL DZERO(WORK(KEIM),NMATIJ(ISYME)) 5645 CALL CCRHS_EINTP(WORK(KEIM),LAMP0,WORK(KEND1),LWRK1, 5646 & 2,ISYAMP,CDUMMY,IDUMMY,IDUMMY,LISTR,IDLSTR) 5647C 5648C contraction with singles Lagrange-multipliers: 5649C 5650 KEND2 = KEND1 + NMATAB(ISYME) 5651 LWRK2 = LWRK - KEND2 5652 IF (LWRK2 .LT. 0) THEN 5653 CALL QUIT('Insufficient work space in cc_r12fmat') 5654 END IF 5655 5656 CALL DZERO(WORK(KEND1),NMATAB(ISYME)) 5657 5658 CALL CCLR_E1C1(THETA1,ZETA1,WORK(KEND1),WORK(KEIM),WORK(KEND2), 5659 & LWRK2,ISYCTR,ISYME,'T') 5660C 5661 IF (LOCDBG) THEN 5662 WRITE(LUPRI,*) "ZETA1:" 5663 CALL CC_PRP(ZETA1,DUMMY,ISYCTR,1,0) 5664 WRITE(LUPRI,*) "E Intermediates:" 5665 CALL CC_PREI(WORK(KEND1),WORK(KEIM),ISYME,1) 5666 write (lupri,*) "Norm^2 (E Intermediates)", 5667 & ddot(nmatab(isyme),work(kend1),1,work(kend1),1), 5668 & ddot(nmatij(isyme),work(keim),1,work(keim),1) 5669 write(lupri,*) "Theta1 after G' term:" 5670 CALL CC_PRP(THETA1,THETA2,ISYRES,1,0) 5671 END IF 5672C 5673 END IF !CC2 5674C 5675C -------------------------------------------------------------- 5676C calculate F'-Term Singles contribution 5677C -------------------------------------------------------------- 5678C 5679C reallocate memory 5680C 5681 KVAJKL = 1 5682 KEND1 = KVAJKL + NVAJKL(ISYAMP) 5683 LWRK1 = LWRK - KEND1 5684 IF (LWRK1 .LT. 0) THEN 5685 CALL QUIT('Insufficient work space in cc_r12fmat') 5686 END IF 5687C 5688C calculate V_{kl}^{alpha jbar}-Intermediate 5689C 5690 IF (.NOT.USEVABKL) THEN 5691 IOPT = 1 5692 CALL CC_R12MKVAMKL0(WORK(KVAJKL),NVAJKL(ISYAMP),IOPT,LAMDHA, 5693 & ISYAMP,WORK(KEND1),LWRK1) 5694 CALL CC_MOFCONR12(LAMDHA,ISYAMP,IDUMMY,IDUMMY,IDUMMY, 5695 & IDUMMY,DUMMY,DUMMY,WORK(KVAJKL),IDUMMY, 5696 & .FALSE.,.TRUE.,.FALSE.,2, 5697 & TIM1,TIM2,TIM3, 5698 & IDUMMY,IDUMMY,IDUMMY,IDUMMY,IDUMMY,IDUMMY, 5699 & WORK(KEND1),LWRK1) 5700C 5701 ELSE 5702C 5703 KVABKL = KEND1 5704 KEND2 = KVABKL + NVABKL(1) 5705 LWRK2 = LWRK - KEND2 5706 IF (LWRK2 .LT. 0) THEN 5707 CALL QUIT('Insuff. work space for V^(albe)_kl in cc_r12fmat') 5708 END IF 5709C 5710 LV = .TRUE. 5711 LVAJKL = .FALSE. 5712 LVIJKL = .FALSE. 5713 CALL CC_R12MKVTF(WORK(KVABKL),WORK(KVAJKL),DUMMY, 5714 & LAMDHA,ISYAMP, 5715 & LV,LVIJKL,LVAJKL,CDUMMY,WORK(KEND2),LWRK2) 5716C 5717 END IF 5718C 5719 if (locdbg) then 5720 write(lupri,*) 'Norm^2 of VAJbarKL in CC_R12FMAT: ', 5721 & ddot(nvajkl(isyamp),work(kvajkl),1,work(kvajkl),1) 5722 write(lupri,*) 'VAJbarKL in CC_R12FMAT:' 5723 do isym2 = 1, nsym 5724 isym1 = muld2h(isyamp,isym2) 5725 write(lupri,*) 'Block isymaj, isymkl: ',isym1,isym2 5726 if ((nt1ao(isym1).eq.0).or.(nmatkl(isym2).eq.0)) then 5727 write(lupri,*) 'This block is empty' 5728 else 5729 call output(work(kvajkl+ivajkl(isym1,isym2)),1,nt1ao(isym1), 5730 & 1,nmatkl(isym2),nt1ao(isym1),nmatkl(isym2),1, 5731 & lupri) 5732 end if 5733 end do 5734 end if 5735C 5736C read in r12 doubles Lagrangian multipliers... 5737C 5738 KZETA12 = KEND1 5739 KEND1 = KZETA12 + NTR12SQ(ISYCTR) 5740 LWRK1 = LWRK - KEND1 5741 5742 CALL CC_R12GETCT(WORK(KZETA12),ISYCTR,2,2.0D0*BRASCL,.FALSE., 5743 & 'T',DUMMY,CDUMMY,DUMMY,LISTL,IDLSTL, 5744 & WORK(KEND1),LWRK1) 5745C 5746C ...and contract: 5747C 5748 CALL CCRHS_GP0(THETA1,ISYRES,WORK(KVAJKL),ISYAMP,LAMH0, 5749 & 1,WORK(KZETA12),ISYCTR,.FALSE.,DUMMY,LOCDBG, 5750 & WORK(KEND1),LWRK1) 5751C 5752C --------------------------------------------------------------- 5753C allocate memory for R12-doubles part of result vector (squared) 5754C --------------------------------------------------------------- 5755C 5756 KR12SCR = 1 5757 KEND1 = KR12SCR + NTR12SQ(ISYRES) 5758 LWRK1 = LWRK - KEND1 5759 IF (LWRK1 .LT. 0) THEN 5760 CALL QUIT('Insufficient work space in cc_r12fmat') 5761 END IF 5762 CALL DZERO(WORK(KR12SCR),NTR12SQ(ISYRES)) 5763C 5764C -------------------------------------------------------------- 5765C calculate G'-Term R12-doubles contribution: 5766C -------------------------------------------------------------- 5767C 5768 IF (LOCDBG) THEN 5769 WRITE (LUPRI,*) 'Call now CCLHTR_GP!' 5770 END IF 5771 CALL CCLHTR_GP(ZETA1,ISYCTR,LAMDPA,ISYAMP,WORK(KR12SCR),ISYRES, 5772 & WORK(KEND1),LWRK1) 5773 IF (LOCDBG) THEN 5774 WRITE(LUPRI,*) 'KR12SCR after CCLHTR_GP' 5775 call cc_prsqr12(WORK(KR12SCR),isyres,'T',1,.false.) 5776 END IF 5777C 5778C ------------------ 5779C add terms for CCSD 5780C ------------------ 5781C 5782 IF (CCSD .OR. CCSDT) THEN 5783C 5784C add B' term 5785C 5786C !read V_(alpha beta)^(kl) from disk 5787 KVABKL = KEND1 5788 KEND2 = KVABKL + NVABKL(1) 5789 LWRK2 = LWRK - KEND2 5790 IF (LWRK2 .LT. 0) THEN 5791 CALL QUIT('Insuff. work space for V^(albe)_kl in CC_R12FMAT') 5792 END IF 5793 lunit = -1 5794 call gpopen(lunit,fvabkl,'unknown',' ','unformatted', 5795 & idummy,.false.) 5796 read (lunit)(work(kvabkl+i-1),i=1,nvabkl(1)) 5797 call gpclose(lunit,'KEEP') 5798C 5799 IOPT = 2 5800 CALL CC_R12MKVIRT(WORK(KVABKL),LAMDPA,ISYAMP,LAMP0,1, 5801 & 'R12VCBDTKL',IOPT,WORK(KEND2),LWRK2) 5802C 5803 !read doubles Lagrangian multipliers 5804 KCTR2 = KEND1 5805 KEND2 = KCTR2 + NT2AM(ISYCTR) 5806 LWRK2 = LWRK - KEND2 5807 IF (LWRK2 .LT. 0) THEN 5808 CALL QUIT('Insuff. work space in CC_R12FMAT') 5809 END IF 5810 IOPT = 2 5811 CALL CC_RDRSP(LISTL,IDLSTL,ISYCTR,IOPT,MODEL,DUMMY,WORK(KCTR2)) 5812C 5813 !calculate contribution 5814 CALL CCRHS_BPP(WORK(KR12SCR),WORK(KCTR2),ISYCTR,.TRUE., 5815 & 'R12VCBDTKL',ISYAMP,WORK(KEND2),LWRK2) 5816C 5817C add A' term 5818C 5819 !read V-intermediate from disk: 5820 KVINT = KEND1 5821 KEND2 = KVINT + NTR12AM(1) 5822 LWRK2 = LWRK - KEND2 5823 IF (LWRK2 .LT. 0) THEN 5824 WRITE(LUPRI,*) 'Available:', LWRK, 'Needed:', KEND2 5825 CALL QUIT('Insufficient memory for V-intermediate '// 5826 & 'in CC_R12FMAT') 5827 ENDIF 5828 lunit = -1 5829 call gpopen(lunit,fccr12v,'old',' ','unformatted', 5830 & idummy,.FALSE.) 5831 9999 read(lunit) ian 5832 read(lunit) (work(kvint+i), i=0, ntr12am(1)-1) 5833 if (ian.ne.ianr12) goto 9999 5834 call gpclose(lunit,'KEEP') 5835C 5836 !read M-intermediate from disk: 5837 KMINT = KEND2 5838 KEND3 = KMINT + N3ORHF(ISYRES) 5839 LWRK3 = LWRK - KEND3 5840 IF (LWRK3 .LT. 0) THEN 5841 CALL QUIT('Insuff. work space in CC_R12FMAT') 5842 END IF 5843 CALL CC_RVEC(LUMIM,FILMIM,LENM,N3ORHF(ISYRES),IINT2, 5844 & WORK(KMINT)) 5845C 5846 CALL CCLHTR_AP(WORK(KR12SCR),WORK(KVINT),WORK(KMINT), 5847 & ISYRES,WORK(KEND3),LWRK3) 5848C 5849C add E' term 5850C 5851 CALL CCLHTR_EP(WORK(KR12SCR),WORK(KVINT),XINT, 5852 & ISYRES,WORK(KEND2),LWRK2) 5853 END IF 5854C 5855C resort result into a symmetry packed triangular matrix and 5856C apply the projection operator 0.5*P_(ij)^(kl): 5857C 5858 IOPT = 1 5859 CALL CCR12PCK2(THETAR12,ISYRES,.TRUE.,WORK(KR12SCR),'T', 5860 & IOPT) 5861 CALL CCLR_DIASCLR12(THETAR12,0.5D0*KETSCL,ISYRES) 5862 CALL DSCAL(NTR12AM(ISYRES),0.5D0,THETAR12,1) 5863 5864 if (locdbg) then 5865 write(lupri,*) 'theta_r12 after packing' 5866 call cc_prpr12(thetar12,isyres,1,.false.) 5867 end if 5868 5869C --------------------------------------------------- 5870C add finite field contributions: 5871C --------------------------------------------------- 5872 IF (NONHF) THEN 5873C allocate memory: 5874 KZETA12 = KEND1 5875 KT12AMP = KZETA12 + NTR12SQ(ISYCTR) 5876 KXINTTRI= KT12AMP + NTR12SQ(ISYAMP) 5877 KXINTSQ = KXINTTRI+ NR12R12P(1) 5878 KSCR = KXINTSQ + NR12R12SQ(1) 5879 KFIELDAO= KSCR + NTR12SQ(ISYRES) 5880 KEND2 = KFIELDAO+ N2BST(1) 5881 LWRK2 = LWRK - KEND2 5882 IF (LWRK2.LT.0) THEN 5883 CALL QUIT('Insufficient work space for R12 in CC_R12FMAT') 5884 END IF 5885 5886C initialize fields: 5887 CALL DZERO(WORK(KFIELDAO),N2BST(1)) 5888 CALL DZERO(WORK(KSCR),NTR12SQ(ISYRES)) 5889 CALL DZERO(WORK(KR12SCR),NTR12SQ(ISYRES)) 5890C 5891 IF (ISYMOP.NE.1) CALL QUIT('ISYMOP .NE. 1 not implemented') 5892 5893C sum up fields: 5894 DO IFLD = 1, NFIELD 5895 IF ( NHFFIELD(IFLD) ) THEN 5896 CALL CC_ONEP(WORK(KFIELDAO),WORK(KEND2),LWRK2, 5897 * EFIELD(IFLD),1,LFIELD(IFLD)) 5898 ELSE IF (.NOT. NHFFIELD(IFLD)) THEN 5899 CALL QUIT('CCR12 response can only handle '// 5900 & 'unrelaxed orbitals (w.r.t. the perturbation)') 5901 END IF 5902 END DO 5903 5904C read r12 Lagrangian multipliers: 5905 CALL CC_R12GETCT(WORK(KZETA12),ISYCTR,2,2.0D0*BRASCL,.FALSE., 5906 & 'N',DUMMY,CDUMMY,DUMMY,LISTL,IDLSTL, 5907 & WORK(KEND2),LWRK2) 5908 5909C read r12 response amplitudes: 5910 CALL CC_R12GETCT(WORK(KT12AMP),ISYAMP,2,KETSCL,.FALSE.,'N', 5911 & DUMMY,CDUMMY,DUMMY,LISTR,IDLSTR,WORK(KEND2),LWRK2) 5912 5913C read r12 overlap matrix: 5914 LUNIT = -1 5915 CALL GPOPEN(LUNIT,FCCR12X,'OLD',' ','UNFORMATTED', 5916 & IDUMMY,.FALSE.) 5917 REWIND(LUNIT) 5918 8888 READ(LUNIT) IAN 5919 READ(LUNIT) (WORK(KXINTTRI-1+I), I=1, NR12R12P(1)) 5920 IF (IAN.NE.IANR12) GOTO 8888 5921 CALL GPCLOSE(LUNIT,'KEEP') 5922 IOPT = 2 5923 CALL CCR12UNPCK2(WORK(KXINTTRI),1,WORK(KXINTSQ),'N',IOPT) 5924C 5925C calculate singles contribution: 5926C 5927 CALL CC_R12ETAA(THETA1,ISYRES,WORK(KZETA12),ISYCTR, 5928 & WORK(KT12AMP),ISYAMP,WORK(KXINTSQ), 5929 & WORK(KFIELDAO),1,LAMP0,LAMH0,.FALSE., 5930 & WORK(KEND2),LWRK2) 5931C 5932C calculate r12 doubles contribution: 5933C 5934 CALL CC_R12XI2A(WORK(KSCR),ISYRES,WORK(KZETA12),ISYCTR, 5935 & WORK(KFIELDAO),1,LAMP0,LAMDHA, 5936 & ISYAMP,'T',WORK(KEND2),LWRK2) 5937 CALL CC_R12XI2B(WORK(KR12SCR),'N',WORK(KXINTSQ),1,'N', 5938 & WORK(KSCR),ISYRES,'N',-ONE) 5939C 5940 IOPT = 1 5941 CALL CCR12PCK2(WORK(KSCR),ISYRES,.FALSE.,WORK(KR12SCR), 5942 & 'N',IOPT) 5943 CALL CCLR_DIASCLR12(WORK(KSCR),0.5D0*KETSCL,ISYRES) 5944 CALL DAXPY(NTR12AM(ISYRES),ONE,WORK(KSCR),1,THETAR12,1) 5945 5946 END IF 5947 5948 IF (LOCDBG) THEN 5949 WRITE(LUPRI,*) "Theta at end of CC_R12FMAT:" 5950 CALL CC_PRP(THETA1,THETA2,ISYRES,1,1) 5951 CALL CC_PRPR12(THETAR12,ISYRES,1,.TRUE.) 5952 END IF 5953 5954 CALL QEXIT('CC_R12FMAT') 5955 5956 RETURN 5957 END 5958 5959*=====================================================================* 5960* END OF SUBROUTINE CC_R12FMAT * 5961*=====================================================================* 5962