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! 18 SUBROUTINE CC3_GMAT(LISTL,IDLSTL,LISTB,IDLSTB, 19 & LISTC,IDLSTC,OMEGA1,OMEGA2,ISYRES,WORK,LWORK, 20 * LUCKJD,FNCKJD,LUTOC,FNTOC, 21 * LU3VI,FN3VI,LUDKBC3,FNDKBC3, 22 * LU3FOPX,FN3FOPX,LU3FOP2X,FN3FOP2X) 23C 24* Purpose: compute triples contribution to G matrix transformation 25* 26* (G T^B T^C)^eff_1,2 = (G T^B T^C)_1,2(CCSD) 27* + (G T^B T^C)_1,2(L_3) 28* 29C (G T^B T^C)_1,2(L_3) = 30C 31C + <L3|[H^BC,T^{0}_{2}],tau_{1}]|HF> 32C 33C + <L3|[H^B,T^{C}_{2}],tau_{1}]|HF> 34C 35C + <L3|[H^C,T^{B}_{2}],tau_{1}]|HF> 36C 37C + <L3|[H^BC,tau_{2}]|HF> 38C 39 IMPLICIT NONE 40C 41#include "priunit.h" 42#include "dummy.h" 43#include "ccsdsym.h" 44#include "inftap.h" 45#include "ccsdinp.h" 46#include "ccorb.h" 47#include "iratdef.h" 48#include "ccinftap.h" 49C 50 INTEGER IDLSTL,IDLSTB,IDLSTC,LWORK,LUCKJD,LUTOC,LU3VI,LUDKBC3 51 INTEGER LU3FOPX,LU3FOP2X 52 INTEGER LU4VBC,LU4VC,LU4VB,LU3SRTR,LUCKJDR,LUDELDR,LUDKBCR 53 INTEGER LUDKBCR4 54 INTEGER ISINT1,ISINT2,ISYMT1,ISYMT2,ISYML,ISYML1 55 INTEGER ISYML2,ISYMB,ISYMB1,ISYMB2,ISYMC1,ISYMC2 56 INTEGER ISINTB,ISINTC,ISYMBC,ISINTBC,ISYMIM,ISYRES 57 INTEGER KFOCKD,KFCKBA,KT1AM,KL1AM,KL2TP,KB1AM,KB2TP 58 INTEGER KC1AM,KC2TP,KT2TP,KLAMDP,KLAMDH,KEND0,LWRK0 59 INTEGER IOPT 60 INTEGER KBIOOOO,KBIOVVO,KBIOOVV,KCIOOOO,KCIOVVO,KCIOOVV 61 INTEGER KBCIOOOO,KBCIOVVO,KBCIOOVV 62 INTEGER ISYMC,ISYMK,KOFF1,KOFF2 63 INTEGER KTROC,KTROC0,KTROC1,KTROC2,KTROC3,KTROC4,KXIAJB 64 INTEGER KEND1,LWRK1,KINTOC,MAXOCC,KEND2,LWRK2,LENGTH 65 INTEGER ISYOPE,IOPTTCME,IOFF 66 INTEGER ISYMD,ISAIJ1,ISYCKB,ISCKB1,ISCKB2,ISCKBB,ISCKBC,ISCKBBC 67 INTEGER KTRVI7,KTRVI8,KTRVI2,KRMAT1,KTRVI0,KTRVI3,KTRVI4,KTRVI5 68 INTEGER KTRVI6,KVVVVB,KVVVVC,KVVVVBC,KEND3,LWRK3,KINTVI 69 INTEGER KEND4,LWRK4 70 INTEGER ISYALJ,ISAIJ2,ISYMBD,ISCKIJ,KSMAT,KQMAT,KDIAG 71 INTEGER KINDSQ,KINDEX,KTMAT,KRMAT2,LENSQ,ISYRES1,KOME1,KOME2 72 INTEGER LUFCK 73C 74COMMENT COMMENT 75COMMENT COMMENT 76 integer kt3am 77COMMENT COMMENT 78COMMENT COMMENT 79C 80C Functions : 81 INTEGER ILSTSYM 82C 83#if defined (SYS_CRAY) 84 REAL OMEGA1(*), OMEGA2(*) 85 REAL WORK(LWORK) 86 REAL XNORMVAL,DDOT 87 REAL HALF , ONE , TWO 88#else 89 DOUBLE PRECISION OMEGA1(*), OMEGA2(*) 90 DOUBLE PRECISION WORK(LWORK) 91 DOUBLE PRECISION XNORMVAL,DDOT 92 DOUBLE PRECISION HALF , ONE , TWO 93#endif 94C 95C 96 CHARACTER*(*) FNCKJD, FNTOC, FN3VI, FNDKBC3, FN3FOPX 97 CHARACTER*(*) FN3FOP2X 98 CHARACTER*1 CDUMMY 99 CHARACTER*3 LISTL, LISTB, LISTC 100 CHARACTER*10 MODEL 101 CHARACTER*13 FN4VC, FN4VB, FN4VBC, FN3SRTR, FNCKJDR 102 CHARACTER*13 FNDELDR, FNDKBCR, FNDKBCR4 103C 104 PARAMETER(HALF = 0.5D0, ONE = 1.0D0, TWO = 2.0D0) 105C 106 CALL QENTER('CC3_GMAT') 107C 108C---------------------------------------------------- 109C Initialise character strings and open files 110C---------------------------------------------------- 111C 112 CDUMMY = ' ' 113C 114 LU4VBC = -1 115 LU4VC = -1 116 LU4VB = -1 117 LU3SRTR = -1 118 LUCKJDR = -1 119 LUDELDR = -1 120 LUDKBCR = -1 121 LUDKBCR4 = -1 122C 123 FN4VBC = 'CC3_GMAT_TMP1' 124 FN4VC = 'CC3_GMAT_TMP2' 125 FN4VB = 'CC3_GMAT_TMP3' 126 FN3SRTR = 'CC3_GMAT_TMP4' 127 FNCKJDR = 'CC3_GMAT_TMP5' 128 FNDELDR = 'CC3_GMAT_TMP6' 129 FNDKBCR = 'CC3_GMAT_TMP7' 130 FNDKBCR4 = 'CC3_GMAT_TMP8' 131C 132 CALL WOPEN2(LU4VBC,FN4VBC,64,0) 133 CALL WOPEN2(LU4VC,FN4VC,64,0) 134 CALL WOPEN2(LU4VB,FN4VB,64,0) 135 CALL WOPEN2(LU3SRTR,FN3SRTR,64,0) 136 CALL WOPEN2(LUCKJDR,FNCKJDR,64,0) 137 CALL WOPEN2(LUDELDR,FNDELDR,64,0) 138 CALL WOPEN2(LUDKBCR,FNDKBCR,64,0) 139 CALL WOPEN2(LUDKBCR4,FNDKBCR4,64,0) 140C 141C------------------------------------------------------------- 142C Set symmetry flags. 143C 144C omega = int1*T2* ISYMB * ISYMC 145C isymres is symmetry of result(omega) 146C isint1 is symmetry of integrals in contraction.(int1) 147C the int1 integrals are one index transformed with B or C or BC 148C isint2 is symmetry of integrals in the triples equation.(int2) 149C isymim is symmetry of S and Q intermediates.(t2*int2) 150C (sym is for all index of S and Q (cbd,klj) 151C thus cklj=b*d*isymim) 152C------------------------------------------------------------- 153C 154 IPRCC = IPRINT 155 ISINT1 = ISYMOP 156 ISINT2 = ISYMOP 157 ISYMT1 = 1 158 ISYMT2 = 1 159 ISYML = ILSTSYM(LISTL,IDLSTL) 160 ISYML1 = ISYML 161 ISYML2 = ISYML 162 ISYMB = ILSTSYM(LISTB,IDLSTB) 163 ISYMB1 = ISYMB 164 ISYMB2 = ISYMB 165 ISYMC = ILSTSYM(LISTC,IDLSTC) 166 ISYMC1 = ISYMC 167 ISYMC2 = ISYMC 168 169 ISINTB = MULD2H(ISINT1,ISYMB) 170 ISINTC = MULD2H(ISINT1,ISYMC) 171 ISYMBC = MULD2H(ISYMB,ISYMC) 172 ISINTBC = MULD2H(ISINT1,ISYMBC) 173 ISYMIM = MULD2H(ISYML,ISYMOP) 174 ISYRES1 = MULD2H(ISYMIM,ISYMBC) 175 IF (ISYRES .NE. ISYRES1) THEN 176 WRITE(LUPRI,*) ' SYMMETRY MISMATCH IN CC3_GMAT' 177 WRITE(LUPRI,*) ' ISYRES =', ISYRES 178 WRITE(LUPRI,*) ' ISYRES1 =', ISYRES1 179 CALL QUIT('Symmetry mismatch IN CC3_GMAT') 180 END IF 181C 182C----------------------------------- 183C Work space allocation 184C----------------------------------- 185C 186 KOME1 = 1 187 KOME2 = KOME1 + NT1AM(ISYRES) 188 KFOCKD = KOME2 + NT2SQ(ISYRES) 189 KFCKBA = KFOCKD + NORBTS 190 KT1AM = KFCKBA + N2BST(ISYMOP) 191 KL1AM = KT1AM + NT1AM(ISYMT1) 192 KL2TP = KL1AM + NT1AM(ISYML1) 193 KB1AM = KL2TP + NT2SQ(ISYML2) 194 KB2TP = KB1AM + NT1AM(ISYMB1) 195 KC1AM = KB2TP + NT2SQ(ISYMB2) 196 KC2TP = KC1AM + NT1AM(ISYMC1) 197 KT2TP = KC2TP + NT2SQ(ISYMC2) 198 KLAMDP = KT2TP + NT2SQ(ISYMT2) 199 KLAMDH = KLAMDP + NLAMDT 200 KEND0 = KLAMDH + NLAMDT 201 LWRK0 = LWORK - KEND0 202C 203 CALL DZERO(WORK(KOME1),NT1AM(ISYRES)) 204 CALL DZERO(WORK(KOME2),NT2AM(ISYRES)) 205COMMENT COMMENT 206COMMENT COMMENT 207 if (.false.) then 208 kt3am = kend0 209 kend0 = kt3am + nt1amx*nt1amx*nt1amx 210 lwrk0 = lwork - kend0 211 call dzero(work(kt3am),nt1amx*nt1amx*nt1amx) 212 endif 213COMMENT COMMENT 214COMMENT COMMENT 215C 216 IF (LWRK0 .LT. 0) THEN 217 WRITE(LUPRI,*) 'Memory available : ',LWORK 218 WRITE(LUPRI,*) 'Memory needed : ',KEND0 219 CALL QUIT('Insufficient space in CC3_GMAT') 220 END IF 221C 222C--------------------------------------------------------- 223C Read the zero'th order amplitudes and calculate 224C the zero'th order Lambda matrices 225C--------------------------------------------------------- 226C 227 IOPT = 1 228 CALL CC_RDRSP('R0',0,1,IOPT,MODEL,WORK(KT1AM),DUMMY) 229C 230 CALL LAMMAT(WORK(KLAMDP),WORK(KLAMDH),WORK(KT1AM), 231 & WORK(KEND0),LWRK0) 232C 233C-------------------------------------- 234C Reorder the t2-amplitudes i T2TP. 235C-------------------------------------- 236C 237 IF (LWORK .LT. NT2SQ(1)) THEN 238 CALL QUIT('Not enough memory to construct T2TP in CC3_GMAT') 239 ENDIF 240C 241 IOPT = 2 242 CALL CC_RDRSP('R0',0,1,IOPT,MODEL, 243 * DUMMY,WORK(KT2TP)) 244 CALL CC_T2SQ(WORK(KT2TP),WORK(KEND0),ISYMT2) 245 CALL CC3_T2TP(WORK(KT2TP),WORK(KEND0),ISYMT2) 246C 247 IF (IPRINT .GT. 55) THEN 248 XNORMVAL = DDOT(NT2SQ(ISYMT2),WORK(KT2TP),1,WORK(KT2TP),1) 249 WRITE(LUPRI,*) 'Norm of T2TP ',XNORMVAL 250 ENDIF 251C 252C-------------------------------------- 253C Reorder the l2-amplitudes i L2TP. 254C-------------------------------------- 255C 256 IF (LWORK .LT. NT2SQ(ISYML2)) THEN 257 CALL QUIT('Not enough memory to construct L2TP in CC3_GMAT') 258 ENDIF 259C 260 IOPT = 3 261 CALL CC_RDRSP(LISTL,IDLSTL,ISYML2,IOPT,MODEL, 262 * WORK(KL1AM),WORK(KL2TP)) 263 CALL CC_T2SQ(WORK(KL2TP),WORK(KEND0),ISYML2) 264 CALL CC3_T2TP(WORK(KL2TP),WORK(KEND0),ISYML2) 265C 266 IF (IPRINT .GT. 55) THEN 267 XNORMVAL = DDOT(NT2SQ(ISYML2),WORK(KL2TP),1,WORK(KL2TP),1) 268 WRITE(LUPRI,*) 'Norm of L2TP ',XNORMVAL 269 ENDIF 270C 271C-------------------------------------- 272C Reorder the B2-amplitudes i B2TP. 273C-------------------------------------- 274C 275 IF (LWORK .LT. NT2SQ(ISYMB2)) THEN 276 CALL QUIT('Not enough memory to construct B2TP in CC3_GMAT') 277 ENDIF 278C 279 IOPT = 3 280 CALL CC_RDRSP(LISTB,IDLSTB,ISYMB2,IOPT,MODEL, 281 & WORK(KB1AM),WORK(KB2TP)) 282 CALL CCLR_DIASCL(WORK(KB2TP),TWO,ISYMB2) 283 CALL CC_T2SQ(WORK(KB2TP),WORK(KEND0),ISYMB2) 284 CALL CC3_T2TP(WORK(KB2TP),WORK(KEND0),ISYMB2) 285C 286 IF (IPRINT .GT. 55) THEN 287 XNORMVAL = DDOT(NT2SQ(ISYMB2),WORK(KB2TP),1,WORK(KB2TP),1) 288 WRITE(LUPRI,*) 'Norm of B2TP ',XNORMVAL 289 ENDIF 290C 291C-------------------------------------- 292C Reorder the C2-amplitudes i C2TP. 293C-------------------------------------- 294C 295 IF (LWORK .LT. NT2SQ(ISYMC2)) THEN 296 CALL QUIT('Not enough memory to construct B2TP in CC3_GMAT') 297 ENDIF 298C 299 IOPT = 3 300 CALL CC_RDRSP(LISTC,IDLSTC,ISYMC2,IOPT,MODEL, 301 & WORK(KC1AM),WORK(KC2TP)) 302 CALL CCLR_DIASCL(WORK(KC2TP),TWO,ISYMC2) 303 CALL CC_T2SQ(WORK(KC2TP),WORK(KEND0),ISYMC2) 304 CALL CC3_T2TP(WORK(KC2TP),WORK(KEND0),ISYMC2) 305C 306 IF (IPRINT .GT. 55) THEN 307 XNORMVAL = DDOT(NT2SQ(ISYMC2),WORK(KC2TP),1,WORK(KC2TP),1) 308 WRITE(LUPRI,*) 'Norm of C2TP ',XNORMVAL 309 ENDIF 310C 311 312C 313C------------------------------------------------------ 314C Calculate the (ck|de)-BCtilde and (ck|lm)-BCtilde 315C------------------------------------------------------ 316C 317 CALL CC3_3BARINT(ISYMB1,LISTB,IDLSTB,ISYMC1,LISTC,IDLSTC, 318 * IDUMMY,CDUMMY,IDUMMY,.FALSE., 319 * WORK(KLAMDP),WORK(KLAMDH),WORK(KEND0),LWRK0, 320 * LU3SRTR,FN3SRTR,LUCKJDR,FNCKJDR) 321C 322 CALL CC3_SORT1(WORK(KEND0),LWRK0,2,ISYMBC,LU3SRTR,FN3SRTR, 323 * LUDELDR,FNDELDR,IDUMMY,CDUMMY,IDUMMY,CDUMMY, 324 * IDUMMY,CDUMMY) 325C 326 CALL CC3_SINT(WORK(KLAMDH),WORK(KEND0),LWRK0,ISYMBC, 327 * LUDELDR,FNDELDR,LUDKBCR,FNDKBCR) 328C 329 CALL CC3_TCME(WORK(KLAMDP),ISYMBC,WORK(KEND0),LWRK0, 330 * IDUMMY,CDUMMY,LUDKBCR,FNDKBCR, 331 * IDUMMY,CDUMMY,IDUMMY,CDUMMY, 332 * IDUMMY,CDUMMY,LUDKBCR4,FNDKBCR4,2) 333C 334C--------------------------------------------------------------- 335C Calculate the g^B integrals for 336C OOOO, OVVO and OOVV integrals. VVVV is stored on file 337C--------------------------------------------------------------- 338C 339 KBIOOOO = KEND0 340 KBIOVVO = KBIOOOO + N3ORHF(ISINTB) 341 KBIOOVV = KBIOVVO + NT2SQ(ISINTB) 342 KEND0 = KBIOOVV + NT2SQ(ISINTB) 343 LWRK0 = LWORK - KEND0 344C 345 IF (LWRK0 .LT. 0) THEN 346 CALL QUIT('Out of memory in CC3_GMAT (g^B[2o2v] kind)') 347 ENDIF 348C 349 CALL DZERO(WORK(KBIOOOO),N3ORHF(ISINTB)) 350 CALL DZERO(WORK(KBIOVVO),NT2SQ(ISINTB)) 351 CALL DZERO(WORK(KBIOOVV),NT2SQ(ISINTB)) 352C 353 CALL CC3_INT(WORK(KBIOOOO),WORK(KBIOOVV),WORK(KBIOVVO), 354 * ISINTB,LU4VB,FN4VB, 355 * .TRUE.,LISTB,IDLSTB,ISYMB1, 356 * .FALSE.,'DUMMY',IDUMMY,IDUMMY, 357 * WORK(KEND0),LWRK0) 358C 359C--------------------------------------------------------------- 360C Calculate the g^C integrals for 361C OOOO, OVVO and OOVV integrals. VVVV is stored on file 362C--------------------------------------------------------------- 363C 364 KCIOOOO = KEND0 365 KCIOVVO = KCIOOOO + N3ORHF(ISINTC) 366 KCIOOVV = KCIOVVO + NT2SQ(ISINTC) 367 KEND0 = KCIOOVV + NT2SQ(ISINTC) 368 LWRK0 = LWORK - KEND0 369C 370 IF (LWRK0 .LT. 0) THEN 371 CALL QUIT('Out of memory in CC3_GMAT (g^C[2o2v] kind)') 372 ENDIF 373C 374 CALL DZERO(WORK(KCIOOOO),N3ORHF(ISINTC)) 375 CALL DZERO(WORK(KCIOVVO),NT2SQ(ISINTC)) 376 CALL DZERO(WORK(KCIOOVV),NT2SQ(ISINTC)) 377C 378 CALL CC3_INT(WORK(KCIOOOO),WORK(KCIOOVV),WORK(KCIOVVO), 379 * ISINTC,LU4VC,FN4VC, 380 * .TRUE.,LISTC,IDLSTC,ISYMC1, 381 * .FALSE.,'DUMMY',IDUMMY,IDUMMY, 382 * WORK(KEND0),LWRK0) 383C 384C--------------------------------------------------------------- 385C Calculate the g^BC integrals for 386C OOOO, OVVO and OOVV integrals. VVVV is stored on file 387C--------------------------------------------------------------- 388C 389 KBCIOOOO = KEND0 390 KBCIOVVO = KBCIOOOO + N3ORHF(ISINTBC) 391 KBCIOOVV = KBCIOVVO + NT2SQ(ISINTBC) 392 KEND0 = KBCIOOVV + NT2SQ(ISINTBC) 393 LWRK0 = LWORK - KEND0 394C 395 IF (LWRK0 .LT. 0) THEN 396 CALL QUIT('Out of memory in CC3_GMAT (g^B[2o2v] kind)') 397 ENDIF 398C 399 CALL DZERO(WORK(KBCIOOOO),N3ORHF(ISINTBC)) 400 CALL DZERO(WORK(KBCIOVVO),NT2SQ(ISINTBC)) 401 CALL DZERO(WORK(KBCIOOVV),NT2SQ(ISINTBC)) 402C 403 CALL CC3_INT(WORK(KBCIOOOO),WORK(KBCIOOVV),WORK(KBCIOVVO), 404 * ISINTBC,LU4VBC,FN4VBC, 405 * .TRUE.,LISTB,IDLSTB,ISYMB1, 406 * .TRUE.,LISTC,IDLSTC,ISYMC1, 407 * WORK(KEND0),LWRK0) 408C 409C--------------------------------------------------------- 410C Read canonical orbital energies and MO coefficients. 411C--------------------------------------------------------- 412C 413 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY, 414 & .FALSE.) 415 REWIND LUSIFC 416C 417 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 418 READ (LUSIFC) 419 READ (LUSIFC) (WORK(KFOCKD+I-1), I=1,NORBTS) 420C 421 CALL GPCLOSE(LUSIFC,'KEEP') 422C 423C--------------------------------------------- 424C Delete frozen orbitals in Fock diagonal. 425C--------------------------------------------- 426C 427 IF (FROIMP .OR. FROEXP) 428 * CALL CCSD_DELFRO(WORK(KFOCKD),WORK(KEND0),LWRK0) 429C 430C----------------------------------------------------- 431C Construct the transformed Fock matrix 432C----------------------------------------------------- 433C 434 LUFCK = -1 435C This AO Fock matrix is constructed from the T1 transformed density 436 CALL GPOPEN(LUFCK,'CC_FCKH','UNKNOWN',' ','UNFORMATTED', 437 * IDUMMY,.FALSE.) 438C This AO Fock matrix is constructed from the CMO transformed density 439C CALL GPOPEN(LUFCK,'CC_FCKREF','UNKNOWN',' ','UNFORMATTED', 440C * IDUMMY,.FALSE.) 441 REWIND(LUFCK) 442 READ(LUFCK)(WORK(KFCKBA + I-1),I = 1,N2BST(ISYMOP)) 443 CALL GPCLOSE(LUFCK,'KEEP' ) 444C 445 IF (IPRINT .GT. 140) THEN 446 CALL AROUND( 'Usual Fock AO matrix' ) 447 CALL CC_PRFCKAO(WORK(KFCKBA),ISYMOP) 448 ENDIF 449C 450 CALL CC_FCKMO(WORK(KFCKBA),WORK(KLAMDP),WORK(KLAMDH), 451 * WORK(KEND0),LWRK0,1,1,1) 452C 453 IF (IPRINT .GT. 50) THEN 454 CALL AROUND( 'In CC3_GMAT: Triples Fock MO matrix' ) 455 CALL CC_PRFCKMO(WORK(KFCKBA),ISYMOP) 456 ENDIF 457C 458C Sort the fock matrix 459C 460C 461 CALL DCOPY(N2BST(ISINT1),WORK(KFCKBA),1,WORK(KEND0),1) 462C 463 DO ISYMC = 1,NSYM 464C 465 ISYMK = MULD2H(ISYMC,ISINT1) 466C 467 DO K = 1,NRHF(ISYMK) 468C 469 DO C = 1,NVIR(ISYMC) 470C 471 KOFF1 = KEND0 + IFCVIR(ISYMK,ISYMC) + 472 * NORB(ISYMK)*(C - 1) + K - 1 473 KOFF2 = KFCKBA + IT1AM(ISYMC,ISYMK) 474 * + NVIR(ISYMC)*(K - 1) + C - 1 475C 476 WORK(KOFF2) = WORK(KOFF1) 477C 478 ENDDO 479 ENDDO 480 ENDDO 481C 482 IF (IPRINT .GT. 50) THEN 483 CALL AROUND('In CC3_GMAT: Triples Fock MO matrix (sort)') 484 CALL CC_PRFCKMO(WORK(KFCKBA),ISYMOP) 485 ENDIF 486C 487C----------------------------- 488C Read occupied integrals. 489C----------------------------- 490C 491C Memory allocation. 492C 493 KTROC = KEND0 494 KTROC0 = KTROC + NTRAOC(ISINT2) 495 KTROC1 = KTROC0 + NTRAOC(ISINT2) 496 KTROC2 = KTROC1 + NTRAOC(ISINT2) 497 KTROC3 = KTROC2 + NTRAOC(ISINT2) 498 KTROC4 = KTROC3 + NTRAOC(ISINTBC) 499 KXIAJB = KTROC4 + NTRAOC(ISINTBC) 500 KEND1 = KXIAJB + NT2AM(ISYMOP) 501 LWRK1 = LWORK - KEND1 502C 503 KINTOC = KEND1 504 MAXOCC = MAX(NTOTOC(ISINT2),NTOTOC(ISINTBC)) 505 KEND2 = KINTOC + MAX(NTOTOC(ISYMOP),MAXOCC) 506 LWRK2 = LWORK - KEND2 507C 508 IF (LWRK2 .LT. 0) THEN 509 WRITE(LUPRI,*) 'Memory available : ',LWORK 510 WRITE(LUPRI,*) 'Memory needed : ',KEND2 511 CALL QUIT('Insufficient space in CC3_LHTR_L3') 512 END IF 513C 514C------------------------ 515C Construct L(ia,jb). 516C------------------------ 517C 518 LENGTH = IRAT*NT2AM(ISYMOP) 519C 520 REWIND(LUIAJB) 521 CALL READI(LUIAJB,LENGTH,WORK(KXIAJB)) 522C 523 ISYOPE = ISYMOP 524 IOPTTCME = 1 525 CALL CCSD_TCMEPK(WORK(KXIAJB),1.0D0,ISYOPE,IOPTTCME) 526C 527 IF ( IPRINT .GT. 55) THEN 528 XNORMVAL = DDOT(NT2AM(ISYMOP),WORK(KXIAJB),1, 529 * WORK(KXIAJB),1) 530 WRITE(LUPRI,*) 'Norm of IAJB ',XNORMVAL 531 ENDIF 532C 533C-------------------------------------------------------------------------- 534C Read in BC-transformed integrals used in contractions and transform. 535C-------------------------------------------------------------------------- 536C 537 IOFF = 1 538 IF (NTOTOC(ISINTBC) .GT. 0) THEN 539 CALL GETWA2(LUCKJDR,FNCKJDR,WORK(KINTOC),IOFF,NTOTOC(ISINTBC)) 540 ENDIF 541C 542 CALL CC3_TROCC(WORK(KINTOC),WORK(KTROC3),WORK(KLAMDP), 543 * WORK(KEND2),LWRK2,ISINTBC) 544C 545 CALL CCFOP_SORT(WORK(KTROC3),WORK(KTROC4),ISINTBC,1) 546C 547 CALL CC3_LSORT1(WORK(KTROC3),ISINTBC,WORK(KEND2),LWRK2,5) 548C 549C------------------------------------------------------------ 550C Read in integrals used in contractions and transform. 551C------------------------------------------------------------ 552C 553 IOFF = 1 554 IF (NTOTOC(ISINT2) .GT. 0) THEN 555 CALL GETWA2(LUCKJD,FNCKJD,WORK(KINTOC),IOFF,NTOTOC(ISINT2)) 556 ENDIF 557C 558 CALL CC3_TROCC(WORK(KINTOC),WORK(KTROC),WORK(KLAMDP), 559 * WORK(KEND2),LWRK2,ISINT2) 560C 561 CALL CCFOP_SORT(WORK(KTROC),WORK(KTROC1),ISINT2,1) 562C 563 CALL CC3_LSORT1(WORK(KTROC),ISINT2,WORK(KEND2),LWRK2,5) 564C 565C------------------------------------------------------------ 566C Read in integrals used in L3 and transform. 567C------------------------------------------------------------ 568C 569 IOFF = 1 570 IF (NTOTOC(ISINT2) .GT. 0) THEN 571 CALL GETWA2(LUTOC,FNTOC,WORK(KINTOC),IOFF,NTOTOC(ISINT2)) 572 ENDIF 573C 574 CALL CC3_TROCC(WORK(KINTOC),WORK(KTROC0),WORK(KLAMDH), 575 * WORK(KEND2),LWRK2,ISINT2) 576C 577C----------------------------------------------------------- 578C Construct 2*C-E of the integrals. 579C----------------------------------------------------------- 580C 581 CALL CCSDT_TCMEOCC(WORK(KTROC0),WORK(KTROC2),ISINT2) 582C 583C------------------------------- 584C Write out norms of arrays. 585C------------------------------- 586C 587 IF (IPRINT .GT. 55) THEN 588 XNORMVAL = DDOT(NTOTOC(ISINT2),WORK(KINTOC),1, 589 * WORK(KINTOC),1) 590 WRITE(LUPRI,*) 'Norm of CKJDEL-INT ',XNORMVAL 591 ENDIF 592C 593 IF (IPRINT .GT. 55) THEN 594 XNORMVAL = DDOT(NTRAOC(ISINT2),WORK(KTROC0),1, 595 * WORK(KTROC0),1) 596 WRITE(LUPRI,*) 'Norm of TROC0 CC3_GMAT : ',XNORMVAL 597 ENDIF 598C 599 IF (IPRINT .GT. 55) THEN 600 XNORMVAL = DDOT(NTRAOC(ISINT2),WORK(KTROC2),1, 601 * WORK(KTROC2),1) 602 WRITE(LUPRI,*) 'Norm of TROC2 CC3_GMAT : ',XNORMVAL 603 ENDIF 604C 605C---------------------------- 606C General loop structure. 607C---------------------------- 608C 609 DO ISYMD = 1,NSYM 610C 611 ISAIJ1 = MULD2H(ISYMD,ISYRES) 612 ISYCKB = MULD2H(ISYMD,ISYMOP) 613 ISCKB1 = MULD2H(ISINT1,ISYMD) 614 ISCKB2 = MULD2H(ISINT2,ISYMD) 615 616 ISCKBB = MULD2H(ISINTB,ISYMD) 617 ISCKBC = MULD2H(ISINTC,ISYMD) 618 ISCKBBC = MULD2H(ISINTBC,ISYMD) 619C 620 IF (IPRINT .GT. 55) THEN 621C 622 WRITE(LUPRI,*) 'In CC3_GMAT: ISAIJ1 :',ISAIJ1 623 WRITE(LUPRI,*) 'In CC3_GMAT: ISYCKB :',ISYCKB 624 WRITE(LUPRI,*) 'In CC3_GMAT: ISCKB2 :',ISCKB2 625 626 WRITE(LUPRI,*) 'In CC3_GMAT: ISCKBB :',ISCKBB 627 WRITE(LUPRI,*) 'In CC3_GMAT: ISCKBC :',ISCKBC 628 WRITE(LUPRI,*) 'In CC3_GMAT: ISCKBBC :',ISCKBBC 629C 630 ENDIF 631C 632C-------------------------- 633C Memory allocation. 634C-------------------------- 635C 636 KTRVI7 = KEND1 637 KTRVI8 = KTRVI7 + NCKATR(ISCKBBC) 638 KTRVI2 = KTRVI8 + NCKATR(ISCKBBC) 639 KRMAT1 = KTRVI2 + NCKATR(ISCKB2) 640 KEND2 = KRMAT1 + NCKI(ISAIJ1) 641 LWRK2 = LWORK - KEND2 642C 643 KTRVI0 = KEND2 644 KTRVI3 = KTRVI0 + NCKATR(ISCKB2) 645 KTRVI4 = KTRVI3 + NCKATR(ISCKB2) 646 KTRVI5 = KTRVI4 + NCKATR(ISCKB2) 647 KTRVI6 = KTRVI5 + NCKATR(ISCKB2) 648 KVVVVB = KTRVI6 + NCKATR(ISCKB2) 649 KVVVVC = KVVVVB + NMAABC(ISCKBB) 650 KVVVVBC = KVVVVC + NMAABC(ISCKBC) 651 KEND3 = KVVVVBC + NMAABC(ISCKBBC) 652 LWRK3 = LWORK - KEND3 653C 654 KINTVI = KEND3 655 KEND4 = KINTVI + MAX(NCKA(ISYMD),NCKA(ISCKB2)) 656 LWRK4 = LWORK - KEND4 657C 658 IF (LWRK4 .LT. 0) THEN 659 WRITE(LUPRI,*) 'Memory available : ',LWORK 660 WRITE(LUPRI,*) 'Memory needed : ',KEND4 661 CALL QUIT('Insufficient space in CC3_LHTR_L3') 662 END IF 663C 664C--------------------- 665C Sum over D 666C--------------------- 667C 668 DO D = 1,NVIR(ISYMD) 669C 670C------------------------------------ 671C Initialize the R1 matrix. 672C------------------------------------ 673C 674 CALL DZERO(WORK(KRMAT1),NCKI(ISAIJ1)) 675C 676C-------------------------------------------- 677C Read in g^B_{vvvv} for a given D 678C-------------------------------------------- 679C 680 IF (NMAABC(ISCKBB) .GT. 0) THEN 681 IOFF = I3VVIR(ISCKBB,ISYMD) 682 * + NMAABC(ISCKBB)*(D-1) 683 * + 1 684 CALL GETWA2(LU4VB,FN4VB,WORK(KVVVVB),IOFF,NMAABC(ISCKBB)) 685 ENDIF 686C 687C-------------------------------------------- 688C Read in g^B_{vvvv} for a given D 689C-------------------------------------------- 690C 691 IF (NMAABC(ISCKBC) .GT. 0) THEN 692 IOFF = I3VVIR(ISCKBC,ISYMD) 693 * + NMAABC(ISCKBC)*(D-1) 694 * + 1 695 CALL GETWA2(LU4VC,FN4VC,WORK(KVVVVC),IOFF,NMAABC(ISCKBC)) 696 ENDIF 697C 698C 699C-------------------------------------------- 700C Read in g^B_{vvvv} for a given D 701C-------------------------------------------- 702C 703 IF (NMAABC(ISCKBBC) .GT. 0) THEN 704 IOFF = I3VVIR(ISCKBBC,ISYMD) 705 * + NMAABC(ISCKBBC)*(D-1) 706 * + 1 707 CALL GETWA2(LU4VBC,FN4VBC,WORK(KVVVVBC),IOFF, 708 * NMAABC(ISCKBBC)) 709 ENDIF 710C 711C-------------------------------------------------------------- 712C Read B-transformed integrals used in contraction 713C-------------------------------------------------------------- 714C 715 IF (NCKATR(ISCKBBC) .GT. 0) THEN 716 IOFF = ICKBD(ISCKBBC,ISYMD) + NCKATR(ISCKBBC)*(D - 1) + 1 717 CALL GETWA2(LUDKBCR,FNDKBCR,WORK(KTRVI7),IOFF, 718 & NCKATR(ISCKBBC)) 719 ENDIF 720C 721 IF (LWRK4 .LT. NCKATR(ISCKBBC)) THEN 722 CALL QUIT('Insufficient space for allocation in '// 723 & 'CC3_GMAT (TRVI7)') 724 END IF 725C 726 CALL CCSDT_SRVIR3(WORK(KTRVI7),WORK(KEND4),ISYMD,D,ISINTBC) 727C 728 IF (NCKATR(ISCKBBC) .GT. 0) THEN 729 IOFF = ICKBD(ISCKBBC,ISYMD) + NCKATR(ISCKBBC)*(D - 1) + 1 730 CALL GETWA2(LUDKBCR4,FNDKBCR4,WORK(KTRVI8),IOFF, 731 & NCKATR(ISCKBBC)) 732 ENDIF 733C 734 IF (LWRK4 .LT. NCKATR(ISCKBBC)) THEN 735 CALL QUIT('Insufficient space for allocation in '// 736 & 'CC3_GMAT (TRVI8)') 737 END IF 738C 739 CALL CCSDT_SRVIR3(WORK(KTRVI8),WORK(KEND4),ISYMD,D,ISINTBC) 740C 741C----------------------------------------------- 742C Integrals used in L3am. 743C----------------------------------------------- 744C 745 IOFF = ICKBD(ISCKB2,ISYMD) + NCKATR(ISCKB2)*(D - 1) + 1 746 IF (NCKATR(ISCKB2) .GT. 0) THEN 747 CALL GETWA2(LUDKBC3,FNDKBC3,WORK(KTRVI0),IOFF, 748 & NCKATR(ISCKB2)) 749 ENDIF 750C 751 CALL CCSDT_SRVIR3(WORK(KTRVI0),WORK(KEND4),ISYMD,D,ISINT2) 752C 753C------------------------------------------------------ 754C Read 2*C-E of integral used for L3 755C------------------------------------------------------ 756C 757 IOFF = ICKBD(ISCKB2,ISYMD) + NCKATR(ISCKB2)*(D - 1) + 1 758 IF (NCKATR(ISCKB2) .GT. 0) THEN 759 CALL GETWA2(LU3FOP2X,FN3FOP2X,WORK(KTRVI4),IOFF, 760 & NCKATR(ISCKB2)) 761 ENDIF 762C 763C----------------------------------------------------------- 764C Sort the integrals for s3am and for t3-bar 765C----------------------------------------------------------- 766C 767 CALL CCSDT_SRTVIR(WORK(KTRVI0),WORK(KTRVI2),WORK(KEND4), 768 * LWRK4,ISYMD,ISINT2) 769C 770 CALL CCSDT_SRTVIR(WORK(KTRVI4),WORK(KTRVI5),WORK(KEND4), 771 * LWRK4,ISYMD,ISINT2) 772C 773C 774 IF (IPRINT .GT. 55) THEN 775 XNORMVAL = DDOT(NCKATR(ISCKB2),WORK(KTRVI0),1, 776 * WORK(KTRVI0),1) 777 WRITE(LUPRI,*) 'Norm of TRVI0 ',XNORMVAL 778 ENDIF 779C 780 IF (IPRINT .GT. 55) THEN 781 XNORMVAL = DDOT(NCKATR(ISCKB2),WORK(KTRVI2),1, 782 * WORK(KTRVI2),1) 783 WRITE(LUPRI,*) 'Norm of TRVI2 ',XNORMVAL 784 ENDIF 785C 786 IF (IPRINT .GT. 55) THEN 787 XNORMVAL = DDOT(NCKATR(ISCKB2),WORK(KTRVI4),1, 788 * WORK(KTRVI4),1) 789 WRITE(LUPRI,*) 'Norm of TRVI4 ',XNORMVAL 790 ENDIF 791C 792 IF (IPRINT .GT. 55) THEN 793 XNORMVAL = DDOT(NCKATR(ISCKB2),WORK(KTRVI5),1, 794 * WORK(KTRVI5),1) 795 WRITE(LUPRI,*) 'Norm of TRVI5 ',XNORMVAL 796 ENDIF 797C 798C------------------------------------------------------ 799C Calculate virtual integrals used in q3am. 800C------------------------------------------------------ 801C 802 IOFF = ICKAD(ISYCKB,ISYMD) + NCKA(ISYCKB)*(D - 1) + 1 803 IF (NCKA(ISYCKB) .GT. 0) THEN 804 CALL GETWA2(LU3VI,FN3VI,WORK(KINTVI),IOFF, 805 * NCKA(ISYCKB)) 806 ENDIF 807C 808 CALL CCSDT_TRVIR(WORK(KINTVI),WORK(KTRVI3),WORK(KLAMDP), 809 * ISYMD,D,ISYMOP,WORK(KEND4),LWRK4) 810C 811 IF (LWRK3 .LT. NCKATR(ISYCKB)) THEN 812 CALL QUIT('Insufficient space for allocation in '// 813 * 'CC3_GMAT (TRVI3)') 814 END IF 815C 816C Can use kend3 since dont need the integrals anymore 817 CALL CCSDT_SRVIR3(WORK(KTRVI3),WORK(KEND3),ISYMD,D,ISINT2) 818C 819C--------------------------------------------------------------- 820C Read virtual integrals used in q3am/u3am for t3-bar. 821C--------------------------------------------------------------- 822C 823 IOFF = ICKBD(ISYCKB,ISYMD) + NCKATR(ISYCKB)*(D - 1) + 1 824 IF (NCKATR(ISYCKB) .GT. 0) THEN 825 CALL GETWA2(LU3FOPX,FN3FOPX,WORK(KTRVI6),IOFF, 826 * NCKATR(ISYCKB)) 827 ENDIF 828C 829 IF (LWRK3 .LT. NCKATR(ISYCKB)) THEN 830 CALL QUIT('Insufficient space for allocation in '// 831 & 'CC3_GMAT (TRVI6)') 832 END IF 833C 834C Can use kend3 since dont need the integrals anymore 835 CALL CCSDT_SRVIR3(WORK(KTRVI6),WORK(KEND4),ISYMD,D,ISINT2) 836C 837 IF (IPRINT .GT. 55) THEN 838 XNORMVAL = DDOT(NCKATR(ISCKB2),WORK(KTRVI6),1, 839 * WORK(KTRVI6),1) 840 WRITE(LUPRI,*) 'Norm of TRVI6 ',XNORMVAL 841 ENDIF 842C 843C--------------------- 844C Calculate. 845C--------------------- 846C 847 DO ISYMB = 1,NSYM 848C 849 ISYALJ = MULD2H(ISYMB,ISYML2) 850 ISAIJ2 = MULD2H(ISYMB,ISYRES) 851 ISYMBD = MULD2H(ISYMB,ISYMD) 852 ISCKIJ = MULD2H(ISYMBD,ISYMIM) 853C 854 IF (IPRINT .GT. 55) THEN 855C 856 WRITE(LUPRI,*) 'In CC3_LHTR_L3: ISYMD :',ISYMD 857 WRITE(LUPRI,*) 'In CC3_LHTR_L3: ISYMB :',ISYMB 858 WRITE(LUPRI,*) 'In CC3_LHTR_L3: ISYALJ:',ISYALJ 859 WRITE(LUPRI,*) 'In CC3_LHTR_L3: ISAIJ2:',ISAIJ2 860 WRITE(LUPRI,*) 'In CC3_LHTR_L3: ISYMBD:',ISYMBD 861 WRITE(LUPRI,*) 'In CC3_LHTR_L3: ISCKIJ:',ISCKIJ 862C 863 ENDIF 864C 865C Can use kend3 since we do not need the integrals anymore. 866 KSMAT = KEND3 867 KQMAT = KSMAT + NCKIJ(ISCKIJ) 868 KDIAG = KQMAT + NCKIJ(ISCKIJ) 869 KINDSQ = KDIAG + NCKIJ(ISCKIJ) 870 KINDEX = KINDSQ + (6*NCKIJ(ISCKIJ) - 1)/IRAT + 1 871 KTMAT = KINDEX + (NCKI(ISYALJ) - 1)/IRAT + 1 872 KRMAT2 = KTMAT + NCKIJ(ISCKIJ) 873 KEND4 = KRMAT2 + NCKI(ISAIJ2) 874 LWRK4 = LWORK - KEND4 875C 876 IF (LWRK4 .LT. 0) THEN 877 WRITE(LUPRI,*) 'Memory available : ',LWORK 878 WRITE(LUPRI,*) 'Memory needed : ',KEND4 879 CALL QUIT('Insufficient space in CC3_LHTR_L3 (inner)') 880 END IF 881C 882C--------------------------------------------- 883C Construct part of the diagonal. 884C--------------------------------------------- 885C 886 CALL CC3_DIAG(WORK(KDIAG),WORK(KFOCKD),ISCKIJ) 887C 888 IF (IPRINT .GT. 55) THEN 889 XNORMVAL = DDOT(NCKIJ(ISCKIJ),WORK(KDIAG),1, 890 * WORK(KDIAG),1) 891 WRITE(LUPRI,*) 'Norm of DIA ',XNORMVAL 892 ENDIF 893 894C 895C------------------------------------- 896C Construct index arrays. 897C------------------------------------- 898C 899 LENSQ = NCKIJ(ISCKIJ) 900 CALL CC3_INDSQ(WORK(KINDSQ),LENSQ,ISCKIJ) 901 CALL CC3_INDEX(WORK(KINDEX),ISYALJ) 902C 903 DO B = 1,NVIR(ISYMB) 904C 905C----------------------------------------- 906C Initialize the R2 matrix. 907C----------------------------------------- 908C 909 CALL DZERO(WORK(KRMAT2),NCKI(ISAIJ2)) 910C 911C--------------------------------------------------------------------------- 912C Calculate the S(ci,bk,dj) matrix for for B,D for L3 913C Scale the smat with a half since have a factor 914C of two in the routine. 915C--------------------------------------------------------------------------- 916C 917C 918 CALL DZERO(WORK(KSMAT),NCKIJ(ISCKIJ)) 919C 920 CALL CCFOP_SMAT(0.0D0,WORK(KL1AM),ISYML1,WORK(KL2TP), 921 * ISYML2,WORK(KTMAT), 922 * WORK(KFCKBA),WORK(KXIAJB),ISINT1, 923 * WORK(KTRVI0),WORK(KTRVI2), 924 * WORK(KTRVI4),WORK(KTRVI5), 925 * WORK(KTROC0),WORK(KTROC2), 926 * ISINT2,WORK(KFOCKD),WORK(KDIAG), 927 * WORK(KSMAT),WORK(KEND4),LWRK4, 928 * WORK(KINDEX),WORK(KINDSQ),LENSQ, 929 * ISYMB,B,ISYMD,D) 930C 931 CALL DSCAL(NCKIJ(ISCKIJ),HALF,WORK(KSMAT),1) 932C 933 CALL T3_FORBIDDEN(WORK(KSMAT),ISYMIM,ISYMB,B,ISYMD,D) 934COMMENT COMMENT 935COMMENT COMMENT 936C call sum_pt3(work(ksmat),isymb,b,isymd,d, 937C * isckij,work(kt3am),1) 938COMMENT COMMENT 939COMMENT COMMENT 940C 941C 942 IF (IPRINT .GT. 55) THEN 943 XNORMVAL = DDOT(NCKIJ(ISCKIJ),WORK(KSMAT),1, 944 * WORK(KSMAT),1) 945 WRITE(LUPRI,*) 'Norm of SMAT-L3 ',XNORMVAL 946 ENDIF 947C 948C------------------------------------------------------------------- 949C Calculate Q(ci,jk) for fixed b,d for t3-bar. 950C------------------------------------------------------------------- 951C 952C 953 CALL DZERO(WORK(KQMAT),NCKIJ(ISCKIJ)) 954C 955 ! g integrals 956c call dzero(WORK(KTRVI3),NCKATR(ISCKB2)) 957 ! L integrals 958c call dzero(WORK(KTRVI6),NCKATR(ISCKB2)) 959C 960 ! g integrals 961c call dzero(WORK(KTROC0),NTRAOC(ISINT2)) 962c ! L integrals 963c call dzero(WORK(KTROC2),NTRAOC(ISINT2)) 964c 965c call dzero(WORK(KXIAJB),NT2AM(ISYMOP)) 966c call dzero(WORK(KFCKBA),NT1AM(1)) 967 CALL CCFOP_QMAT(0.0D0,WORK(KL1AM),ISYML1, 968 * WORK(KL2TP),ISYML2, 969 * WORK(KTMAT),WORK(KFCKBA), 970 * WORK(KXIAJB),ISINT1,WORK(KTRVI3), 971 * WORK(KTRVI6),WORK(KTROC0), 972 * WORK(KTROC2),ISINT2,WORK(KFOCKD), 973 * WORK(KDIAG),WORK(KQMAT), 974 * WORK(KEND4),LWRK4,WORK(KINDSQ), 975 * LENSQ,ISYMB,B,ISYMD,D) 976C 977 CALL DSCAL(NCKIJ(ISCKIJ),HALF,WORK(KQMAT),1) 978C 979 CALL T3_FORBIDDEN(WORK(KQMAT),ISYMIM,ISYMB,B,ISYMD,D) 980COMMENT COMMENT COMMENT 981COMMENT COMMENT COMMENT 982c call sum_pt3(work(kqmat),isymb,b,isymd,d, 983c * isckij,work(kt3am),2) 984COMMENT COMMENT COMMENT 985COMMENT COMMENT COMMENT 986C 987C 988 IF (IPRINT .GT. 55) THEN 989 XNORMVAL = DDOT(NCKIJ(ISCKIJ),WORK(KQMAT),1, 990 * WORK(KQMAT),1) 991 WRITE(LUPRI,*) 'Norm of QMAT-L3 ',XNORMVAL 992 ENDIF 993C 994C----------------------------------------------------------------------- 995C Calculate the contributions to omega2 996C <L3|[H^BC,tau_{2}]|HF> 997C----------------------------------------------------------------------- 998c WRITE(lupri,*)'ome2 before CFOP_CONOCC' 999c call output(WORK(KOME2),1,nt1am(1),1,nt1am(1), 1000c * nt1am(1),nt1am(1),1,lupri) 1001C 1002 CALL CCFOP_CONVIR(WORK(KRMAT2),WORK(KSMAT), 1003 * WORK(KQMAT),WORK(KTMAT),ISYMIM, 1004 * WORK(KTRVI7),WORK(KTRVI8),ISINTBC, 1005 * WORK(KEND4),LWRK4,WORK(KINDSQ), 1006 * LENSQ,ISYMB,B,ISYMD,D) 1007c write(lupri,*)'ISYMB,B,ISYMD,D,ISINTBC',ISYMB,B,ISYMD,D,ISINTBC 1008c xnormval = ddot(NCKIJ(ISCKIJ),WORK(KSMAT),1,WORK(KSMAT),1) 1009c write(lupri,*)'KSMAT',xnormval 1010c xnormval = ddot(NCKIJ(ISCKIJ),WORK(KQMAT),1,WORK(KQMAT),1) 1011c write(lupri,*)'KQMAT',xnormval 1012c xnormval = ddot(NCKATR(ISCKBBC),WORK(KTRVI7),1,WORK(KTRVI7),1) 1013c write(lupri,*)'KTRVI7,ISINTBC',xnormval,ISINTBC 1014c xnormval = ddot(NCKATR(ISCKBBC),WORK(KTRVI8),1,WORK(KTRVI8),1) 1015c write(lupri,*)'KTRVI8',xnormval 1016c xnormval = ddot(NTRAOC(ISINTBC),WORK(KTROC3),1,WORK(KTROC3),1) 1017c write(lupri,*)'KTROC3',xnormval 1018c xnormval = ddot(NTRAOC(ISINTBC),WORK(KTROC4),1,WORK(KTROC4),1) 1019c write(lupri,*)'KTROC4',xnormval 1020C 1021 CALL CCFOP_CONOCC(WORK(KOME2),WORK(KRMAT1), 1022 * WORK(KRMAT2),WORK(KSMAT), 1023 * WORK(KTMAT),ISYMIM, 1024 * WORK(KTROC3),WORK(KTROC4),ISINTBC, 1025 * WORK(KEND4),LWRK4,WORK(KINDSQ), 1026 * LENSQ,ISYMB,B,ISYMD,D) 1027c XNORMVAL = DDOT(NT2AM(ISYRES),WORK(KOME2),1, 1028c * WORK(KOME2),1) 1029c WRITE(LUPRI,*) 'Norm Rho22-after CCFOP_CONOCC',XNORMVAL 1030c WRITE(lupri,*)'ome2 after CFOP_CONOCC' 1031c call output(WORK(KOME2),1,nt1am(1),1,nt1am(1), 1032c * nt1am(1),nt1am(1),1,lupri) 1033 1034C 1035C------------------------------------------------------------------------ 1036C Calculate the L3 contribution to omega1 1037C------------------------------------------------------------------------ 1038C 1039 CALL DSCAL(NCKIJ(ISCKIJ),-ONE,WORK(KSMAT),1) 1040 CALL DSCAL(NCKIJ(ISCKIJ),-ONE,WORK(KQMAT),1) 1041C 1042C <L3|[H^BC,T^{0}_{2}],tau_{1}]|HF> 1043c write(lupri,*)'ISYMB,B,ISYMD,D,ISINTBC' 1044c write(lupri,*)ISYMB,B,ISYMD,D,ISINTBC 1045C 1046c XNORMVAL = DDOT(N3ORHF(ISINTBC),WORK(KBCIOOOO),1, 1047c * WORK(KBCIOOOO),1) 1048c write(lupri,*)'KBCIOOOO',xnormval 1049c XNORMVAL = DDOT(NT2SQ(ISINTBC),WORK(KBCIOVVO),1, 1050c * WORK(KBCIOVVO),1) 1051c write(lupri,*)'KBCIOVVO',xnormval 1052c XNORMVAL = DDOT(NT2SQ(ISINTBC),WORK(KBCIOOVV),1, 1053c * WORK(KBCIOOVV),1) 1054c write(lupri,*)'KBCIOOVV',xnormval 1055c XNORMVAL = DDOT(NMAABC(ISCKBBC),WORK(KVVVVBC),1, 1056c * WORK(KVVVVBC),1) 1057c write(lupri,*)'KVVVVBC',xnormval 1058 !<L3|[H^BC,T^{0}_{2}],tau_{1}]|HF> 1059 CALL CC3_L3_OMEGA1(WORK(KOME1),ISYRES,WORK(KSMAT), 1060 * WORK(KQMAT),WORK(KTMAT),ISYMIM, 1061 * WORK(KBCIOOOO),WORK(KBCIOVVO), 1062 * WORK(KBCIOOVV),WORK(KVVVVBC), 1063 * ISINTBC, 1064 * WORK(KT2TP),ISYMT2,WORK(KEND4), 1065 * LWRK4,LENSQ,WORK(KINDSQ), 1066 * ISYMB,B,ISYMD,D) 1067C 1068c WRITE(lupri,*)'ome1 after CC3_L3_OMEGA2_NODDY t20 cont' 1069c call output(WORK(KOME1),1,nvir(1),1,nrhf(1), 1070c * nvir(1),nrhf(1),1,lupri) 1071 1072 IF (IPRINT .GT. 55) THEN 1073 XNORMVAL = DDOT(NT1AM(ISYRES),OMEGA1,1,OMEGA1,1) 1074 WRITE(LUPRI,*) 'Norm CC3_GMAT after CC3_L3_OMEGA1', 1075 * XNORMVAL 1076 ENDIF 1077C 1078C <L3|[H^B,T^{C}_{2}],tau_{1}]|HF> 1079C 1080C 1081 CALL CC3_L3_OMEGA1(WORK(KOME1),ISYRES,WORK(KSMAT), 1082 * WORK(KQMAT),WORK(KTMAT),ISYMIM, 1083 * WORK(KBIOOOO),WORK(KBIOVVO), 1084 * WORK(KBIOOVV),WORK(KVVVVB),ISINTB, 1085 * WORK(KC2TP),ISYMC2,WORK(KEND4), 1086 * LWRK4,LENSQ,WORK(KINDSQ), 1087 * ISYMB,B,ISYMD,D) 1088c WRITE(lupri,*)'ome1 after CC3_L3_OMEGA2_NODDY t2c cont' 1089c call output(WORK(KOME1),1,nvir(1),1,nrhf(1), 1090c * nvir(1),nrhf(1),1,lupri) 1091 1092C 1093 IF (IPRINT .GT. 55) THEN 1094 XNORMVAL = DDOT(NT1AM(ISYRES),OMEGA1,1,OMEGA1,1) 1095 WRITE(LUPRI,*) 'Norm CC3_GMAT after CC3_L3_OMEGA1', 1096 * XNORMVAL 1097 END IF 1098 1099C 1100 1101C <L3|[H^C,T^{B}_{2}],tau_{1}]|HF> 1102C 1103 CALL CC3_L3_OMEGA1(WORK(KOME1),ISYRES,WORK(KSMAT), 1104 * WORK(KQMAT),WORK(KTMAT),ISYMIM, 1105 * WORK(KCIOOOO),WORK(KCIOVVO), 1106 * WORK(KCIOOVV),WORK(KVVVVC),ISINTC, 1107 * WORK(KB2TP),ISYMB2,WORK(KEND4), 1108 * LWRK4,LENSQ,WORK(KINDSQ), 1109 * ISYMB,B,ISYMD,D) 1110c WRITE(lupri,*)'ome1 after CC3_L3_OMEGA2_NODDY t2b cont' 1111c call output(WORK(KOME1),1,nvir(1),1,nrhf(1), 1112c * nvir(1),nrhf(1),1,lupri) 1113 1114C 1115 IF (IPRINT .GT. 55) THEN 1116 XNORMVAL = DDOT(NT1AM(ISYRES),OMEGA1,1,OMEGA1,1) 1117 WRITE(LUPRI,*) 'Norm CC3_GMAT after CC3_L3_OMEGA1', 1118 * XNORMVAL 1119 ENDIF 1120C 1121C------------------------------------------- 1122C Accumulate R2 in Omega2 1123C------------------------------------------- 1124C 1125 CALL CC3_RACC(WORK(KOME2),WORK(KRMAT2),ISYMB,B, 1126 * ISYRES) 1127c WRITE(lupri,*)'ome2 after cc3_racc rmat2' 1128c call output(WORK(KOME2),1,nt1am(1),1,nt1am(1), 1129c * nt1am(1),nt1am(1),1,lupri) 1130c WRITE(lupri,*)'ome2 triangular form ' 1131c call outpak(WORK(KOME2),nt1amx,1,lupri) 1132C 1133 IF (IPRINT .GT. 55) THEN 1134 XNORMVAL = DDOT(NT2AM(ISYRES),WORK(KOME2),1, 1135 * WORK(KOME2),1) 1136 WRITE(LUPRI,*) 'Norm of Rho22-after CC3_RACC', 1137 * XNORMVAL 1138 ENDIF 1139C 1140 IF (IPRINT .GT. 220) THEN 1141 CALL AROUND('After CC3_RACC: ') 1142 CALL CC_PRP(DUMMY,OMEGA2,ISYRES,0,1) 1143 ENDIF 1144C 1145 ENDDO ! B 1146 ENDDO ! ISYMB 1147C 1148C--------------------------------------- 1149C Accumulate R1 in Omega2. 1150C--------------------------------------- 1151C 1152 CALL CC3_RACC(WORK(KOME2),WORK(KRMAT1),ISYMD,D,ISYRES) 1153c WRITE(lupri,*)'ome2 after cc3_racc rmat1' 1154c call output(WORK(KOME2),1,nt1am(1),1,nt1am(1), 1155c * nt1am(1),nt1am(1),1,lupri) 1156C 1157 IF (IPRINT .GT. 55) THEN 1158 XNORMVAL = DDOT(NT2AM(ISYRES),WORK(KOME2),1, 1159 * WORK(KOME2),1) 1160 WRITE(LUPRI,*) 'Norm of Rho22-after CC3_RACC-2',XNORMVAL 1161 ENDIF 1162C 1163 IF (IPRINT .GT. 220) THEN 1164 CALL AROUND('After CC3_RACC-2: ') 1165 CALL CC_PRP(DUMMY,WORK(KOME2),ISYRES,0,1) 1166 ENDIF 1167C 1168 ENDDO ! D 1169 ENDDO ! ISYMD 1170C 1171COMMENT COMMENT COMMENT COMMENT COMMENT 1172COMMENT COMMENT COMMENT COMMENT COMMENT 1173COMMENT COMMENT COMMENT COMMENT COMMENT 1174c write(lupri,*) 'L3 Amplitudes in cc3_Gmat : ' 1175C call dscal(nt1amx*nt1amx*nt1amx,-1.0D0,work(kt3am),1) 1176c call print_pt3(work(kt3am),isckij,3) 1177COMMENT COMMENT COMMENT COMMENT COMMENT 1178COMMENT COMMENT COMMENT COMMENT COMMENT 1179COMMENT COMMENT COMMENT COMMENT COMMENT 1180C------------------------------- 1181C Close and delete files 1182C------------------------------- 1183C 1184 CALL WCLOSE2(LU4VBC,FN4VBC,'DELETE') 1185 CALL WCLOSE2(LU4VC,FN4VC,'DELETE') 1186 CALL WCLOSE2(LU4VB,FN4VB,'DELETE') 1187 CALL WCLOSE2(LU3SRTR,FN3SRTR,'DELETE') 1188 CALL WCLOSE2(LUCKJDR,FNCKJDR,'DELETE') 1189 CALL WCLOSE2(LUDELDR,FNDELDR,'DELETE') 1190 CALL WCLOSE2(LUDKBCR,FNDKBCR,'DELETE') 1191 CALL WCLOSE2(LUDKBCR4,FNDKBCR4,'DELETE') 1192C 1193COMMENT COMMENT COMMENT COMMENT 1194C write(lupri,*) 'OMEGA1EFF in cc3_Gmat before added' 1195C do i = 1,nt1am(isyres) 1196C write(lupri,*) i,OMEGA1EFF(i) 1197C end do 1198COMMENT COMMENT COMMENT COMMENT 1199C write(lupri,*) 'OMEGA2EFF in cc3_Gmat before added' 1200C call outpak(OMEGA2EFF,NT2AM(isyres),1,lupri) 1201COMMENT COMMENT COMMENT COMMENT 1202 1203C 1204C WRITE(LUPRI,*)' OME1 FINAL' 1205C CALL PRINT_MATAI(WORK(KOME1),ISYRES) 1206C XNORMVAL = DDOT(NT1AM(ISYRES),WORK(KOME1),1, 1207C * WORK(KOME1),1) 1208C WRITE(LUPRI,*) 'Norm of OME1',XNORMVAL 1209C 1210 DO I = 1,NT1AM(ISYRES) 1211 OMEGA1(I) = OMEGA1(I) + WORK(KOME1-1+I) 1212 END DO 1213C 1214c WRITE(LUPRI,*)' OMEGA1 FINAL' 1215c CALL PRINT_MATAI(OMEGA1,ISYRES) 1216c XNORMVAL = DDOT(NT1AM(ISYRES),OMEGA1,1, 1217c * OMEGA1,1) 1218c WRITE(LUPRI,*) 'Norm of OMEGA1',XNORMVAL 1219 1220c WRITE(lupri,*)'ome2 triangular form final ' 1221c call outpak(WORK(KOME2),nt1amx,1,lupri) 1222C XNORMVAL = DDOT(NT2AM(ISYRES),WORK(KOME2),1, 1223C * WORK(KOME2),1) 1224C WRITE(LUPRI,*) 'Norm of ome2 ',XNORMVAL 1225C 1226 DO I = 1,NT2AM(ISYRES) 1227 OMEGA2(I) = OMEGA2(I) + WORK(KOME2-1+I) 1228 END DO 1229C 1230c WRITE(lupri,*)'omega2 triangular form final ' 1231c call outpak(OMEGA2,nt1amx,1,lupri) 1232C XNORMVAL = DDOT(NT2AM(ISYRES),OMEGA2,1, 1233C * OMEGA2,1) 1234C WRITE(LUPRI,*) 'Norm of OMEGA2',XNORMVAL 1235C 1236COMMENT COMMENT COMMENT COMMENT 1237C write(lupri,*) 'OMEGA1EFF in cc3_Gmat after added' 1238C do i = 1,nt1am(isyres) 1239C write(lupri,*) i,OMEGA1EFF(i) 1240C end do 1241COMMENT COMMENT COMMENT COMMENT 1242C write(lupri,*) 'OMEGA2EFF in cc3_Gmat after added' 1243C call outpak(OMEGA2EFF,NT2AM(isyres),1,lupri) 1244COMMENT COMMENT COMMENT COMMENT 1245 1246C 1247C------------- 1248C End 1249C------------- 1250C 1251 CALL QEXIT('CC3_GMAT') 1252C 1253 RETURN 1254C 1255 1 FORMAT(7X,'Time used in',2X,A12,F12.2,' seconds') 1256C 1257 END 1258C /* Deck cc3_3barint */ 1259 SUBROUTINE CC3_3BARINT(ISYMX,LISTX,IDLSTX, 1260 * ISYMY,LISTY,IDLSTY, 1261 * ISYMZ1,LISTZ,IDLSTZ,Z1EXIST, 1262 * XLAMDP,XLAMDH,WORK,LWORK, 1263 * LU3SRTR,FN3SRTR,LUCKJDR,FNCKJDR) 1264C 1265C 1266C Interface to calculate 1267C the (ck|de)-X,Y,Z bar and (ck|lm)-X,Y,Z bar integrals 1268C 1269C e and l index is normal T1-transformed (and can be transformed outside) 1270C 1271C (alfa beta /gamma del) is transformed to (ck/d del) and (ck/del m) 1272C 1273C if Z1EXIST eq true 1274C one indes transform with z1am 1275C else 1276C use XLAMDP,XLAMDH as transformation matrix 1277C endif 1278C 1279C The atomic integrals (alfa beta /gamma del) have symmetry isymop 1280C 1281 IMPLICIT NONE 1282C 1283#include "ccsdsym.h" 1284#include "priunit.h" 1285#include "dummy.h" 1286#include "maxorb.h" 1287#include "maxash.h" 1288#include "mxcent.h" 1289#include "aovec.h" 1290#include "iratdef.h" 1291#include "ccorb.h" 1292#include "ccisao.h" 1293#include "r12int.h" 1294#include "blocks.h" 1295#include "ccsdinp.h" 1296#include "distcl.h" 1297#include "cbieri.h" 1298#include "eritap.h" 1299C#include "cclr.h" 1300C 1301 INTEGER ISYMX,ISYMY,ISYMZ,ISYMZ1,LWORK,LU3SRTR,LUCKJDR 1302 INTEGER KLAMPX,KLAMHX,KLAMPY,KLAMHY,KLAMPZ,KEND1,LWRK1 1303 INTEGER KCCFB1, KINDXB, KODCL1, KODCL2, KODBC1, KODBC2 1304 INTEGER KRDBC1, KRDBC2, KODPP1, KODPP2, KRDPP1, KRDPP2 1305 INTEGER KFREE,LFREE 1306 INTEGER NTOSYM,NTOT,ILLL,NUMDIS,KRECNR,IDEL2,IDEL,ISYMD 1307 INTEGER ISYDIS,ISYAIK,ISYCKJ,KXINT,KCKJ,KEND2,LWRK2 1308 INTEGER KLAMHZ,ISYMD1 1309 INTEGER IDLSTX,IDLSTY,IDLSTZ 1310 INTEGER IDUM 1311C 1312 INTEGER INDEXA(MXCORB_CC) 1313#if defined (SYS_CRAY) 1314 REAL XLAMDP(*),XLAMDH(*) 1315 REAL WORK(LWORK) 1316 real ddot 1317#else 1318 DOUBLE PRECISION XLAMDP(*),XLAMDH(*) 1319 DOUBLE PRECISION WORK(LWORK) 1320 DOUBLE PRECISION ddot 1321#endif 1322C 1323 CHARACTER*3 LISTX, LISTY, LISTZ 1324 CHARACTER*(*) FN3SRTR, FNCKJDR 1325 CHARACTER*1 CDUMMY 1326C 1327 LOGICAL Z1EXIST 1328C 1329 CALL QENTER('CC3_3BARINT') 1330C 1331 CDUMMY = ' ' 1332C 1333 KLAMPX = 1 1334 KLAMHX = KLAMPX + NLAMDT 1335 KLAMPY = KLAMHX + NLAMDT 1336 KLAMHY = KLAMPY + NLAMDT 1337 KLAMPZ = KLAMHY + NLAMDT 1338 KLAMHZ = KLAMPZ + NLAMDT 1339 KEND1 = KLAMHZ + NLAMDT 1340 LWRK1 = LWORK - KEND1 1341C 1342C 1343C-------------------------------------- 1344C Calculate lambda-X, lambda-Y, lambda-Z 1345C-------------------------------------- 1346C 1347c CALL CCLR_LAMTRA(XLAMDP,WORK(KLAMPX),XLAMDH, 1348c * WORK(KLAMHX),X1AM,ISYMX) 1349c CALL CCLR_LAMTRA(XLAMDP,WORK(KLAMPY),XLAMDH, 1350c * WORK(KLAMHY),Y1AM,ISYMY) 1351 CALL GET_LAMBDAX(WORK(KLAMPX),WORK(KLAMHX),LISTX,IDLSTX, 1352 * ISYMX,XLAMDP,XLAMDH,WORK(KEND1),LWRK1) 1353 CALL GET_LAMBDAX(WORK(KLAMPY),WORK(KLAMHY),LISTY,IDLSTY, 1354 * ISYMY,XLAMDP,XLAMDH,WORK(KEND1),LWRK1) 1355 1356 IF (Z1EXIST) THEN 1357 ISYMZ = ISYMZ1 1358 CALL GET_LAMBDAX(WORK(KLAMPZ),WORK(KLAMHZ),LISTZ,IDLSTZ, 1359 * ISYMZ,XLAMDP,XLAMDH,WORK(KEND1),LWRK1) 1360C CALL CCLR_LAMTRA(XLAMDP,WORK(KLAMPZ),XLAMDH, 1361C * WORK(KLAMHZ),Z1AM,ISYMZ) 1362 ELSE 1363 CALL DZERO(WORK(KLAMPZ),NLAMDT) 1364 CALL DZERO(WORK(KLAMHZ),NLAMDT) 1365 ISYMZ = 1 1366 CALL DCOPY(NLAMDT,XLAMDP,1,WORK(KLAMPZ),1) 1367 CALL DCOPY(NLAMDT,XLAMDH,1,WORK(KLAMHZ),1) 1368 END IF 1369C 1370C-------------------------------- 1371C Calculate integrals 1372C-------------------------------- 1373C 1374 IF (DIRECT) THEN 1375 IF (HERDIR) THEN 1376 CALL HERDI1(WORK(KEND1),LWRK1,IPRERI) 1377 ELSE 1378 KCCFB1 = KEND1 1379 KINDXB = KCCFB1 + MXPRIM*MXCONT 1380 KEND1 = KINDXB + (8*MXSHEL*MXCONT + 1)/IRAT 1381 LWRK1 = LWORK - KEND1 1382 CALL ERIDI1(KODCL1,KODCL2,KODBC1,KODBC2,KRDBC1,KRDBC2, 1383 & KODPP1,KODPP2,KRDPP1,KRDPP2, 1384 & KFREE,LFREE,KEND1,WORK(KCCFB1),WORK(KINDXB), 1385 & WORK(KEND1),LWRK1,IPRERI) 1386 KEND1 = KFREE 1387 LWRK1 = LFREE 1388 ENDIF 1389 NTOSYM = 1 1390 ELSE 1391 NTOSYM = NSYM 1392 ENDIF 1393C 1394 DO ISYMD1 = 1,NTOSYM 1395C 1396 IF (DIRECT) THEN 1397 IF (HERDIR) THEN 1398 NTOT = MAXSHL 1399 ELSE 1400 NTOT = MXCALL 1401 ENDIF 1402 ELSE 1403 NTOT = NBAS(ISYMD1) 1404 ENDIF 1405C 1406 DO ILLL = 1,NTOT 1407C 1408C------------------------------------------ 1409C If direct calculate the integrals. 1410C------------------------------------------ 1411C 1412 IF (DIRECT) THEN 1413C 1414 IF (HERDIR) THEN 1415 CALL HERDI2(WORK(KEND1),LWRK1,INDEXA,ILLL,NUMDIS, 1416 & IPRERI) 1417 ELSE 1418 CALL ERIDI2(ILLL,INDEXA,NUMDIS,0,0, 1419 & WORK(KODCL1),WORK(KODCL2), 1420 & WORK(KODBC1),WORK(KODBC2), 1421 & WORK(KRDBC1),WORK(KRDBC2), 1422 & WORK(KODPP1),WORK(KODPP2), 1423 & WORK(KRDPP1),WORK(KRDPP2), 1424 & WORK(KCCFB1),WORK(KINDXB), 1425 & WORK(KEND1), LWRK1,IPRERI) 1426 ENDIF 1427C 1428 KRECNR = KEND1 1429 KEND1 = KRECNR + (NBUFX(0) - 1)/IRAT + 1 1430 LWRK1 = LWORK - KEND1 1431 IF (LWRK1 .LT. 0) THEN 1432 CALL QUIT('Insufficient core in CC3_3BARINT') 1433 END IF 1434C 1435 ELSE 1436 NUMDIS = 1 1437 ENDIF 1438C 1439C-------------------------------------------------- 1440C Loop over number of distributions in disk. 1441C-------------------------------------------------- 1442C 1443 DO IDEL2 = 1,NUMDIS 1444C 1445 IF (DIRECT) THEN 1446 IDEL = INDEXA(IDEL2) 1447 IF (NOAUXB) THEN 1448 IDUM = 1 1449 CALL IJKAUX(IDEL,IDUM,IDUM,IDUM) 1450 END IF 1451 ISYMD = ISAO(IDEL) 1452 ELSE 1453 IDEL = IBAS(ISYMD1) + ILLL 1454 ISYMD = ISYMD1 1455 ENDIF 1456C 1457 ISYDIS = MULD2H(ISYMD,ISYMOP) 1458 ISYAIK = MULD2H(ISYDIS,ISYMOP) 1459 ISYCKJ = ISYDIS 1460C 1461C------------------------------------------ 1462C Work space allocation no. 2. 1463C------------------------------------------ 1464C 1465 KXINT = KEND1 1466 KCKJ = KXINT + NDISAO(ISYDIS) 1467 KEND2 = KCKJ + NCKI(ISYCKJ) 1468 LWRK2 = LWORK - KEND2 1469C 1470 IF (LWRK2 .LT. 0) THEN 1471 WRITE(LUPRI,*) 'Need : ',KEND2,'Available : ',LWORK 1472 CALL QUIT('Insufficient space in CC3_GMAT ') 1473 ENDIF 1474C 1475 CALL DZERO(WORK(KCKJ),NCKI(ISYCKJ)) 1476C 1477C----------------------------------------- 1478C Read in batch of integrals. 1479C----------------------------------------- 1480C 1481 CALL CCRDAO(WORK(KXINT),IDEL,IDEL2,WORK(KEND2),LWRK2, 1482 * WORK(KRECNR),DIRECT) 1483C 1484C----------------------------------------------------------- 1485C Calculate transformed X,Y,Z bar integrals 1486C----------------------------------------------------------- 1487C 1488 CALL CC3_BARXYZINT(WORK(KXINT),XLAMDP,XLAMDH 1489 * ,WORK(KLAMPX),WORK(KLAMHX),ISYMX 1490 * ,WORK(KLAMPY),WORK(KLAMHY),ISYMY 1491 * ,WORK(KLAMPZ),WORK(KLAMHZ),ISYMZ, 1492 * WORK(KEND2),LWRK2, 1493 * IDEL,ISYMD,LU3SRTR,FN3SRTR, 1494 * LUCKJDR,FNCKJDR) 1495c write(lupri,*)'IDEL,ISYMD',IDEL,ISYMD 1496c write(lupri,*)'KXINT',ddot(NDISAO(ISYDIS),WORK(KXINT), 1497c * 1,WORK(KXINT),1) 1498C 1499 ENDDO ! IDEL2 1500 ENDDO ! ILLL 1501 ENDDO ! ISYMD1 1502C 1503 CALL QEXIT('CC3_3BARINT') 1504 RETURN 1505 END 1506C /* Deck cc3_barxyzint */ 1507 SUBROUTINE CC3_BARXYZINT(XINT,XLAMDP,XLAMDH,XLAMPX,XLAMHX,ISYMX, 1508 * XLAMPY,XLAMHY,ISYMY,XLAMPZ,XLAMHZ,ISYMZ, 1509 * WORK,LWORK,IDEL,ISYDEL, 1510 * LUOUT1,FNOUT1,LUOUT2,FNOUT2) 1511C 1512C Purpose: Calculate X,Y,Z BAR integrals used in CC3. 1513C 1514 IMPLICIT NONE 1515C 1516#include "priunit.h" 1517#include "ccinftap.h" 1518#include "ccorb.h" 1519#include "ccsdsym.h" 1520C 1521 INTEGER ISYMX,ISYMY,ISYMZ,LWORK,IDEL,ISYDEL,LUOUT1,LUOUT2 1522 INTEGER KXCKD,KXCKJ,ID,LENGTH,IOFF 1523 INTEGER KEND1,LWRK1 1524 INTEGER ISYMXY,ISYMXYZ,ISYCKD,ISYCKJ 1525#if defined (SYS_CRAY) 1526 REAL XINT(*), XLAMDP(*),XLAMDH(*),XLAMPX(*),XLAMHX(*) 1527 REAL XLAMPY(*),XLAMHY(*),XLAMPZ(*),XLAMHZ(*) 1528 REAL WORK(LWORK) 1529 real xnormval,ddot 1530#else 1531 DOUBLE PRECISION XINT(*), XLAMDP(*),XLAMDH(*),XLAMPX(*),XLAMHX(*) 1532 DOUBLE PRECISION XLAMPY(*),XLAMHY(*),XLAMPZ(*),XLAMHZ(*) 1533 DOUBLE PRECISION WORK(LWORK) 1534 DOUBLE PRECISION xnormval,ddot 1535#endif 1536 1537C 1538 CHARACTER*(*) FNOUT1, FNOUT2 1539C 1540 CALL QENTER('CC3_BARXYZINT') 1541C 1542c xnormval = ddot(nlamdt,XLAMDP,1,XLAMDP,1) 1543c write(lupri,*)'XLAMDP',xnormval 1544c xnormval = ddot(nlamdt,XLAMDH,1,XLAMDH,1) 1545c write(lupri,*)'XLAMDH',xnormval 1546c xnormval = ddot(nglmdt(isymx),XLAMPX,1,XLAMPX,1) 1547c write(lupri,*)'XLAMPX',xnormval,isymx 1548c xnormval = ddot(nglmdt(isymx),XLAMHX,1,XLAMHX,1) 1549c write(lupri,*)'XLAMHX',xnormval,isymx 1550c xnormval = ddot(nglmdt(isymY),XLAMPY,1,XLAMPY,1) 1551c write(lupri,*)'XLAMPY',xnormval,isymY 1552c xnormval = ddot(nglmdt(isymY),XLAMHY,1,XLAMHY,1) 1553c write(lupri,*)'XLAMHY',xnormval,isymY 1554c xnormval = ddot(nglmdt(isymZ),XLAMPZ,1,XLAMPZ,1) 1555c write(lupri,*)'XLAMPZ',xnormval,isymz 1556c xnormval = ddot(nglmdt(isymz),XLAMHZ,1,XLAMHZ,1) 1557c write(lupri,*)'XLAMHZ',xnormval,isymz 1558C 1559 ISYMXY = MULD2H(ISYMX,ISYMY) 1560 ISYMXYZ = MULD2H(ISYMXY,ISYMZ) 1561 ISYCKD = MULD2H(ISYMXYZ,MULD2H(ISYDEL,ISYMOP)) 1562 ISYCKJ = ISYCKD 1563C 1564C--------------------------------- 1565C Allocation of work space. 1566C--------------------------------- 1567C 1568 KXCKD = 1 1569 KXCKJ = KXCKD + NCKATR(ISYCKD) 1570 KEND1 = KXCKJ + NCKI(ISYCKJ) 1571 LWRK1 = LWORK - KEND1 1572C 1573 IF (LWRK1 .LT. 0) THEN 1574 CALL QUIT('Insufficient core in CC3_BARXYZINT') 1575 ENDIF 1576C 1577C---------------------------------------- 1578C Calculate transformed integrals. 1579C---------------------------------------- 1580C 1581 CALL DZERO(WORK(KXCKD),NCKATR(ISYCKD)) 1582 CALL DZERO(WORK(KXCKJ),NCKI(ISYCKJ)) 1583C 1584 CALL CC3_CKD1(XINT,WORK(KXCKD),WORK(KXCKJ), 1585 * XLAMPX(NT1AO(ISYMX)+1),XLAMHY,XLAMPZ(NT1AO(ISYMZ)+1), 1586 * XLAMHZ,ISYMX,ISYMY,ISYMZ,ISYMZ, 1587 * WORK(KEND1),LWRK1,IDEL,ISYDEL) 1588C 1589 CALL CC3_CKD1(XINT,WORK(KXCKD),WORK(KXCKJ), 1590 * XLAMPX(NT1AO(ISYMX)+1),XLAMHZ,XLAMPY(NT1AO(ISYMY)+1), 1591 * XLAMHY,ISYMX,ISYMZ,ISYMY,ISYMY, 1592 * WORK(KEND1),LWRK1,IDEL,ISYDEL) 1593C 1594 CALL CC3_CKD1(XINT,WORK(KXCKD),WORK(KXCKJ), 1595 * XLAMPY(NT1AO(ISYMY)+1),XLAMHX,XLAMPZ(NT1AO(ISYMZ)+1), 1596 * XLAMHZ,ISYMY,ISYMX,ISYMZ,ISYMZ, 1597 * WORK(KEND1),LWRK1,IDEL,ISYDEL) 1598C 1599 CALL CC3_CKD1(XINT,WORK(KXCKD),WORK(KXCKJ), 1600 * XLAMPZ(NT1AO(ISYMZ)+1),XLAMHX,XLAMPY(NT1AO(ISYMY)+1), 1601 * XLAMHY,ISYMZ,ISYMX,ISYMY,ISYMY, 1602 * WORK(KEND1),LWRK1,IDEL,ISYDEL) 1603C 1604 CALL CC3_CKD1(XINT,WORK(KXCKD),WORK(KXCKJ), 1605 * XLAMPY(NT1AO(ISYMY)+1),XLAMHZ,XLAMPX(NT1AO(ISYMX)+1), 1606 * XLAMHX,ISYMY,ISYMZ,ISYMX,ISYMX, 1607 * WORK(KEND1),LWRK1,IDEL,ISYDEL) 1608C 1609 CALL CC3_CKD1(XINT,WORK(KXCKD),WORK(KXCKJ), 1610 * XLAMPZ(NT1AO(ISYMZ)+1),XLAMHY,XLAMPX(NT1AO(ISYMX)+1), 1611 * XLAMHX,ISYMZ,ISYMY,ISYMX,ISYMX, 1612 * WORK(KEND1),LWRK1,IDEL,ISYDEL) 1613C 1614C-------------------------------- 1615C Write to disk (ck|d alpha). 1616C-------------------------------- 1617C 1618 ID = IDEL - IBAS(ISYDEL) 1619C 1620 LENGTH = NCKATR(ISYCKD) 1621C 1622 IOFF = ICKDAO(ISYCKD,ISYDEL) + NCKATR(ISYCKD)*(ID - 1) + 1 1623C 1624 IF (LENGTH .GT. 0) THEN 1625 CALL PUTWA2(LUOUT1,FNOUT1,WORK(KXCKD),IOFF,LENGTH) 1626 ENDIF 1627c write(lupri,*)'id,isydel',id,isydel 1628C xnormval = ddot(NCKATR(ISYCKD),WORK(KXCKD),1,WORK(KXCKD),1) 1629C write(lupri,*)'KXCKD',xnormval 1630C 1631 LENGTH = NCKI(ISYCKJ) 1632C 1633 IOFF = ICKID(ISYCKJ,ISYDEL) + NCKI(ISYCKJ)*(ID - 1) + 1 1634C 1635 IF (LENGTH .GT. 0) THEN 1636 CALL PUTWA2(LUOUT2,FNOUT2,WORK(KXCKJ),IOFF,LENGTH) 1637 ENDIF 1638C xnormval = ddot(NCKI(ISYCKJ),WORK(KXCKJ),1,WORK(KXCKJ),1) 1639C write(lupri,*)'KXCKJ',xnormval 1640C 1641 CALL QEXIT('CC3_BARXYZINT') 1642 RETURN 1643 END 1644