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 19C /* Deck CC3_XISD */ 20 SUBROUTINE CC3_XISD(XI1,XI2,XI1EFF,XI2EFF,ISYRES, 21 * FOCKY,ISYFKY,FREQY,ICAU,LCAUCHY,LABELB, 22 * XLAMDP,XLAMDH,FOCK0, 23 * LUCKJD,FNCKJD,LUDKBC,FNDKBC,LUDELD,FNDELD, 24 * LUTOC,FNTOC,LU3VI,FN3VI,LU3VI2,FN3VI2, 25 * WORK,LWORK) 26C 27C Written by F. Pawlowski and P Jorgensen , Spring 2002. 28C 29C Partioning the triples part of the right hand side amplitude 30C equation \xi^Y into the singles and doubles space. 31C 32C Initially construct 33C 34C t3(ai,Bj,Ck) = S^BC(aikj) + U^BC(aikj) + S^CB(aijk) + U^CB(aijk) 35C 36C W^BC(aikj) = sum_d X(ad)*t3(di,Bj,Ck) - sum_l X(li)*t3(al,Bj,Ck) 37C - sum_ld X(ld)*t2(ai,Cl)*t2(Bj,dk) 38C - sum_ld X(ld)*t2(ai,Bl)*t2(Ck,dj) 39C 40C Kasper Hald, Summer 2002 41C Generalized to calculate the entire T3^Y and 42C calculate the contributions for F. 43C add ( -<mu3 | [[H,T1^Y],T2] | HF> -<mu3| [H,T2^Y] |HF> ) 44C 45C 46C Filip Pawlowski, Spring 2004, Aarhus, generalized to treat Cauchy vectors. 47C 48 IMPLICIT NONE 49C 50 LOGICAL LCAUCHY 51C 52 INTEGER ISYRES, ISINT1, ISINT2, ISYMT3, ISYMW 53 INTEGER ISYFKY, ISYMT2, IOPT, KFOCKD, KOMG1 54 INTEGER KOMG22, KFCKBA, KT2AM, KEND0, LWRK0, LWORK 55 INTEGER ISYMC, ISYMK, ISYMD, ISYMB 56 INTEGER KOFF1, KOFF2 57 INTEGER KTROC0, KTROC02, KEND1, LWRK1 58 INTEGER KINTOC, LWRK2, KEND2 59 INTEGER IOFF, ISAIJ1, ISAIJ2, KRMAT1, KRMAT2, KFCKY 60 INTEGER ISCKB1, ISCKB2 61 INTEGER KTRVI, KTRVI3, KTRVI1, KTRVI2, KTRVI0, KEND3, LWRK3 62 INTEGER KINTVI, KEND4, LWRK4 63 INTEGER ISYALJ, ISYALJ2, ISYMBD, ISCKIJ, ISCKD2 64 INTEGER KSMAT, KSMAT3, KUMAT, KUMAT3, KDIAG, KINDSQ 65 INTEGER KINDEX, KINDEX2, KTMAT, KTRVI8, KTRVI9 66 INTEGER KTRVI10, KEND5, LWRK5, LENSQ, LUFCK 67 INTEGER LUCKJD, LUDKBC, LUDELD, LUTOC, LU3VI, LU3VI2 68 INTEGER KWMAT, ISWMAT, LENSQW, KINDSQW 69 INTEGER KXIAJB, KDIAGW, KTROC, KTROC1 70 INTEGER LENGTH, ISYOPE, IOPTTCME, KRBJIA 71 INTEGER KT1B, KT2B 72 INTEGER LU3SRTR, LUCKJDR, LUDELDR, LUDKBCR, LUDKBCR4 73 INTEGER KTROC0Y, KTRVI0Y, KTRVI8Y 74C 75 INTEGER ICAU 76C 77 integer kx3am 78C 79 INTEGER IDLSTB,ISYMBB,NCAU,KC1,KC2,MCAU,IDLSTBM,KOCCC1,ISYMBBD 80 INTEGER KVIRDC1,ISYMBBB,KVIRBC1,KOFFOCC,KOFFINTD,KOFFINTB 81C 82 INTEGER ISYMLSTB 83C 84Cfunctions 85 INTEGER ILRCAMP 86 87C 88 REAL*8 XI1(*), XI2(*), XI1EFF(*), XI2EFF(*) 89 REAL*8 FOCKY(*), FOCK0(*), XLAMDP(*), XLAMDH(*) 90 REAL*8 XINT, XTROC0, XTRVI0, XTRVI2, XDIA, XSMAT, XT2TP 91 REAL*8 XTROC02, XUMAT, FREQY, DDOT, XIAJB 92 REAL*8 XXI1EFF, XXI2EFF, FREQB, TWO 93 REAL*8 WORK(LWORK) 94 CHARACTER*10 MODEL 95 CHARACTER FNCKJD*8, FNDKBC*8, FNDELD*8 96 CHARACTER FNTOC*8, FN3VI2*8, LABELB*8 97 CHARACTER FN3VI*6 98 CHARACTER CDUMMY*1 99 CHARACTER*16 FNCKJDR,FNDKBCR 100 CHARACTER*13 FN3SRTR,FNDELDR 101C 102 PARAMETER(FN3SRTR = 'CC3_XISD_TMP1', FNDELDR = 'CC3_XISD_TMP2') 103C 104#include "priunit.h" 105#include "dummy.h" 106#include "iratdef.h" 107#include "ccsdsym.h" 108#include "inftap.h" 109#include "ccinftap.h" 110#include "ccorb.h" 111#include "ccsdinp.h" 112#include "ccr1rsp.h" 113#include "ccrc1rsp.h" 114C 115 PARAMETER(TWO = 2.0D0, CDUMMY = ' ') 116C 117 CALL QENTER('CC3_XISD') 118C 119 !Check value of ICAU 120 IF (LCAUCHY .AND. ICAU.LT.-1) THEN 121 WRITE(LUPRI,*)'ICAU = ', ICAU 122 WRITE(LUPRI,*)'No CC3 cauchy implementation for ICAU.LT.-1' 123 CALL QUIT('Illegal value of ICAU in CC3_XISD') 124 END IF 125C 126 !Set variables associated with cauchy vectors 127 IF (LCAUCHY .AND. ICAU.GE.0) THEN 128 ISYMLSTB = 1 129 IDLSTB = ILRCAMP(LABELB,ICAU,ISYMLSTB) 130 ISYMBB = ISYLRC(IDLSTB) 131c FREQB = ZERO 132 LABELB = LRCLBL(IDLSTB) 133 NCAU = ILRCAU(IDLSTB) 134C 135 !Check symmetry 136 IF (ISYMBB .NE. ISYFKY) THEN 137 WRITE(LUPRI,*)'Symmetry of the perturbation operator: ',ISYFKY 138 WRITE(LUPRI,*)'Symmetry of the Cauchy vector: ',ISYMBB 139 CALL QUIT('Symmetry inconsistency in CC3_XISD (Cauchy vec)') 140 END IF 141C 142 END IF !LCAUCHY 143C 144C-------------------------- 145C Save and set flags. 146C-------------------------- 147C 148C Set symmetry flags. 149C 150C 151C ISYMT2 is symmetry of T2TP 152C ISYFKY is symmetry of perturbation operator 153C ISINT2 is symmetry of integrals in triples equation (ISINT2=1) 154C ISINT1 is symmetry of integrals in contraction (ISINT1=1) 155C ISYMT3 is symmetry of S and U intermediate 156C ISYMW is symmetry of W intermediate 157C ISYRES is symmetry of result vector (here the same as ISYFKY) 158C ISYMW = ISYFKY*ISYMT3 159C ISYRES = ISYMT3*ISYFKY*ISINT1 160C 161C------------------------------------------------------------- 162C 163 IPRCC = IPRINT 164 ISINT1 = 1 165 ISINT2 = 1 166 ISYMT2 = 1 167 ISYMT3 = MULD2H(ISINT2,ISYMT2) 168 ISYMW = MULD2H(ISYMT3,ISYFKY) 169 ISYRES = MULD2H(ISYMW,ISINT1) 170C 171C--------------------------------------------------------- 172C Work allocation 173C--------------------------------------------------------- 174C 175 KT2AM = 1 176 KRBJIA = KT2AM + NT2SQ(ISYMT2) 177 KFOCKD = KRBJIA + NT2SQ(ISYRES) 178 KFCKBA = KFOCKD + NORBTS 179 KFCKY = KFCKBA + N2BST(ISINT1) 180 KEND0 = KFCKY + N2BST(ISYFKY) 181 LWRK0 = LWORK - KEND0 182C 183 IF (LCAUCHY) THEN 184 KC1 = KEND0 185 KC2 = KC1 + (NCAU+1)*NT1AM(ISYMBB) 186 KEND0 = KC2 + (NCAU+1)*NT2SQ(ISYMBB) 187 LWRK0 = LWORK - KEND0 188 END IF 189C 190 IF (LWRK0.LT.0) THEN 191 WRITE(LUPRI,*)'Memory available :',LWORK 192 WRITE(LUPRI,*)'Memory needed :',KEND0 193 CALL QUIT('Insufficient memory in CC3_XISD(c1)') 194 END IF 195C 196 CALL DZERO(WORK(KRBJIA),NT2SQ(ISYRES)) 197 CALL DZERO(XI1EFF,NT1AM(ISYRES)) 198C 199C-------------------------------------- 200C Read in t2 , square it and reorder 201C-------------------------------------- 202C 203 IOPT = 2 204 CALL CC_RDRSP('R0',0,ISYMT2,IOPT,MODEL,DUMMY,WORK(KEND0)) 205 CALL CC_T2SQ(WORK(KEND0),WORK(KT2AM),ISYMT2) 206 IF (LWORK .LT. NT2SQ(ISYMT2)) THEN 207 CALL QUIT('Not enough memory to construct T2TP in CC3_XISD') 208 ENDIF 209C 210 CALL DCOPY(NT2SQ(ISYMT2),WORK(KT2AM),1,WORK(KEND0),1) 211 CALL CC3_T2TP(WORK(KT2AM),WORK(KEND0),ISYMT2) 212C 213 IF (IPRINT .GT. 55) THEN 214 XT2TP = DDOT(NT2SQ(ISYMT2),WORK(KT2AM),1,WORK(KT2AM),1) 215 WRITE(LUPRI,*) 'Norm of T2TP ',XT2TP 216 ENDIF 217C 218C 219 IF (LCAUCHY .AND. ICAU.GE.0) THEN 220C 221 IF (LWRK0.LT.NT2SQ(ISYMBB)) THEN 222 WRITE(LUPRI,*)'Memory available :',LWRK0 223 WRITE(LUPRI,*)'Memory needed :',NT2SQ(ISYMBB) 224 CALL QUIT('Insufficient memory in CC3_XISD(c2)') 225 END IF 226C 227C-------------------------------------------------------------------------- 228C Loop over cauchy vectors in order to 229C 1) get C1 and C2 cauchy vectors 230C 2) get C1-transformed integrals needed in <mu3|[[H^0,C1],T2^0]|HF> 231C-------------------------------------------------------------------------- 232C 233 DO MCAU = 0, NCAU 234C 235 !Open temporary files 236 LU3SRTR = -1 237 LUDELDR = -1 238 CALL WOPEN2(LU3SRTR,FN3SRTR,64,0) 239 CALL WOPEN2(LUDELDR,FNDELDR,64,0) 240C 241 !Get the list for current MCAU 242 IDLSTBM = ILRCAMP(LABELB,MCAU,ISYMBB) 243C 244 !Consistency check 245 IF (MCAU.EQ.NCAU .AND. IDLSTBM.NE.IDLSTB) THEN 246 CALL QUIT('Inconsistency in Cauchy loop in CC3_XISD') 247 END IF 248C 249 !Read in C1 and C2 Cauchy vectors for each MCAU 250 IOPT = 3 251 KOFF1 = KC1 + MCAU*NT1AM(ISYMBB) 252 KOFF2 = KC2 + MCAU*NT2SQ(ISYMBB) 253 CALL CC_RDRSP('RC ',IDLSTBM,ISYMBB,IOPT,MODEL,WORK(KOFF1), 254 * WORK(KOFF2)) 255 CALL CCLR_DIASCL(WORK(KOFF2),TWO,ISYMBB) 256 CALL CC_T2SQ(WORK(KOFF2),WORK(KEND0),ISYMBB) 257 CALL CC3_T2TP(WORK(KOFF2),WORK(KEND0),ISYMBB) 258C 259 !Generate file names for C1-transformed integrals 260 CALL GET_FILENAME(FNCKJDR,FNDKBCR,MCAU) 261C 262 !Open the files for C1-transformed integrals 263 LUCKJDR = -1 264 LUDKBCR = -1 265 CALL WOPEN2(LUCKJDR,FNCKJDR,64,0) 266 CALL WOPEN2(LUDKBCR,FNDKBCR,64,0) 267C 268 !Construct C1-transformed integrals and store on file 269 CALL CC3_BARINT(WORK(KOFF1),ISYMBB,XLAMDP, 270 * XLAMDH,WORK(KEND0),LWRK0, 271 * LU3SRTR,FN3SRTR,LUCKJDR,FNCKJDR) 272C 273 CALL CC3_SORT1(WORK(KEND0),LWRK0,2,ISYMBB,LU3SRTR,FN3SRTR, 274 * LUDELDR,FNDELDR,IDUMMY,CDUMMY,IDUMMY,CDUMMY, 275 * IDUMMY,CDUMMY) 276C 277 CALL CC3_SINT(XLAMDH,WORK(KEND0),LWRK0,ISYMBB, 278 * LUDELDR,FNDELDR,LUDKBCR,FNDKBCR) 279 280C 281 !Close the files for C1-transformed integrals 282 CALL WCLOSE2(LUCKJDR,FNCKJDR,'KEEP') 283 CALL WCLOSE2(LUDKBCR,FNDKBCR,'KEEP') 284C 285 !Close and delete temporary files 286 CALL WCLOSE2(LU3SRTR,FN3SRTR,'DELETE') 287 CALL WCLOSE2(LUDELDR,FNDELDR,'DELETE') 288C 289 END DO !MCAU 290C 291 END IF!LCAUCHY 292 293C 294C--------------------------------------------------------- 295C Read canonical orbital energies. 296C--------------------------------------------------------- 297C 298C-------------------------------------- 299C Read in orbital energies 300C-------------------------------------- 301C 302 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY, 303 & .FALSE.) 304 REWIND LUSIFC 305C 306 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 307 READ (LUSIFC) 308 READ (LUSIFC) (WORK(KFOCKD+I-1), I=1,NORBTS) 309C 310 CALL GPCLOSE(LUSIFC,'KEEP') 311C 312C--------------------------------------------- 313C Delete frozen orbitals in Fock diagonal. 314C--------------------------------------------- 315C 316 IF (FROIMP .OR. FROEXP) 317 * CALL CCSD_DELFRO(WORK(KFOCKD),WORK(KEND0),LWRK0) 318C 319C If we want to sum the T3 amplitudes 320 if (.false.) then 321 kx3am = kend0 322 kend0 = kx3am + nt1amx*nt1amx*nt1amx 323 call dzero(work(kx3am),nt1amx*nt1amx*nt1amx) 324 lwrk0 = lwork - kend0 325 endif 326C 327 IF (LWRK0 .LT. 0) THEN 328 WRITE(LUPRI,*) 'Memory available : ',LWORK 329 WRITE(LUPRI,*) 'Memory needed : ',KEND0 330 CALL QUIT('Insufficient space in CC3_XISD') 331 END IF 332C 333C----------------------------------------------------- 334C Construct the T1 transformed Fock matrix 335C----------------------------------------------------- 336C 337 LUFCK = -1 338C This AO Fock matrix is constructed from the T1 transformed density 339 CALL GPOPEN(LUFCK,'CC_FCKH','UNKNOWN',' ','UNFORMATTED', 340 * IDUMMY,.FALSE.) 341 REWIND(LUFCK) 342 READ(LUFCK)(WORK(KFCKBA + I-1),I = 1,N2BST(ISINT1)) 343 CALL GPCLOSE(LUFCK,'KEEP' ) 344C 345 IF (IPRINT .GT. 140) THEN 346 CALL AROUND( 'Usual Fock AO matrix' ) 347 CALL CC_PRFCKAO(WORK(KFCKBA),1) 348 ENDIF 349C 350 ! SCF Fock matrix in transformed using Lamda vector 351 CALL CC_FCKMO(WORK(KFCKBA),XLAMDP,XLAMDH, 352 * WORK(KEND0),LWRK0,1,1,1) 353C 354 IF (IPRINT .GT. 50) THEN 355 CALL AROUND( 'In CC3_XISD: T1 transformed Fock matrix' ) 356 CALL CC_PRFCKMO(WORK(KFCKBA),ISINT1) 357 ENDIF 358C 359C Sort the fock matrix 360C 361 CALL DCOPY(N2BST(ISINT1),WORK(KFCKBA),1,WORK(KEND0),1) 362C 363 DO ISYMC = 1,NSYM 364C 365 ISYMK = MULD2H(ISYMC,ISINT1) 366C 367 DO K = 1,NRHF(ISYMK) 368C 369 DO C = 1,NVIR(ISYMC) 370C 371 KOFF1 = KEND0 + IFCVIR(ISYMK,ISYMC) + 372 * NORB(ISYMK)*(C - 1) + K - 1 373 KOFF2 = KFCKBA + IT1AM(ISYMC,ISYMK) 374 * + NVIR(ISYMC)*(K - 1) + C - 1 375C 376 WORK(KOFF2) = WORK(KOFF1) 377C 378 ENDDO 379 ENDDO 380 ENDDO 381C 382 IF (IPRINT .GT. 50) THEN 383 CALL AROUND('In CC3_XISD: T1 transf Fock matrix (sort(ck))') 384 CALL CC_PRFCKMO(WORK(KFCKBA),1) 385 ENDIF 386C 387C---------------------------------------------- 388C Sort FOCKY to (ck) 389C---------------------------------------------- 390C 391 CALL DCOPY(N2BST(ISYFKY),FOCKY,1,WORK(KEND0),1) 392C 393 DO ISYMC = 1,NSYM 394C 395 ISYMK = MULD2H(ISYMC,ISYFKY) 396C 397 DO K = 1,NRHF(ISYMK) 398C 399 DO C = 1,NVIR(ISYMC) 400C 401 KOFF1 = KEND0 + IFCVIR(ISYMK,ISYMC) + 402 * NORB(ISYMK)*(C - 1) + K - 1 403 KOFF2 = KFCKY + IT1AM(ISYMC,ISYMK) 404 * + NVIR(ISYMC)*(K - 1) + C - 1 405C 406 WORK(KOFF2) = WORK(KOFF1) 407C 408 ENDDO 409 ENDDO 410 ENDDO 411C 412 IF (IPRINT .GT. 50) THEN 413 CALL AROUND('In CC3_XISD: T1^Y transf Fock matrix (sort(ck))') 414 CALL CC_PRFCKMO(WORK(KFCKY),1) 415 ENDIF 416C 417C----------------------------- 418C Read occupied integrals. 419C----------------------------- 420C 421C Memory allocation. 422C 423 KTROC = KEND0 424 KTROC1 = KTROC + NTRAOC(ISINT1) 425 KTROC0 = KTROC1 + NTRAOC(ISINT1) 426 KTROC02= KTROC0 + NTRAOC(ISINT2) 427 KXIAJB = KTROC02+ NTRAOC(ISINT2) 428 KEND1 = KXIAJB + NT2AM(ISINT1) 429 LWRK1 = LWORK - KEND1 430C 431 KINTOC = KEND1 432 KEND2 = KINTOC + MAX(NTOTOC(ISINT2),NTOTOC(ISINT1)) 433 LWRK2 = LWORK - KEND2 434C 435 IF (LWRK2 .LT. 0) THEN 436 WRITE(LUPRI,*) 'Memory available : ',LWORK 437 WRITE(LUPRI,*) 'Memory needed : ',KEND2 438 CALL QUIT('Insufficient space in CC3_XISD') 439 END IF 440C 441C------------------------ 442C Construct L(ia,jb). 443C------------------------ 444C 445 LENGTH = IRAT*NT2AM(ISINT1) 446C 447 REWIND(LUIAJB) 448 CALL READI(LUIAJB,LENGTH,WORK(KXIAJB)) 449 ISYOPE = ISINT1 450 IOPTTCME = 1 451 CALL CCSD_TCMEPK(WORK(KXIAJB),1.0D0,ISYOPE,IOPTTCME) 452C 453C 454 IF ( IPRINT .GT. 55) THEN 455 XIAJB = DDOT(NT2AM(ISINT1),WORK(KXIAJB),1, 456 * WORK(KXIAJB),1) 457 WRITE(LUPRI,*) 'Norm of IAJB ',XIAJB 458 ENDIF 459C 460C------------------------ 461C Occupied integrals. 462C 463C Read in integrals for SMAT etc. 464C----------------------- 465C 466 IOFF = 1 467 IF (NTOTOC(ISINT2) .GT. 0) THEN 468 CALL GETWA2(LUCKJD,FNCKJD,WORK(KINTOC),IOFF,NTOTOC(ISINT2)) 469 ENDIF 470C 471C---------------------------------- 472C Write out norms of Integrals. 473C---------------------------------- 474C 475 IF (IPRINT .GT. 55) THEN 476 XINT = DDOT(NTOTOC(ISINT2),WORK(KINTOC),1, 477 * WORK(KINTOC),1) 478 WRITE(LUPRI,*) 'Norm of CCSDT_OC-INT ',XINT 479 ENDIF 480C 481C---------------------------------------------------------------------- 482C Transform (ai|j delta) integrals to (ai|j k) and sort as (ij,k,a) 483C---------------------------------------------------------------------- 484C here the xlamdp transformation need to be used 485C (delta is particle index) 486C 487 CALL CC3_TROCC(WORK(KINTOC),WORK(KTROC0),XLAMDP, 488 * WORK(KEND2),LWRK2,ISINT2) 489 490C 491C---------------------------------------------------------------------- 492C (ai|j k) sorted as (ij,k,a) 493C---------------------------------------------------------------------- 494C 495 CALL CCFOP_SORT(WORK(KTROC0),WORK(KTROC02),ISINT2,1) 496 497C 498C------------------------------- 499C Write out norms of arrays. 500C------------------------------- 501C 502C 503 IF (IPRINT .GT. 55) THEN 504 XINT = DDOT(NTOTOC(ISINT2),WORK(KINTOC),1, 505 * WORK(KINTOC),1) 506 WRITE(LUPRI,*) 'Norm of CKJDEL-INT ',XINT 507 ENDIF 508C 509 IF (IPRINT .GT. 55) THEN 510 XTROC0 = DDOT(NTRAOC(ISINT2),WORK(KTROC0),1, 511 * WORK(KTROC0),1) 512 WRITE(LUPRI,*) 'Norm of TROC0 ',XTROC0 513 ENDIF 514C 515 IF (IPRINT .GT. 55) THEN 516 XTROC02 = DDOT(NTRAOC(ISINT2),WORK(KTROC02),1, 517 * WORK(KTROC02),1) 518 WRITE(LUPRI,*) 'Norm of TROC02 ',XTROC02 519 ENDIF 520C 521C 522C---------------------------------- 523C 524C------------------------ 525C Occupied integrals. 526C 527C Read in integrals for WMAT etc. 528C----------------------- 529C 530C 531 IOFF = 1 532 IF (NTOTOC(ISYMOP) .GT. 0) THEN 533 CALL GETWA2(LUTOC,FNTOC,WORK(KINTOC),IOFF,NTOTOC(ISINT1)) 534 ENDIF 535C 536C---------------------------------- 537C Write out norms of Integrals. 538C---------------------------------- 539C 540 IF (IPRINT .GT. 55) THEN 541 XINT = DDOT(NTOTOC(ISYMOP),WORK(KINTOC),1, 542 * WORK(KINTOC),1) 543 WRITE(LUPRI,*) 'Norm of CCSDT_OC-INT ',XINT 544 ENDIF 545C 546C 547C---------------------------------------------------------------------- 548C Transform (ia|j delta) integrals to (ia|j k) and sort as (i,j,k,a) 549C---------------------------------------------------------------------- 550C 551 CALL CCSDT_TROCC(WORK(KINTOC),WORK(KTROC),XLAMDH, 552 * WORK(KEND2),LWRK2) 553C 554C---------------------------------------------------------------------- 555C sort (i,j,k,a) as (a,i,j,k) 556C---------------------------------------------------------------------- 557C 558 559 CALL CCSDT_SRTOC2(WORK(KTROC),WORK(KTROC1),ISINT1, 560 * WORK(KEND2),LWRK2) 561C 562 IF (LCAUCHY .AND. ICAU.GE.0) THEN 563C 564 !allocate space for ALL C1-transformed occupied integrals 565 KOCCC1 = KEND1 566 KEND1 = KOCCC1 + (NCAU+1)*NTRAOC(ISYMBB) 567 LWRK1 = LWORK - KEND1 568C 569 IF (LWRK1 .LT. 0) THEN 570 WRITE(LUPRI,*) 'Memory available : ',LWORK 571 WRITE(LUPRI,*) 'Memory needed : ',KEND1 572 CALL QUIT('Insufficient space in CC3_XISD(c3)') 573 END IF 574C 575 DO MCAU = 0,NCAU 576C 577C-------------------------------------------------------------- 578C Read in C1-transformed occupied integrals for each MCAU 579C------------------------------------------------------------- 580C 581 !Generate file name for C1-transformed integrals 582 !(FNDKBCR is not needed here) 583 CALL GET_FILENAME(FNCKJDR,FNDKBCR,MCAU) 584C 585 !Open the file for C1-transformed integrals 586 LUCKJDR = -1 587 CALL WOPEN2(LUCKJDR,FNCKJDR,64,0) 588C 589 !Get the offset for a given MCAU 590 KOFF1 = KOCCC1 + MCAU*NTRAOC(ISYMBB) 591C 592 !Read in the integrals 593 CALL INTOCC_T3X(LUCKJDR,FNCKJDR,XLAMDP,ISYMBB, 594 * WORK(KOFF1),WORK(KEND1),LWRK1) 595C 596 !Close and delete the file for C1-transformed occupied integrals 597 CALL WCLOSE2(LUCKJDR,FNCKJDR,'DELETE') 598C 599 END DO !MCAU 600C 601 END IF !LCAUCHY 602C 603C 604C---------------------------- 605C General loop structure. 606C---------------------------- 607C 608 DO ISYMD = 1,NSYM 609C 610 ISAIJ1 = MULD2H(ISYMD,ISYRES) 611 ISCKB2 = MULD2H(ISINT2,ISYMD) 612 ISCKB1 = MULD2H(ISINT1,ISYMD) 613C 614 IF (IPRINT .GT. 55) THEN 615C 616 WRITE(LUPRI,*) 'In CC3_XI: ISCKB2 :',ISCKB2 617C 618 ENDIF 619 620C 621C-------------------------- 622C Memory allocation. 623C-------------------------- 624C 625 IF (LCAUCHY .AND. ICAU.GE.0) THEN 626C 627 !get symmetry of C1-transformed virtual integrals with fixed D 628 ISYMBBD = MULD2H(ISYMBB,ISYMD) 629C 630 !allocation for ALL C1-transformed virtual integrals with fixed D 631 KVIRDC1 = KEND1 632 KEND1 = KVIRDC1 + (NCAU+1)*NCKATR(ISYMBBD) 633 LWRK1 = LWORK - KEND1 634C 635 IF (LWRK1 .LT. 0) THEN 636 WRITE(LUPRI,*) 'Memory available : ',LWORK 637 WRITE(LUPRI,*) 'Memory needed : ',KEND1 638 CALL QUIT('Insufficient space in CC3_XISD(c4)') 639 END IF 640C 641 END IF !LCAUCHY 642C 643 KTRVI = KEND1 644 KTRVI1 = KTRVI + NCKATR(ISCKB1) 645 KTRVI3 = KTRVI1 + NCKATR(ISCKB1) 646 KTRVI2 = KTRVI3 + NCKATR(ISCKB2) 647 KRMAT1 = KTRVI2 + NCKATR(ISCKB2) 648 KEND2 = KRMAT1 + NCKI(ISAIJ1) 649 LWRK2 = LWORK - KEND2 650C 651 KTRVI0 = KEND2 652 KEND3 = KTRVI0 + NCKATR(ISCKB2) 653 LWRK3 = LWORK - KEND3 654C 655 KINTVI = KEND3 656 KEND4 = KINTVI + MAX(NCKA(ISCKB2),NCKA(ISCKB1)) 657 LWRK4 = LWORK - KEND4 658C 659 IF (LWRK4 .LT. 0) THEN 660 WRITE(LUPRI,*) 'Memory available : ',LWORK 661 WRITE(LUPRI,*) 'Memory needed : ',KEND4 662 CALL QUIT('Insufficient space in CC3_XISD') 663 END IF 664C 665C--------------------- 666C Sum over D 667C--------------------- 668C 669 DO D = 1,NVIR(ISYMD) 670C 671C--------------------------------------- 672C Initialise 673C--------------------------------------- 674C 675 CALL DZERO(WORK(KRMAT1),NCKI(ISAIJ1)) 676C 677C------------------------------------------------------------ 678C Read and transform integrals used in contraction 679C with W intermediate. 680C------------------------------------------------------------ 681C 682 IOFF = ICKAD(ISCKB1,ISYMD) + NCKA(ISCKB1)*(D - 1) + 1 683 IF (NCKA(ISCKB1) .GT. 0) THEN 684 CALL GETWA2(LU3VI2,FN3VI2,WORK(KINTVI),IOFF, 685 & NCKA(ISCKB1)) 686 ENDIF 687 CALL CCSDT_TRVIR(WORK(KINTVI),WORK(KTRVI),XLAMDP, 688 * ISYMD,D,ISYMOP,WORK(KEND4),LWRK4) 689C 690 IOFF = ICKAD(ISCKB1,ISYMD) + NCKA(ISCKB1)*(D - 1) + 1 691 IF (NCKA(ISCKB1) .GT. 0) THEN 692 CALL GETWA2(LU3VI,FN3VI,WORK(KINTVI),IOFF, 693 & NCKA(ISCKB1)) 694 ENDIF 695 CALL CCSDT_TRVIR(WORK(KINTVI),WORK(KTRVI1),XLAMDP, 696 * ISYMD,D,ISYMOP,WORK(KEND4),LWRK4) 697C 698C----------------------------------------------- 699C Read virtual integrals used in first s3am. 700C----------------------------------------------- 701C 702 IOFF = ICKBD(ISCKB2,ISYMD) + NCKATR(ISCKB2)*(D - 1) + 1 703 IF (NCKATR(ISCKB2) .GT. 0) THEN 704 CALL GETWA2(LUDKBC,FNDKBC,WORK(KTRVI0),IOFF, 705 & NCKATR(ISCKB2)) 706 ENDIF 707C 708C----------------------------------------------------------- 709C Sort the integrals for s3am 710C----------------------------------------------------------- 711C 712 CALL CCSDT_SRTVIR(WORK(KTRVI0),WORK(KTRVI2),WORK(KEND4), 713 * LWRK4,ISYMD,ISINT2) 714C 715C 716 IF (IPRINT .GT. 55) THEN 717 XTRVI0= DDOT(NCKATR(ISCKB2),WORK(KTRVI0),1, 718 * WORK(KTRVI0),1) 719 WRITE(LUPRI,*) 'Norm of TRVI0 ',XTRVI0 720 ENDIF 721C 722 IF (IPRINT .GT. 55) THEN 723 XTRVI2= DDOT(NCKATR(ISCKB2),WORK(KTRVI2),1, 724 * WORK(KTRVI2),1) 725 WRITE(LUPRI,*) 'Norm of TRVI2 ',XTRVI2 726 ENDIF 727 728 729C 730C------------------------------------------------------ 731C Read virtual integrals used in first U. 732C------------------------------------------------------ 733C 734 IOFF = ICKAD(ISCKB2,ISYMD) + NCKA(ISCKB2)*(D - 1) + 1 735 IF (NCKA(ISCKB2) .GT. 0) THEN 736 CALL GETWA2(LUDELD,FNDELD,WORK(KINTVI),IOFF, 737 & NCKA(ISCKB2)) 738 739 ENDIF 740C 741 CALL CCSDT_TRVIR(WORK(KINTVI),WORK(KTRVI3),XLAMDH, 742 * ISYMD,D,ISINT2,WORK(KEND4),LWRK4) 743C 744 745 IF (LCAUCHY .AND. ICAU.GE.0) THEN 746C 747 DO MCAU = 0, NCAU 748C 749C--------------------------------------------------------------------------- 750C Read in C1-transformed virt int with fixed D for each MCAU 751C--------------------------------------------------------------------------- 752C 753 !Generate file names for C1-transformed integrals 754 !(FNCKJDR is not needed here) 755 CALL GET_FILENAME(FNCKJDR,FNDKBCR,MCAU) 756C 757 !Open the files for C1-transformed integrals 758 LUDKBCR = -1 759 CALL WOPEN2(LUDKBCR,FNDKBCR,64,0) 760C 761 !Get the offset for a given MCAU 762 KOFF1 = KVIRDC1 + MCAU*NCKATR(ISYMBBD) 763C 764 !Read in the integrals 765 IOFF = ICKBD(ISYMBBD,ISYMD)+NCKATR(ISYMBBD)*(D-1)+1 766 IF (NCKATR(ISYMBBD) .GT. 0) THEN 767 CALL GETWA2(LUDKBCR,FNDKBCR,WORK(KOFF1),IOFF, 768 & NCKATR(ISYMBBD)) 769 ENDIF 770C 771 !Close the file for C1-transformed virtual integrals 772 !(and keep it: it will be needed in B-loop) 773 CALL WCLOSE2(LUDKBCR,FNDKBCR,'KEEP') 774C 775 END DO !MCAU 776C 777 END IF!LCAUCHY 778C 779C-------------------------------------------------------- 780C 781 DO ISYMB = 1,NSYM 782C 783 ISAIJ2 = MULD2H(ISYMB,ISYRES) 784 ISYALJ = MULD2H(ISYMB,ISYMT2) 785 ISYALJ2 = MULD2H(ISYMD,ISYMT2) 786 ISYMBD = MULD2H(ISYMB,ISYMD) 787 ISCKIJ = MULD2H(ISYMBD,ISYMT3) 788 ISCKD2 = MULD2H(ISINT2,ISYMB) 789 ISWMAT = MULD2H(ISYFKY,ISCKIJ) 790C 791 IF (IPRINT .GT. 55) THEN 792C 793 WRITE(LUPRI,*) 'In CC3_XISD : ISYMD :',ISYMD 794 WRITE(LUPRI,*) 'In CC3_XISD : ISYMB :',ISYMB 795 WRITE(LUPRI,*) 'In CC3_XISD : ISYALJ:',ISYALJ 796 WRITE(LUPRI,*) 'In CC3_XISD : ISYMBD:',ISYMBD 797 WRITE(LUPRI,*) 'In CC3_XISD : ISCKIJ:',ISCKIJ 798 WRITE(LUPRI,*) 'In CC3_XISD : ISWMAT:',ISWMAT 799C 800 ENDIF 801C 802 !We can use KEND3 since we do not need KINTVI array from D-loop 803 !anymore 804C 805 IF (LCAUCHY .AND. ICAU.GE.0) THEN 806C 807 !get symmetry of C1-transformed virt int with fixed B 808 ISYMBBB = MULD2H(ISYMBB,ISYMB) 809C 810 !allocation for ALL C1-transformed virt int with fixed B 811 KVIRBC1 = KEND3 812 KEND3 = KVIRBC1 + (NCAU+1)*NCKATR(ISYMBBB) 813 LWRK3 = LWORK - KEND3 814C 815 IF (LWRK3 .LT. 0) THEN 816 WRITE(LUPRI,*) 'Memory available : ',LWORK 817 WRITE(LUPRI,*) 'Memory needed : ',KEND3 818 CALL QUIT('Insufficient space in CC3_XISD(c5)') 819 END IF 820C 821 END IF !LCAUCHY 822C 823 KSMAT = KEND3 824 KSMAT3 = KSMAT + NCKIJ(ISCKIJ) 825 KUMAT = KSMAT3 + NCKIJ(ISCKIJ) 826 KUMAT3 = KUMAT + NCKIJ(ISCKIJ) 827 KDIAG = KUMAT3 + NCKIJ(ISCKIJ) 828 KDIAGW = KDIAG + NCKIJ(ISCKIJ) 829 KWMAT = KDIAGW + NCKIJ(ISWMAT) 830 KRMAT2 = KWMAT + NCKIJ(ISWMAT) 831 KINDSQW = KRMAT2 + NCKI(ISAIJ2) 832 KINDSQ = KINDSQW + (6*NCKIJ(ISWMAT) - 1)/IRAT + 1 833 KINDEX = KINDSQ + (6*NCKIJ(ISCKIJ) - 1)/IRAT + 1 834 KINDEX2 = KINDEX + (NCKI(ISYALJ) - 1)/IRAT + 1 835 KTMAT = KINDEX2 + (NCKI(ISYALJ2) - 1)/IRAT + 1 836 KTRVI8 = KTMAT + MAX(NCKIJ(ISCKIJ),NCKIJ(ISWMAT)) 837 KTRVI9 = KTRVI8 + NCKATR(ISCKD2) 838 KTRVI10 = KTRVI9 + NCKATR(ISCKD2) 839 KEND4 = KTRVI10 + NCKATR(ISCKD2) 840 LWRK4 = LWORK - KEND4 841C 842 KINTVI = KEND4 843 KEND5 = KINTVI + NCKA(ISCKD2) 844 LWRK5 = LWORK - KEND5 845C 846 IF (LWRK5 .LT. 0) THEN 847 WRITE(LUPRI,*) 'Memory available : ',LWORK 848 WRITE(LUPRI,*) 'Memory needed : ',KEND5 849 CALL QUIT('Insufficient space in CC3_XISD ') 850 END IF 851C 852C--------------------------------------------- 853C Construct part of the diagonal. 854C--------------------------------------------- 855C 856 CALL CC3_DIAG(WORK(KDIAG),WORK(KFOCKD),ISCKIJ) 857 CALL CC3_DIAG(WORK(KDIAGW),WORK(KFOCKD),ISWMAT) 858C 859 IF (IPRINT .GT. 55) THEN 860 XDIA = DDOT(NCKIJ(ISCKIJ),WORK(KDIAG),1, 861 * WORK(KDIAG),1) 862 WRITE(LUPRI,*) 'Norm of DIA ',XDIA 863 XDIA = DDOT(NCKIJ(ISCKIJ),WORK(KDIAGW),1, 864 * WORK(KDIAGW),1) 865 WRITE(LUPRI,*) 'Norm of DIA_W',XDIA 866 ENDIF 867C 868C------------------------------------- 869C Construct index arrays. 870C------------------------------------- 871C 872 LENSQ = NCKIJ(ISCKIJ) 873 CALL CC3_INDSQ(WORK(KINDSQ),LENSQ,ISCKIJ) 874 CALL CC3_INDEX(WORK(KINDEX),ISYALJ) 875 CALL CC3_INDEX(WORK(KINDEX2),ISYALJ2) 876 LENSQW = NCKIJ(ISWMAT) 877 CALL CC3_INDSQ(WORK(KINDSQW),LENSQW,ISWMAT) 878C 879 DO B = 1,NVIR(ISYMB) 880C 881C--------------------------------------- 882C Initialise 883C--------------------------------------- 884C 885 CALL DZERO(WORK(KRMAT2),NCKI(ISAIJ2)) 886 CALL DZERO(WORK(KWMAT),NCKIJ(ISWMAT)) 887C 888C----------------------------------------------- 889C Read virtual integrals used in second s3am. 890C----------------------------------------------- 891C 892 IOFF = ICKBD(ISCKD2,ISYMB) + 893 & NCKATR(ISCKD2)*(B - 1) + 1 894 IF (NCKATR(ISCKD2) .GT. 0) THEN 895 CALL GETWA2(LUDKBC,FNDKBC,WORK(KTRVI8),IOFF, 896 & NCKATR(ISCKD2)) 897 ENDIF 898C 899C----------------------------------------------------------- 900C Sort the integrals for s3am 901C----------------------------------------------------------- 902C 903 CALL CCSDT_SRTVIR(WORK(KTRVI8), 904 * WORK(KTRVI9),WORK(KEND4), 905 * LWRK4,ISYMB,ISINT2) 906C 907C 908C---------------------------------------------------------- 909C Read virtual integrals used in second U 910C---------------------------------------------------------- 911C 912 IOFF = ICKAD(ISCKD2,ISYMB) + NCKA(ISCKD2)*(B - 1) + 1 913 IF (NCKA(ISCKD2) .GT. 0) THEN 914 CALL GETWA2(LUDELD,FNDELD,WORK(KINTVI),IOFF, 915 * NCKA(ISCKD2)) 916 ENDIF 917C 918 CALL CCSDT_TRVIR(WORK(KINTVI),WORK(KTRVI10),XLAMDH, 919 * ISYMB,B,ISINT2,WORK(KEND5),LWRK5) 920C 921C 922 IF (LCAUCHY .AND. ICAU.GE.0) THEN 923C 924 DO MCAU = 0, NCAU 925C 926C--------------------------------------------------------------------------- 927C Read in C1-transformed virt int with fixed B for each MCAU 928C--------------------------------------------------------------------------- 929C 930 !Generate file names for C1-transformed integrals 931 !(FNCKJDR is not needed here) 932 CALL GET_FILENAME(FNCKJDR,FNDKBCR,MCAU) 933C 934 !Open the files for C1-transformed integrals 935 LUDKBCR = -1 936 CALL WOPEN2(LUDKBCR,FNDKBCR,64,0) 937C 938 !Get the offset for a given MCAU 939 KOFF1 = KVIRBC1 + MCAU*NCKATR(ISYMBBB) 940C 941 !Read in the integrals 942 IOFF = ICKBD(ISYMBBB,ISYMB) 943 * +NCKATR(ISYMBBB)*(B-1)+1 944 IF (NCKATR(ISYMBBB) .GT. 0) THEN 945 CALL GETWA2(LUDKBCR,FNDKBCR,WORK(KOFF1),IOFF, 946 & NCKATR(ISYMBBB)) 947 ENDIF 948C 949 !Close and keep the file for C1-transformed virt int 950 !...or delte it if you don't need it anymore 951 952 IF ( (ISYMD.EQ.NSYM) .AND. (ISYMB.EQ.NSYM) 953 * .AND. (D.EQ.NVIR(ISYMD)) 954 * .AND. (B.EQ.NVIR(ISYMB))) THEN 955C 956 CALL WCLOSE2(LUDKBCR,FNDKBCR,'DELETE') 957 ELSE 958 CALL WCLOSE2(LUDKBCR,FNDKBCR,'KEEP') 959 END IF 960C 961 END DO !MCAU 962C 963 END IF!LCAUCHY 964 965C 966C------------------------------------------------------------------------ 967C Calculate the S(ci,bk,dj) matrix for T3 for B,D. 968C------------------------------------------------------------------- 969C 970 CALL CC3_SMAT(0.0D0,WORK(KT2AM),ISYMT2, 971 * WORK(KTMAT),WORK(KTRVI0), 972 * WORK(KTRVI2),WORK(KTROC0),ISINT2, 973 * WORK(KFOCKD),WORK(KDIAG), 974 * WORK(KSMAT),WORK(KEND4),LWRK4, 975 * WORK(KINDEX),WORK(KINDSQ),LENSQ, 976 * ISYMB,B,ISYMD,D) 977C 978 CALL T3_FORBIDDEN(WORK(KSMAT),ISYMT3,ISYMB,B,ISYMD,D) 979C 980c call sum_pt3(work(ksmat),isymb,b,isymd,d, 981c * isckij,work(kx3am),1) 982C 983 IF (IPRINT .GT. 55) THEN 984 XSMAT = DDOT(NCKIJ(ISCKIJ),WORK(KSMAT),1, 985 * WORK(KSMAT),1) 986 WRITE(LUPRI,*) 'Norm of SMAT ',XSMAT 987 ENDIF 988C 989C------------------------------------------------------------------- 990C Calculate the S(ci,bk,dj) matrix for T3 for D,B. 991C------------------------------------------------------------------- 992C 993 CALL CC3_SMAT(0.0D0,WORK(KT2AM),ISYMT2, 994 * WORK(KTMAT),WORK(KTRVI8), 995 * WORK(KTRVI9),WORK(KTROC0),ISINT2, 996 * WORK(KFOCKD),WORK(KDIAG), 997 * WORK(KSMAT3),WORK(KEND4),LWRK4, 998 * WORK(KINDEX2),WORK(KINDSQ),LENSQ, 999 * ISYMD,D,ISYMB,B) 1000C 1001 CALL T3_FORBIDDEN(WORK(KSMAT3),ISYMT3,ISYMD,D,ISYMB,B) 1002C 1003c call sum_pt3(work(ksmat3),isymd,d,isymb,b, 1004c * isckij,work(kx3am),1) 1005C 1006 IF (IPRINT .GT. 55) THEN 1007 XSMAT = DDOT(NCKIJ(ISCKIJ),WORK(KSMAT3),1, 1008 * WORK(KSMAT3),1) 1009 WRITE(LUPRI,*) 'Norm of SMAT3 ',XSMAT 1010 ENDIF 1011C 1012C--------------------------------------------------------------------------- 1013C Calculate U(ci,jk) for fixed b,d. 1014C-------------------------------------------------- 1015C 1016 CALL CC3_UMAT(0.0D0,WORK(KT2AM),ISYMT2,WORK(KTRVI3), 1017 * WORK(KTROC02),ISINT2,WORK(KFOCKD), 1018 * WORK(KDIAG),WORK(KUMAT),WORK(KTMAT), 1019 * WORK(KEND4),LWRK4,WORK(KINDSQ),LENSQ, 1020 * ISYMB,B,ISYMD,D) 1021C 1022 CALL T3_FORBIDDEN(WORK(KUMAT),ISYMT3,ISYMB,B,ISYMD,D) 1023C 1024c call sum_pt3(work(kumat),isymb,b,isymd,d, 1025c * isckij,work(kx3am),3) 1026C 1027 IF (IPRINT .GT. 55) THEN 1028 XUMAT = DDOT(NCKIJ(ISCKIJ),WORK(KUMAT),1, 1029 * WORK(KUMAT),1) 1030 WRITE(LUPRI,*) 'Norm of UMAT ',XUMAT 1031 ENDIF 1032C 1033C-------------------------------------------------- 1034C 1035C 1036C-------------------------------------------------- 1037C Calculate U(ci,jk) for fixed d,b. 1038C-------------------------------------------------- 1039C 1040 CALL CC3_UMAT(0.0D0,WORK(KT2AM),ISYMT2,WORK(KTRVI10), 1041 * WORK(KTROC02),ISINT2,WORK(KFOCKD), 1042 * WORK(KDIAG),WORK(KUMAT3),WORK(KTMAT), 1043 * WORK(KEND4),LWRK4,WORK(KINDSQ),LENSQ, 1044 * ISYMD,D,ISYMB,B) 1045C 1046 CALL T3_FORBIDDEN(WORK(KUMAT3),ISYMT3,ISYMD,D,ISYMB,B) 1047C 1048c call sum_pt3(work(kumat3),isymd,d,isymb,b, 1049c * isckij,work(kx3am),3) 1050C 1051 IF (IPRINT .GT. 55) THEN 1052 XUMAT = DDOT(NCKIJ(ISCKIJ),WORK(KUMAT3),1, 1053 * WORK(KUMAT3),1) 1054 WRITE(LUPRI,*) 'Norm of UMAT3 ',XUMAT 1055 ENDIF 1056C 1057C 1058C------------------------------------------------- 1059C Calculate the term <mu2|[V,T3]|HF> 1060C------------------------------------------------- 1061C 1062 IF (.NOT.(LCAUCHY .AND. (ICAU.GE.0))) THEN 1063 CALL CCFOP_ONEL(XI2EFF,WORK(KRMAT1),WORK(KRMAT2), 1064 * WORK(KFCKY),WORK(KSMAT),WORK(KTMAT), 1065 * ISYMT3,ISYFKY,WORK(KINDSQ), 1066 * LENSQ,WORK(KEND4),LWRK4, 1067 * ISYMB,B,ISYMD,D) 1068 END IF 1069C 1070C-------------------------------------------------- 1071C Sum up S and U intermediates to get T3 BD amplitudes 1072C-------------------------------------------------- 1073C 1074 CALL CC3_T3BD(ISCKIJ,WORK(KSMAT),WORK(KSMAT3), 1075 * WORK(KUMAT),WORK(KUMAT3),WORK(KTMAT), 1076 * WORK(KINDSQ),LENSQ) 1077c 1078c call sum_pt3(work(ktmat),isymb,b,isymd,d, 1079c * 1,work(kx3am),3) 1080 1081C 1082 IF (IPRINT .GT. 55) THEN 1083 XUMAT = DDOT(NCKIJ(ISCKIJ),WORK(KTMAT),1, 1084 * WORK(KTMAT),1) 1085 WRITE(LUPRI,*) 'Norm of TMAT^BD ',XUMAT 1086 ENDIF 1087C 1088C------------------------------------------------------ 1089C Based on T3 BD amplitudes calculate W BD intermediates 1090C------------------------------------------------------ 1091C 1092C------------------------------------------------------ 1093C Calculate the term <mu3|[Y,T3]|HF> virtual contribution 1094C added in W^BD(KWMAT) 1095C------------------------------------------------------ 1096C 1097 CALL WBD_V(WORK(KTMAT),ISCKIJ,FOCKY,ISYFKY, 1098 * WORK(KWMAT), 1099 * ISWMAT,WORK(KEND4),LWRK4) 1100C 1101C------------------------------------------------------ 1102C Calculate the term <mu3|[Y,T3]|HF> occupied contribution 1103C added in W^BD(KWMAT) 1104C------------------------------------------------------ 1105C 1106 CALL WBD_O(WORK(KTMAT),ISCKIJ,FOCKY,ISYFKY, 1107 * WORK(KWMAT), 1108 * ISWMAT,WORK(KEND4),LWRK4) 1109C 1110C------------------------------------------------------ 1111C Calculate the term <mu3|[[Y,T2],T2]|HF> 1112C added in W^BD(KWMAT) 1113C------------------------------------------------------ 1114C 1115 CALL WBD_T2(.FALSE.,B,ISYMB,D,ISYMD, 1116 * WORK(KT2AM),ISYMT2,WORK(KT2AM),ISYMT2, 1117 * FOCKY,ISYFKY, 1118 * WORK(KINDSQW),LENSQW,WORK(KWMAT), 1119 * ISWMAT,WORK(KEND4),LWRK4) 1120C 1121* call sum_pt3(work(kwmat),isymb,b,isymd,d, 1122* * iswmat,work(kx3am),4) 1123C 1124C 1125 IF (LCAUCHY .AND. ICAU.GE.0) THEN 1126C 1127 DO MCAU = 0, NCAU 1128C 1129 !Get offset for C2 vec for a given MCAU 1130 KOFF2 = KC2 + MCAU*NT2SQ(ISYMBB) 1131C 1132 !Calculate <mu3|[H^0,C2]|HF> cont to C3 cauchy vec 1133 CALL WBD_GROUND(WORK(KOFF2),ISYMBB,WORK(KTMAT), 1134 * WORK(KTRVI0),WORK(KTRVI8), 1135 * WORK(KTROC0),ISINT2,WORK(KWMAT), 1136 * WORK(KEND4),LWRK4, 1137 * WORK(KINDSQW),LENSQW, 1138 * ISYMB,B,ISYMD,D) 1139 1140C 1141 !Get offset for C1-trans occ int for a given MCAU 1142 KOFFOCC = KOCCC1 + MCAU*NTRAOC(ISYMBB) 1143C 1144 !Get offset for C1-trans virt int with fixed D for MCAU 1145 KOFFINTD = KVIRDC1 + MCAU*NCKATR(ISYMBBD) 1146 !Get offset for C1-trans virt int with fixed B for MCAU 1147 KOFFINTB = KVIRBC1 + MCAU*NCKATR(ISYMBBB) 1148C 1149 1150 !Calculate <mu3|[[H^0,C1],T2^0]|HF> 1151 CALL WBD_GROUND(WORK(KT2AM),ISYMT2,WORK(KTMAT), 1152 * WORK(KOFFINTD),WORK(KOFFINTB), 1153 * WORK(KOFFOCC),ISYMBB,WORK(KWMAT), 1154 * WORK(KEND4),LWRK4, 1155 * WORK(KINDSQW),LENSQW, 1156 * ISYMB,B,ISYMD,D) 1157C 1158* call sum_pt3(work(kwmat),isymb,b,isymd,d, 1159* * iswmat,work(kx3am),4) 1160 1161 !Divide by the energy difference 1162 !and remove the forbidden elements 1163 CALL WBD_DIA(B,ISYMB,D,ISYMD,FREQY, 1164 * ISWMAT,WORK(KWMAT), 1165 * WORK(KDIAGW),WORK(KFOCKD)) 1166C 1167 CALL T3_FORBIDDEN(WORK(KWMAT),ISYFKY,ISYMB,B, 1168 * ISYMD,D) 1169C 1170 END DO !MCAU 1171C 1172 !trun the sign C_T <- -C_T 1173 CALL DSCAL(NCKIJ(ISWMAT),-1.0D0,WORK(KWMAT),1) 1174C 1175 END IF !LCAUCHY 1176C 1177 1178* call sum_pt3(work(kwmat),isymb,b,isymd,d, 1179* * iswmat,work(kx3am),4) 1180 1181 !Divide by the energy difference and 1182 !remove the forbidden elements 1183 CALL WBD_DIA(B,ISYMB,D,ISYMD,FREQY, 1184 * ISWMAT,WORK(KWMAT), 1185 * WORK(KDIAGW),WORK(KFOCKD)) 1186C 1187 CALL T3_FORBIDDEN(WORK(KWMAT),ISYFKY,ISYMB,B, 1188 * ISYMD,D) 1189C 1190* call sum_pt3(work(kwmat),isymb,b,isymd,d, 1191* * iswmat,work(kx3am),4) 1192C 1193C------------------------------------------------------ 1194C Calculate the term <mu1|[H,W^BD(3)]|HF> 1195C added in XI1EFF 1196C------------------------------------------------------ 1197C 1198 CALL CC3_CY1(XI1EFF,ISYRES,WORK(KWMAT),ISWMAT, 1199 * WORK(KTMAT),WORK(KXIAJB), 1200 * ISINT1,WORK(KINDSQW),LENSQW, 1201 * WORK(KEND4),LWRK4, 1202 * ISYMB,B,ISYMD,D) 1203 1204C 1205C------------------------------------------------------ 1206C Calculate the term <mu2|[H,W^BD(3)]|HF> ( Fock matrix cont ) 1207C added in XI2EFF 1208C------------------------------------------------------ 1209C 1210 CALL CC3_CY2F(XI2EFF,ISYRES,WORK(KWMAT),ISWMAT, 1211 * WORK(KTMAT), FOCK0,ISINT1,WORK(KINDSQW), 1212 * LENSQW,WORK(KEND4),LWRK4, 1213 * ISYMB,B,ISYMD,D) 1214c 1215C------------------------------------------------------ 1216C Calculate the term <mu2|[H,W^BD(3)]|HF> ( Occupied cont ) 1217C added in XI2EFF 1218C------------------------------------------------------ 1219C 1220 CALL CC3_CY2O(XI2EFF,ISYRES,WORK(KWMAT),ISWMAT, 1221 * WORK(KTMAT),WORK(KTROC), 1222 * WORK(KTROC1),ISINT1,WORK(KEND4),LWRK4, 1223 * WORK(KINDSQW),LENSQW, 1224 * ISYMB,B,ISYMD,D) 1225C 1226C 1227C------------------------------------------------------ 1228C Calculate the term <mu2|[H,W^BD(3)]|HF> ( Vccupied cont ) 1229C added in XI2EFF 1230C------------------------------------------------------ 1231C 1232 CALL CC3_CY2V(XI2EFF,ISYRES,WORK(KRBJIA),WORK(KWMAT), 1233 * ISWMAT,WORK(KTMAT),WORK(KTRVI), 1234 * WORK(KTRVI1),ISINT1,WORK(KEND4),LWRK4, 1235 * WORK(KINDSQW),LENSQW, 1236 * ISYMB,B,ISYMD,D) 1237C 1238C-------------------------------------- 1239C Accumulate RMAT2 in XI2EFF 1240C-------------------------------------- 1241C 1242 IF (.NOT.(LCAUCHY .AND. (ICAU.GE.0))) THEN 1243 CALL CC3_RACC(XI2EFF,WORK(KRMAT2),ISYMB,B,ISYRES) 1244 END IF 1245c 1246 END DO 1247 END DO 1248C 1249C-------------------------------------- 1250C Accumulate RMAT1 in XI2EFF 1251C-------------------------------------- 1252C 1253 IF (.NOT.(LCAUCHY .AND. (ICAU.GE.0))) THEN 1254 CALL CC3_RACC(XI2EFF,WORK(KRMAT1),ISYMD,D,ISYRES) 1255 END IF 1256C 1257 END DO 1258 END DO 1259* write(lupri,*) 'TETAXKL in CC3_XISD : ' 1260* write(lupri,*) 'C3 in CC3_XISD, ncau = ',ncau 1261* call print_pt3(work(kx3am),1,4) 1262C 1263C 1264C------------------------------------------------------ 1265C Accumulate RBJIA from <mu2|[H,W^BD(3)]|HF> ( Vccupied cont ) 1266C in XI2EFF 1267C------------------------------------------------------ 1268C 1269 CALL CC3_RBJIA(XI2EFF,ISYRES,WORK(KRBJIA)) 1270c 1271c write(lupri,*)'XI1EFF in CC3_XISD' 1272c call output(XI1EFF,1,nvir(1),1,nrhf(1),nvir(1),nrhf(1),1,lupri) 1273c 1274c write(lupri,*)'XI2EFF in CC3_XISD' 1275c call outpak(XI2EFF,NT1AM(isyres),1,lupri) 1276C 1277 IF (IPRINT .GT. 55) THEN 1278 XXI2EFF = DDOT(NT2AM(ISYRES),XI2EFF,1,XI2EFF,1) 1279 WRITE(LUPRI,*) 'Norm of XI2EFF final before added ',XXI2EFF 1280 ENDIF 1281C 1282 DO I = 1,NT2AM(ISYRES) 1283 XI2EFF(I) = XI2EFF(I) + XI2(I) 1284 END DO 1285C 1286 IF (IPRINT .GT. 55) THEN 1287 XXI2EFF = DDOT(NT2AM(ISYRES),XI2EFF,1,XI2EFF,1) 1288 WRITE(LUPRI,*) 'Norm of XI2EFF final, xi2eff + xi2 ',XXI2EFF 1289 ENDIF 1290C 1291 DO I = 1,NT1AM(ISYRES) 1292 XI1EFF(I) = XI1EFF(I) + XI1(I) 1293 END DO 1294 1295* write(lupri,*)'xi1 at end of CC3_XISD' 1296* call PRINT_MATAI(XI1,ISYRES) 1297C 1298* write(lupri,*)'xi1eff at end of CC3_XISD' 1299* call PRINT_MATAI(XI1eff,ISYRES) 1300C 1301 IF (IPRINT .GT. 55) THEN 1302 XXI1EFF = DDOT(NT1AM(ISYRES),XI1EFF,1,XI1EFF,1) 1303 WRITE(LUPRI,*) 'Norm of XI1EFF final, xi1eff + xi1 ',XXI1EFF 1304 ENDIF 1305C 1306C--------------------------------------------------------------------- 1307C It might have happened that 'CC3_CAUINT_V*' files have not been 1308C deleted up there. Do it now! 1309C--------------------------------------------------------------------- 1310C 1311 IF (LCAUCHY .AND. ICAU.GE.0) THEN 1312 CALL DELETE_FILES('CC3_CAUINT_V*') 1313 END IF 1314C 1315C------------- 1316C End 1317C------------- 1318C 1319 CALL QEXIT('CC3_XISD ') 1320C 1321 RETURN 1322C 1323 END 1324C 1325C /* Deck cc3_t3bd */ 1326 SUBROUTINE CC3_T3BD(ISSMAT,SMAT,SMAT2, 1327 * UMAT,UMAT2,TMAT, 1328 * INDSQ,LENSQ) 1329C 1330C Written by P. Jorgensen and F. Pawlowski, Spring 2002. 1331C 1332C Sum up S and U intermediates for fixed b and d and 1333C write on TMAT 1334C 1335C tmat^{bd}(ckij) = smat^{bd}(ckij) + umat^{bd}(ckij) 1336C * + smat^{db}(ckji) + umat^{db}(ckji) 1337C 1338C ISSMAT is symmetry of S matrix (=symmetry of TMAT) 1339C 1340C 1341 IMPLICIT NONE 1342C 1343#include "ccsdsym.h" 1344C 1345 INTEGER ISSMAT, LENGTH, LENSQ, INDSQ(LENSQ,6) 1346C 1347 REAL*8 SMAT(*), SMAT2(*) 1348 REAL*8 UMAT(*), UMAT2(*), TMAT(*) 1349C 1350 CALL QENTER('CC3_T3BD') 1351C 1352 LENGTH = NCKIJ(ISSMAT) 1353 1354C 1355C--------------------------- 1356C Sum up intermediates 1357C--------------------------- 1358C 1359 DO I = 1, LENGTH 1360 TMAT(I) = 1361 * SMAT(I) 1362 * + UMAT(I) 1363 * + SMAT2(INDSQ(I,3)) 1364 * + UMAT2(INDSQ(I,3)) 1365 ENDDO 1366C 1367C--------------------- 1368C End. 1369C--------------------- 1370C 1371 CALL QEXIT('CC3_T3BD') 1372C 1373 RETURN 1374 END 1375 1376 1377c /* deck print_t3bd */ 1378c 1379 subroutine print_t3bd(tmat,isymim,b,c,iopt) 1380c 1381c Print t3^{BC}_{aikj} amplitudes which have been 1382c summed up (outside) to tmat array. (IOPT = 1) 1383c 1384c Print W^{BC}_{aikj} "amplitudes" which have been 1385c summed up (outside) to wmat array. (IOPT = 2) 1386c 1387c isymim is symmetry of (aikj) 1388c 1389c OBS! B and C are FIXED and coming from OUTSIDE 1390c 1391c 1392c P. Jorgensen and F. Pawlowski, Spring 2002. 1393c 1394 IMPLICIT NONE 1395c 1396#include "priunit.h" 1397#include "ccsdsym.h" 1398#include "ccorb.h" 1399C 1400 INTEGER ISYMIM, KOFF1, KOFF4, KOFF5, KOFF6, KOFF7 1401 INTEGER ISYMA, ISYMI, ISYMJ, ISYMK 1402 INTEGER KH, ISYIJK, ISYMJK, ISYMAI, ISYAIK 1403 INTEGER IOPT 1404C 1405 REAL*8 tmat(*) 1406C 1407 CALL QENTER('print_t3bd') 1408C 1409C---------------------------------------- 1410C Print the triples amplitudes. 1411C---------------------------------------- 1412C 1413 DO ISYMA = 1, NSYM 1414C 1415 KOFF1 = 0 1416 DO KH = 1, ISYMA-1 1417 KOFF1 = KOFF1 + NVIR(KH) 1418 ENDDO 1419C 1420 ISYIJK = MULD2H(ISYMIM,ISYMA) 1421C 1422 DO ISYMI = 1, NSYM 1423C 1424 KOFF4 = 0 1425 DO KH = 1, ISYMI-1 1426 KOFF4 = KOFF4 + NRHF(KH) 1427 ENDDO 1428C 1429 ISYMJK = MULD2H(ISYIJK,ISYMI) 1430 ISYMAI = MULD2H(ISYMA,ISYMI) 1431C 1432 DO ISYMJ = 1, NSYM 1433C 1434 KOFF5 = 0 1435 DO KH = 1, ISYMJ-1 1436 KOFF5 = KOFF5 + NRHF(KH) 1437 ENDDO 1438C 1439 ISYMK = MULD2H(ISYMJK,ISYMJ) 1440 ISYAIK = MULD2H(ISYMAI,ISYMK) 1441C 1442 KOFF6 = 0 1443 DO KH = 1, ISYMK-1 1444 KOFF6 = KOFF6 + NRHF(KH) 1445 ENDDO 1446C 1447 DO A = 1, NVIR(ISYMA) 1448 DO I = 1, NRHF(ISYMI) 1449 DO J = 1, NRHF(ISYMJ) 1450 DO K = 1, NRHF(ISYMK) 1451C 1452 KOFF7 = ISAIKJ(ISYAIK,ISYMJ) 1453 * + NCKI(ISYAIK)*(J - 1) 1454 * + ISAIK(ISYMAI,ISYMK) + NT1AM(ISYMAI)*(K-1) 1455 * + IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I-1) + A 1456 1457C 1458 IF (ABS(TMAT(KOFF7)) 1459 * .GT. 1.0D-12) THEN 1460 IF (IOPT .EQ. 1) THEN 1461 write(lupri,1) 'T3BD(',a+koff1,',',b,',', 1462 * c,',',i+koff4,',', 1463 * j+koff5,',',k+koff6,') = ', 1464 * TMAT(KOFF7) 1465 ELSE IF (IOPT .EQ. 2) THEN 1466 write(lupri,1) 'WBD(',a+koff1,',',b,',', 1467 * c,',',i+koff4,',', 1468 * j+koff5,',',k+koff6,') = ', 1469 * TMAT(KOFF7) 1470 ELSE 1471 call quit('Illegal value for IOPT in print_t3bd ') 1472 ENDIF 1473C 1474 ENDIF 1475C 1476 ENDDO ! K 1477 ENDDO ! J 1478 ENDDO ! I 1479 ENDDO ! A 1480 ENDDO ! ISYMJ 1481 ENDDO ! ISYMI 1482C 1483 ENDDO ! ISYMA 1484C 1485 CALL QEXIT('print_t3bd') 1486C 1487 1 FORMAT(1X,A6,I3,A1,I3,A1,I3,A1,I3,A1,I3,A1,I3,A4,E20.10) 1488 RETURN 1489 END 1490C 1491C /* Deck wbd_v */ 1492 SUBROUTINE WBD_V(TMAT,ISTMAT,FOCKY,ISYFKY, 1493 * WMAT,ISWMAT,WRK,LWRK) 1494C 1495C WBD(a,i,k,j) = WBD(a,i,k,j) + sum (f) focky(a,f)*tmatBD(f,i,k,j) 1496C 1497C 1498C Written by P. Jorgensen and F. Pawlowski, Spring 2002. 1499C 1500 1501C 1502 IMPLICIT NONE 1503C 1504 INTEGER LWRK, KFCAF, KEND0, LWRK0, KOFF1, KOFF2 1505 INTEGER NF, KOFFY, KOFFT, KOFFW 1506 INTEGER ISTMAT, ISYFKY, ISWMAT, ISFIKJ, ISYFIK 1507 INTEGER ISYMA, ISYAI, ISYAIK, NA 1508 INTEGER ISYMF, ISYMJ, ISYMK, ISYMI, ISYFI 1509C 1510 REAL*8 TMAT(*), FOCKY(*), WMAT(*), WRK(*) 1511 REAL*8 HALF, ONE 1512C 1513#include "priunit.h" 1514#include "dummy.h" 1515#include "iratdef.h" 1516#include "ccsdsym.h" 1517#include "inftap.h" 1518#include "ccinftap.h" 1519#include "ccorb.h" 1520#include "ccsdinp.h" 1521C 1522 PARAMETER (HALF = 0.5D0, ONE = 1.0D0) 1523C 1524 CALL QENTER('WBD_V') 1525C 1526C 1527C RESORT VIR-VIR FOCKY ELEMENTS (A,F) 1528C 1529C 1530 KFCAF = 1 1531 KEND0 = KFCAF + NMATAB(ISYFKY) 1532 LWRK0 = LWRK - KEND0 1533C 1534 IF (LWRK0 .LT. 0) THEN 1535 WRITE(LUPRI,*) 'Memory available : ',LWRK0 1536 WRITE(LUPRI,*) 'Memory needed : ',KEND0 1537 CALL QUIT('Insufficient space in WBD_V') 1538 END IF 1539C 1540 DO ISYMF = 1,NSYM 1541 ISYMA = MULD2H(ISYMF,ISYFKY) 1542 DO F = 1,NVIR(ISYMF) 1543 KOFF1 = IFCVIR(ISYMA,ISYMF) + NORB(ISYMA)*(F - 1) 1544 * + NRHF(ISYMA) + 1 1545 KOFF2 = KFCAF + IMATAB(ISYMA,ISYMF) + NVIR(ISYMA)*(F - 1) 1546 CALL DCOPY(NVIR(ISYMA),FOCKY(KOFF1),1,WRK(KOFF2),1) 1547 END DO 1548 END DO 1549C 1550C CARRY OUT MATRIX MULTIPLICATION 1551C WBD(a,i,k,j) = WBD(a,i,k,j) + sum (f) focky(a,f)*tmatBD(f,i,k,j) 1552C 1553 ISFIKJ = ISTMAT 1554 DO ISYMJ = 1,NSYM 1555 ISYFIK =MULD2H(ISYMJ,ISFIKJ) 1556 DO J = 1,NRHF(ISYMJ) 1557 DO ISYMK = 1,NSYM 1558 ISYFI = MULD2H(ISYMK,ISYFIK) 1559 DO K = 1,NRHF(ISYMK) 1560 DO ISYMI = 1,NSYM 1561 ISYMF = MULD2H(ISYFI,ISYMI) 1562 ISYMA = MULD2H(ISYFKY,ISYMF) 1563 ISYAIK = MULD2H(ISWMAT,ISYMJ) 1564 ISYAI = MULD2H(ISYAIK,ISYMK) 1565 NA = MAX(1,NVIR(ISYMA)) 1566 NF = MAX(1,NVIR(ISYMF)) 1567 KOFFY = KFCAF + IMATAB(ISYMA,ISYMF) 1568 KOFFT = ISAIKJ(ISYFIK,ISYMJ) 1569 * + NCKI(ISYFIK)*(J-1) 1570 * + ISAIK(ISYFI,ISYMK) 1571 * + NT1AM(ISYFI)*(K-1) 1572 * + IT1AM(ISYMF,ISYMI) + 1 1573 KOFFW = ISAIKJ(ISYAIK,ISYMJ) 1574 * + NCKI(ISYAIK)*(J-1) 1575 * + ISAIK(ISYAI,ISYMK) 1576 * + NT1AM(ISYAI)*(K-1) 1577 * + IT1AM(ISYMA,ISYMI) + 1 1578C 1579C SYMMETRY BETWEEN BJ AND CK INTRODUCE A FACTOR TWO 1580C DENOTE t3 IS CALCULATED WITH NEGATIVE SIGN 1581C 1582 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI), 1583 * NVIR(ISYMF),-ONE,WRK(KOFFY),NA, 1584 * TMAT(KOFFT),NF,ONE,WMAT(KOFFW),NA) 1585 END DO 1586 END DO 1587 END DO 1588 END DO 1589 END DO 1590C 1591 CALL QEXIT('WBD_V') 1592C 1593 RETURN 1594 END 1595C 1596C /* Deck wbd_o */ 1597 SUBROUTINE WBD_O(TMAT,ISTMAT,FOCKY,ISYFKY, 1598 * WMAT,ISWMAT,WRK,LWRK) 1599C 1600C WBD(a,i,k,j) = WBD(a,i,k,j) - sum (l) tmatBD(a,l,k,j)*focky(l,i) 1601C 1602C 1603C Written by P. Jorgensen and F. Pawlowski, Spring 2002. 1604C 1605 1606 IMPLICIT NONE 1607C 1608 INTEGER LWRK, KFCLI, KEND0, LWRK0, KOFF1, KOFF2 1609 INTEGER NL, KOFFY, KOFFT, KOFFW 1610 INTEGER ISTMAT, ISYFKY, ISWMAT, ISALKJ 1611 INTEGER ISYMA, ISYAI, ISYAIK, ISYALK, ISYAL, NA 1612 INTEGER ISYMJ, ISYMK, ISYMI, ISYML, ISYFI 1613C 1614 REAL*8 TMAT(*), FOCKY(*), WMAT(*), WRK(*) 1615 REAL*8 MHALF, ONE 1616C 1617#include "priunit.h" 1618#include "dummy.h" 1619#include "iratdef.h" 1620#include "ccsdsym.h" 1621#include "inftap.h" 1622#include "ccinftap.h" 1623#include "ccorb.h" 1624#include "ccsdinp.h" 1625C 1626 PARAMETER (MHALF = -0.5D0, ONE = 1.0D0) 1627C 1628 CALL QENTER('WBD_O') 1629C 1630 1631C 1632C RESORT OCC-OCC FOCKY ELEMENTS (L,I) 1633C 1634C 1635 KFCLI = 1 1636 KEND0 = KFCLI + NMATIJ(ISYFKY) 1637 LWRK0 = LWRK - KEND0 1638C 1639 IF (LWRK0 .LT. 0) THEN 1640 WRITE(LUPRI,*) 'Memory available : ',LWRK0 1641 WRITE(LUPRI,*) 'Memory needed : ',KEND0 1642 CALL QUIT('Insufficient space in WBD_O') 1643 END IF 1644C 1645 DO ISYMI = 1,NSYM 1646 ISYML = MULD2H(ISYMI,ISYFKY) 1647 DO I = 1,NRHF(ISYMI) 1648 KOFF1 = IFCRHF(ISYML,ISYMI) + NORB(ISYML)*(I - 1) + 1 1649 KOFF2 = KFCLI + IMATIJ(ISYML,ISYMI) + NRHF(ISYML)*(I - 1) 1650 CALL DCOPY(NRHF(ISYML),FOCKY(KOFF1),1,WRK(KOFF2),1) 1651 END DO 1652 END DO 1653C 1654C CARRY OUT MATRIX MULTIPLICATION 1655C WBD(a,i,k,j) = WBD(a,i,k,j) - sum (l) tmatBD(a,l,k,j)*focky(l,i) 1656C 1657 ISALKJ = ISTMAT 1658 DO ISYMJ = 1,NSYM 1659 ISYALK =MULD2H(ISYMJ,ISALKJ) 1660 DO J = 1,NRHF(ISYMJ) 1661 DO ISYMK = 1,NSYM 1662 ISYAL = MULD2H(ISYMK,ISYALK) 1663 DO K = 1,NRHF(ISYMK) 1664 DO ISYML = 1,NSYM 1665 ISYMA = MULD2H(ISYAL,ISYML) 1666 ISYMI = MULD2H(ISYFKY,ISYML) 1667 ISYAIK = MULD2H(ISWMAT,ISYMJ) 1668 ISYAI = MULD2H(ISYAIK,ISYMK) 1669 NA = MAX(1,NVIR(ISYMA)) 1670 NL = MAX(1,NRHF(ISYML)) 1671 KOFFY = KFCLI + IMATIJ(ISYML,ISYMI) 1672 KOFFT = ISAIKJ(ISYALK,ISYMJ) 1673 * + NCKI(ISYALK)*(J-1) 1674 * + ISAIK(ISYAL,ISYMK) 1675 * + NT1AM(ISYAL)*(K-1) 1676 * + IT1AM(ISYMA,ISYML) + 1 1677 KOFFW = ISAIKJ(ISYAIK,ISYMJ) 1678 * + NCKI(ISYAIK)*(J-1) 1679 * + ISAIK(ISYAI,ISYMK) 1680 * + NT1AM(ISYAI)*(K-1) 1681 * + IT1AM(ISYMA,ISYMI) + 1 1682C 1683C SYMMETRY BETWEEN BJ AND CK INTRODUCE A FACTOR TWO 1684C DENOTE t3 IS CALCULATED WITH NEGATIVE SIGN 1685C 1686 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI), 1687 * NRHF(ISYML),ONE,TMAT(KOFFT),NA, 1688 * WRK(KOFFY),NL,ONE,WMAT(KOFFW),NA) 1689 END DO 1690 END DO 1691 END DO 1692 END DO 1693 END DO 1694C 1695 CALL QEXIT('WBD_O') 1696C 1697 RETURN 1698 END 1699C /* Deck wbd_t2 */ 1700 SUBROUTINE WBD_T2(T2XNET2Z,B,ISYMB,D,ISYMD, 1701 * T2TPX,ISYMT2X,T2TPZ,ISYMT2Z, 1702 * FOCKY,ISYFKY, 1703 * INDSQ,LENSQ,WMAT,ISWMAT,WRK,LWRK) 1704* 1705********************************************************************** 1706* 1707* Calculate the term <mu3|[[Y,T2X],T2Z]|HF> as the contribution to 1708* WBC intermmediate. 1709* 1710* If T2XNET2Z = .true. then calculate: 1711* 1712* what is returned is: 1713* 1714* WBC(a,i,k,j) = WBC(a,i,k,j) - 1715* P(bj,ck) sum(l,d) focky(l,d)*( t2X(ai,cl)*t2Z(bj,dk)+t2Z(ai,cl)*t2X(bj,dk)) 1716* 1717* ELSE 1718* 1719* <mu3|[[Y,T2],T2]|HF> = P(aibjck) 2 sum(l,d) focky(l,d)*(t2X(ai,cl)*t2Z(bj,dk)) 1720* 1721* what is returned by the routine is WBC contribution corresponding to 1722* (ai,bj,ck) + (ai,ck,bj) permutation. 1723* 1724* Totally is returned: (note that 2 factor is absorbed usually by 1/2 1725* in the formula; for example in T3X !!!) 1726* 1727* WBC(a,i,k,j) = WBC(a,i,k,j) - 1728* P(bj,ck) sum(l,d) focky(l,d)*( t2X(ai,cl)*t2Z(bj,dk) ) 1729* 1730* 1731* F. Pawlowski, P. Jorgensen, Aarhus 14-May-2003. 1732* 1733* This is just a little driver which calls WBD_T2_1. 1734* 1735********************************************************************** 1736* 1737 IMPLICIT NONE 1738C 1739#include "priunit.h" 1740#include "ccsdsym.h" 1741C 1742 LOGICAL T2XNET2Z 1743C 1744 INTEGER LENSQ, INDSQ(LENSQ,6), LWRK 1745 INTEGER ISYMB, ISYMD, ISYMT2X, ISYFKY, ISWMAT, ISYMT2Z 1746C 1747 REAL*8 T2TPX(*), FOCKY(*), WMAT(*), WRK(*), T2TPZ(*) 1748 REAL*8 FAC,HALF,ONE 1749C 1750 PARAMETER (HALF = 0.5D0, ONE = 1.0D0 ) 1751C 1752 CALL QENTER('WBDT2') 1753C 1754 IF (T2XNET2Z) THEN 1755 FAC = ONE 1756 ELSE 1757 FAC = ONE 1758 END IF 1759C 1760 CALL WBD_T2_1(FAC,B,ISYMB,D,ISYMD, 1761 * T2TPX,ISYMT2X,T2TPZ,ISYMT2Z, 1762 * FOCKY,ISYFKY, 1763 * INDSQ,LENSQ,WMAT,ISWMAT,WRK,LWRK) 1764C 1765 IF (T2XNET2Z) THEN 1766C 1767 CALL WBD_T2_1(FAC,B,ISYMB,D,ISYMD, 1768 * T2TPZ,ISYMT2Z,T2TPX,ISYMT2X, 1769 * FOCKY,ISYFKY, 1770 * INDSQ,LENSQ,WMAT,ISWMAT,WRK,LWRK) 1771C 1772 END IF 1773C 1774 CALL QEXIT('WBDT2') 1775C 1776 RETURN 1777 END 1778C /* Deck wbd_t2_1 */ 1779 SUBROUTINE WBD_T2_1(FAC,B,ISYMB,D,ISYMD,T2TPX,ISYMT2X,T2TPZ, 1780 * ISYMT2Z, 1781 * FOCKY,ISYFKY, 1782 * INDSQ,LENSQ,WMAT,ISWMAT,WRK,LWRK) 1783C 1784C WBD(a,i,k,j) = WBD(a,i,k,j) - 1785C sum (f,l) focky(l,f)*( t2X(ai,dl)*t2Z(fk,bj) + t2X(ai,bl)*t2Z(fj,dk) ) 1786C 1787C 1788C 1789C Written by P. Jorgensen and F. Pawlowski, Spring 2002. 1790C 1791C Modified to handle two different symmetries of T2, 14-May-2003. 1792 1793 IMPLICIT NONE 1794C 1795#include "priunit.h" 1796#include "dummy.h" 1797#include "iratdef.h" 1798#include "ccsdsym.h" 1799#include "inftap.h" 1800#include "ccinftap.h" 1801#include "ccorb.h" 1802#include "ccsdinp.h" 1803C 1804 INTEGER LENSQ 1805 INTEGER INDSQ(LENSQ,6) 1806 INTEGER LWRK, KFCLF, KEND0, LWRK0, KOFF1, KOFF2, KTB, KEND1, LWRK1 1807 INTEGER NL, NF, KOFFY, KOFFT2, KOFFT, KOFFW, KTD, KW 1808 INTEGER ISYMB, ISYMD, ISYMT2X, ISYFKY, ISWMAT 1809 INTEGER ISYAIL, ISYAI, ISYAIK, NA, NAI, LENGTH 1810 INTEGER ISYMF, ISYML, ISYFKJ, ISYTB, ISYMJ, ISYFK, ISYMK, ISYLK 1811 INTEGER ISYFJK, ISYTD, ISYLJ, ISYFJ, ISYAIJ 1812 INTEGER ISYMT2Z 1813C 1814 REAL*8 T2TPX(*), FOCKY(*), WMAT(*), WRK(*) 1815 REAL*8 HALF, ONE, ZERO, FAC 1816 REAL*8 T2TPZ(*) 1817C 1818 PARAMETER (HALF = 0.5D0, ONE = 1.0D0, ZERO = 0.0D0) 1819C 1820 CALL QENTER('WT21') 1821C 1822C 1823C RESORT VIR-OCC FOCKY ELEMENTS (l,f) 1824C 1825C 1826 KW = 1 1827 KFCLF = KW + NCKIJ(ISWMAT) 1828 KEND0 = KFCLF + NT1AM(ISYFKY) 1829 LWRK0 = LWRK - KEND0 1830 CALL DZERO(WRK(KW),NCKIJ(ISWMAT)) 1831C 1832 IF (LWRK0 .LT. 0) THEN 1833 WRITE(LUPRI,*) 'Memory available : ',LWRK0 1834 WRITE(LUPRI,*) 'Memory needed : ',KEND0 1835 CALL QUIT('Insufficient space in WBD_T2_1 (1)') 1836 END IF 1837C 1838 DO ISYMF = 1,NSYM 1839 ISYML = MULD2H(ISYMF,ISYFKY) 1840 DO L = 1,NRHF(ISYML) 1841 DO F = 1,NVIR(ISYMF) 1842 KOFF1 = IFCVIR(ISYML,ISYMF) + NORB(ISYML)*(F - 1) + L 1843 KOFF2 = KFCLF + IT1AMT(ISYML,ISYMF) 1844 * + NRHF(ISYML)*(F - 1) + L -1 1845C 1846 WRK(KOFF2) = FOCKY(KOFF1) 1847C 1848 END DO 1849 END DO 1850 END DO 1851C 1852C calculate first t2 contribution to W matrix 1853C 1854C construct tZB(l,k,j) = sum (f) focky(l,f)*t2tpZ(f,k,j,B) 1855C 1856 ISYFKJ = MULD2H(ISYMT2Z,ISYMB) 1857 ISYTB = MULD2H(ISYFKY,ISYFKJ) 1858 KTB = KEND0 1859 KEND1 = KTB + NMAIJK(ISYTB) 1860 LWRK1 = LWRK - KEND1 1861C 1862 CALL DZERO(WRK(KTB),NMAIJK(ISYTB)) 1863C 1864 IF (LWRK1 .LT. 0) THEN 1865 WRITE(LUPRI,*) 'Memory available : ',LWRK1 1866 WRITE(LUPRI,*) 'Memory needed : ',KEND1 1867 CALL QUIT('Insufficient space in WBD_T2_1 (2)') 1868 END IF 1869C 1870 DO ISYMJ = 1,NSYM 1871 ISYFK = MULD2H(ISYFKJ,ISYMJ) 1872 DO J = 1,NRHF(ISYMJ) 1873 DO ISYMK = 1,NSYM 1874 ISYMF = MULD2H(ISYFK,ISYMK) 1875 ISYML = MULD2H(ISYFKY,ISYMF) 1876 ISYLK = MULD2H(ISYML,ISYMK) 1877 NL = MAX(1,NRHF(ISYML)) 1878 NF = MAX(1,NVIR(ISYMF)) 1879 KOFFY = KFCLF + IT1AMT(ISYML,ISYMF) 1880 KOFFT2 = IT2SP(ISYFKJ,ISYMB) + NCKI(ISYFKJ)*(B-1) 1881 * + ISAIK(ISYFK,ISYMJ) + NT1AM(ISYFK)*(J-1) 1882 * + IT1AM(ISYMF,ISYMK) + 1 1883 KOFFT = KTB + IMAIJK(ISYLK,ISYMJ) + NMATIJ(ISYLK)*(J-1) 1884 * + IMATIJ(ISYML,ISYMK) 1885C 1886 CALL DGEMM('N','N',NRHF(ISYML),NRHF(ISYMK), 1887 * NVIR(ISYMF),ONE,WRK(KOFFY),NL, 1888 * T2TPZ(KOFFT2),NF,ONE,WRK(KOFFT),NL) 1889C 1890 END DO 1891 END DO 1892 END DO 1893C 1894C WBD(a,i,k,j) = WBD(a,i,k,j) - 1895C sum (f,l) focky(l,f)* t2X(ai,Dl)*t2Z(fk,Bj) 1896C = WBD(a,i,k,j) - 1897C sum(l) t2tpX(a,i,l,D) * tZB(l,k,j) 1898C 1899 ISYAIL = MULD2H(ISYMT2X,ISYMD) 1900 DO ISYMJ = 1,NSYM 1901 ISYLK = MULD2H(ISYTB,ISYMJ) 1902 DO J = 1,NRHF(ISYMJ) 1903 DO ISYMK = 1,NSYM 1904 ISYML = MULD2H(ISYLK,ISYMK) 1905 ISYAI = MULD2H(ISYAIL,ISYML) 1906 ISYAIK = MULD2H(ISYAI,ISYMK) 1907 NAI = MAX(1,NT1AM(ISYAI)) 1908 NL = MAX(1,NRHF(ISYML)) 1909 KOFFT2 = IT2SP(ISYAIL,ISYMD) + NCKI(ISYAIL)*(D-1) 1910 * + ISAIK(ISYAI,ISYML) + 1 1911 KOFFT = KTB + IMAIJK(ISYLK,ISYMJ) + NMATIJ(ISYLK)*(J-1) 1912 * + IMATIJ(ISYML,ISYMK) 1913 KOFFW = ISAIKJ(ISYAIK,ISYMJ) + NCKI(ISYAIK)*(J-1) 1914 * + ISAIK(ISYAI,ISYMK) + 1 1915 CALL DGEMM('N','N',NT1AM(ISYAI),NRHF(ISYMK), 1916 * NRHF(ISYML),-FAC,T2TPX(KOFFT2),NAI, 1917 * WRK(KOFFT),NL,ONE,WMAT(KOFFW),NAI) 1918 1919C 1920 END DO 1921 END DO 1922 END DO 1923C 1924C calculate second t2 contribution to W matrix 1925C 1926C construct tZD(l,j,k) = sum (f) focky(l,f)*t2tpZ(f,j,k,D) 1927C 1928 ISYFJK = MULD2H(ISYMT2Z,ISYMD) 1929 ISYTD = MULD2H(ISYFKY,ISYFJK) 1930 KTD = KEND0 1931 KEND1 = KTD + NMAIJK(ISYTD) 1932 LWRK1 = LWRK - KEND1 1933C 1934 CALL DZERO(WRK(KTD),NMAIJK(ISYTD)) 1935C 1936 IF (LWRK1 .LT. 0) THEN 1937 WRITE(LUPRI,*) 'Memory available : ',LWRK1 1938 WRITE(LUPRI,*) 'Memory needed : ',KEND1 1939 CALL QUIT('Insufficient space in WBD_T2_1 (3)') 1940 END IF 1941C 1942 1943 DO ISYMK = 1,NSYM 1944 ISYFJ = MULD2H(ISYFJK,ISYMK) 1945 DO K = 1,NRHF(ISYMK) 1946 DO ISYMJ = 1,NSYM 1947 ISYMF = MULD2H(ISYFJ,ISYMJ) 1948 ISYML = MULD2H(ISYFKY,ISYMF) 1949 ISYLJ = MULD2H(ISYML,ISYMJ) 1950 NL = MAX(1,NRHF(ISYML)) 1951 NF = MAX(1,NVIR(ISYMF)) 1952 KOFFY = KFCLF + IT1AMT(ISYML,ISYMF) 1953 KOFFT2 = IT2SP(ISYFJK,ISYMD) + NCKI(ISYFJK)*(D-1) 1954 * + ISAIK(ISYFJ,ISYMK) + NT1AM(ISYFJ)*(K-1) 1955 * + IT1AM(ISYMF,ISYMJ) + 1 1956 KOFFT = KTD + IMAIJK(ISYLJ,ISYMK) + NMATIJ(ISYLJ)*(K-1) 1957 * + IMATIJ(ISYML,ISYMJ) 1958 CALL DGEMM('N','N',NRHF(ISYML),NRHF(ISYMJ), 1959 * NVIR(ISYMF),ONE,WRK(KOFFY),NL, 1960 * T2TPZ(KOFFT2),NF,ONE,WRK(KOFFT),NL) 1961C 1962 END DO 1963 END DO 1964 END DO 1965C 1966C WBD(a,i,k,j) = WBD(a,i,k,j) - 1967C sum (f,l) focky(l,f)*t2X(ai,Bl)*t2Z(fj,Dk) ) 1968C = WBD(a,i,k,j) - 1969C sum(l) t2tpX(a,i,l,B) * tZD(l,j,k) 1970C 1971 ISYAIL = MULD2H(ISYMT2X,ISYMB) 1972 DO ISYMK = 1,NSYM 1973 ISYLJ = MULD2H(ISYTD,ISYMK) 1974 DO K = 1,NRHF(ISYMK) 1975 DO ISYMJ = 1,NSYM 1976 ISYML = MULD2H(ISYLJ,ISYMJ) 1977 ISYAI = MULD2H(ISYAIL,ISYML) 1978 ISYAIJ = MULD2H(ISYAI,ISYMJ) 1979 NAI = MAX(1,NT1AM(ISYAI)) 1980 NL = MAX(1,NRHF(ISYML)) 1981 KOFFT2 = IT2SP(ISYAIL,ISYMB) + NCKI(ISYAIL)*(B-1) 1982 * + ISAIK(ISYAI,ISYML) + 1 1983 KOFFT = KTD + IMAIJK(ISYLJ,ISYMK) + NMATIJ(ISYLJ)*(K-1) 1984 * + IMATIJ(ISYML,ISYMJ) 1985 KOFFW = KW + ISAIKJ(ISYAIJ,ISYMK) + NCKI(ISYAIJ)*(K-1) 1986 * + ISAIK(ISYAI,ISYMJ) 1987 CALL DGEMM('N','N',NT1AM(ISYAI),NRHF(ISYMJ), 1988 * NRHF(ISYML),-FAC,T2TPX(KOFFT2),NAI, 1989 * WRK(KOFFT),NL,ONE,WRK(KOFFW),NAI) 1990 1991C 1992 END DO 1993 END DO 1994 END DO 1995C 1996C change order aijk to aikj 1997C 1998 DO I = 1,NCKIJ(ISWMAT) 1999 WMAT(I) = WMAT(I) + WRK(INDSQ(I,3)) 2000 END DO 2001C 2002 CALL QEXIT('WT21') 2003C 2004 RETURN 2005 END 2006 2007 2008C /* Deck cc3_cY1 */ 2009 SUBROUTINE CC3_CY1(OMEGA1,ISOMG1,WMAT,ISWMAT,TMAT, 2010 * XIAJB,ISYINT,INDSQ,LENSQ,WORK,LWORK, 2011 * ISYMIB,IB,ISYMID,ID) 2012C 2013C P. Jorgensen and F. Pawlowski. Spring 2002 2014C 2015 IMPLICIT NONE 2016C 2017C Calculate Omega1 2018C 2019C 2020C General symmetry: ISWMAT is symmetry of WMAT and TMAT intermdiates. 2021C (exclusive isymb,isymd) 2022C ISYINT = XIAJB 2023C ISOMG1 = ISWMAT*ISYINT 2024C 2025C Omega1(ai) 2026C 2027C = sum_{bjck} ( W^{abc}_{ijk} - W^{abc}_{kji} ) * L_{jbkc} 2028C (1): 2029C = sum_{jk} ( W^{BC}_{aikj} - W^{BC}_{akij} ) * L_{CkBj} 2030C (2): 2031C + sum_{bjk} ( W^{CA}_{bjik} - W^{AC}_{bjki} ) * L_{bjCk} 2032C (3): 2033C + sum_{cjk} ( W^{AB}_{ckji} - W^{AB}_{cijk} ) * L_{ckBj} 2034C 2035C 2036C 2037 INTEGER ISOMG1, ISWMAT, ISYINT, LENSQ, LWORK 2038 INTEGER ISYMIB, IB, ISYMID, ID, INDEX, ISYBJK 2039 INTEGER ISYBC, ISYMB, ISYMC, ISYKJ, ISYAI, ISYMJ, ISYBJ 2040 INTEGER ISYMK, NKJ, NBJ, NCK, NCKBJ, NTOTAI, KOFF1, NBJK 2041 INTEGER ISYMA, ISYMJK, ISYMI, ISYCK, NBJCK, NTOTA, NTOBJK 2042 INTEGER KOFF2, ISYCKJ, NCKJ, NTOCKJ, KSTART, KEND, LEND 2043 INTEGER INDSQ(LENSQ,6) 2044 REAL*8 OMEGA1(*), WMAT(*), TMAT(*), XIAJB(*), WORK(*) 2045 REAL*8 ONE, ZERO, TWO 2046 2047C 2048#include "priunit.h" 2049#include "ccsdsym.h" 2050#include "ccorb.h" 2051#include "ccsdinp.h" 2052C 2053 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 2054C 2055 INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 2056C 2057 CALL QENTER('CC3_CY1') 2058C 2059C---------------------------------- 2060C Calculate (1): 2061C sum_{jk} ( W^{BC}_{aikj} - W^{BC}_{akij} ) * L_{CkBj} 2062C---------------------------------- 2063C 2064C 2065C maximum work space 2066C 2067 KSTART = 1 2068 KEND = KSTART + NRHFT*NRHFT*NVIRT 2069 LEND = LWORK - KEND 2070C 2071 IF (LEND .LT. 0) THEN 2072 WRITE(LUPRI,*) 'Memory available : ',LWORK 2073 WRITE(LUPRI,*) 'Memory needed : ',KEND 2074 CALL QUIT('Insufficient space in CC3_CY1') 2075 END IF 2076 2077 B = IB 2078 C = ID 2079 ISYMB = ISYMIB 2080 ISYMC = ISYMID 2081C 2082C T(aikj) = W^{BC}_{aikj} - W^{BC}_{akij} ) 2083C 2084 DO I = 1,LENSQ 2085 TMAT(I) = WMAT(I) 2086 * - WMAT(INDSQ(I,1)) 2087 ENDDO 2088C 2089 IF (NSYM .GT. 1) THEN 2090 CALL CC_GATHER(LENSQ,WORK(KEND),TMAT,INDSQ(1,6)) 2091 CALL DCOPY(LENSQ,WORK(KEND),1,TMAT,1) 2092 ENDIF 2093C 2094 ISYBC = MULD2H(ISYMB,ISYMC) 2095 ISYKJ = MULD2H(ISYBC,ISYINT) 2096 ISYAI = ISOMG1 2097C 2098C Construct LCB(k,j) = L(Ck,Bj) 2099C --------------------------- 2100C 2101 DO 400 ISYMJ = 1,NSYM 2102C 2103 ISYMK = MULD2H(ISYMJ,ISYKJ) 2104 ISYCK = MULD2H(ISYMC,ISYMK) 2105 ISYBJ = MULD2H(ISYMB,ISYMJ) 2106C 2107 DO 410 J = 1,NRHF(ISYMJ) 2108C 2109 NBJ = IT1AM(ISYMB,ISYMJ) + NVIR(ISYMB)*(J - 1) + B 2110C 2111 DO 420 K = 1,NRHF(ISYMK) 2112C 2113 NKJ = IMATIJ(ISYMK,ISYMJ)+ NRHF(ISYMK)*(J - 1) + K 2114 NCK = IT1AM(ISYMC,ISYMK) + NVIR(ISYMC)*(K - 1) + C 2115C 2116 NCKBJ = IT2AM(ISYCK,ISYBJ) + INDEX(NCK,NBJ) 2117C 2118 WORK(NKJ) = XIAJB(NCKBJ) 2119C 2120 420 CONTINUE 2121 410 CONTINUE 2122 400 CONTINUE 2123C 2124 NTOTAI = MAX(NT1AM(ISYAI),1) 2125C 2126 KOFF1 = ISAIKL(ISYAI,ISYKJ) + 1 2127C 2128 CALL DGEMV('N',NT1AM(ISYAI),NMATIJ(ISYKJ),-ONE,TMAT(KOFF1), 2129 * NTOTAI,WORK,1,ONE,OMEGA1,1) 2130C 2131C 2132C Calculate (2): 2133C 2134C Omega1(ai) = Omega1(ai) 2135C + sum_{jk} ( W^{AC}_{bjki} - W^{AC}_{bjik} ) * L_{bjCk} 2136C 2137C + sum_{bjk} ( W^{CA}_{bjik} - W^{CA}_{bjki} ) * L_{bjCk} 2138C 2139 C = IB 2140 A = ID 2141c write(lupri,*)' start term 2 a,id',a,id 2142C 2143 ISYMC = ISYMIB 2144 ISYMA = ISYMID 2145C 2146C 2147C T(bjki) = W^{CA}_{bjik} - W^{AC}_{bjki} ) 2148C 2149 DO I = 1,LENSQ 2150C 2151 TMAT(I) = - WMAT(I) 2152 * + WMAT(INDSQ(I,3)) 2153C 2154 ENDDO 2155C 2156C Construct LC(bj,k) = L(bj,Ck) 2157C --------------------------- 2158C 2159 ISYBJK = MULD2H(ISYINT,ISYMC) 2160 ISYMI = MULD2H(ISYBJK,ISWMAT) 2161 DO 100 ISYMK = 1,NSYM 2162C 2163 ISYCK = MULD2H(ISYMC,ISYMK) 2164 ISYBJ = MULD2H(ISYMK,ISYBJK) 2165C 2166 DO 110 K = 1,NRHF(ISYMK) 2167C 2168 NCK = IT1AM(ISYMC,ISYMK) + NVIR(ISYMC)*(K - 1) + C 2169C 2170 DO 120 NBJ = 1,NT1AM(ISYBJ) 2171C 2172 NBJCK = IT2AM(ISYBJ,ISYCK) + INDEX(NBJ,NCK) 2173 NBJK = ICKI(ISYBJ,ISYMK) 2174 * + NT1AM(ISYBJ)*(K - 1) + NBJ 2175C 2176 WORK(NBJK) = XIAJB(NBJCK) 2177C 2178 120 CONTINUE 2179 110 CONTINUE 2180 100 CONTINUE 2181C 2182 NTOTA = MAX(NVIR(ISYMA),1) 2183 NTOBJK = MAX(NCKI(ISYBJK),1) 2184C 2185 KOFF1 = ISAIKJ(ISYBJK,ISYMI) + 1 2186 KOFF2 = IT1AM(ISYMA,ISYMI) + A 2187C 2188 CALL DGEMV('T',NCKI(ISYBJK),NRHF(ISYMI),-ONE,TMAT(KOFF1), 2189 * NTOBJK,WORK,1,ONE,OMEGA1(KOFF2),NTOTA) 2190C 2191C Calculate (3): 2192C 2193C Omega1(ai) = Omega1(ai) 2194C + sum_{jk} ( W^{AB}_{ckji} - W^{AB}_{cijk} ) * L_{ckBj} 2195C 2196C 2197 A = IB 2198 B = ID 2199C 2200 ISYMA = ISYMIB 2201 ISYMB = ISYMID 2202C 2203C 2204C T(ckji) = W^{AB}_{ckji} - W^{AB}_{cijk} 2205C 2206 DO I = 1,LENSQ 2207C 2208 TMAT(I) = WMAT(I) 2209 * - WMAT(INDSQ(I,5)) 2210C 2211 ENDDO 2212C 2213C Construct LB(ck,j) = L(ck,Bj) 2214C ---------------------------- 2215C 2216 ISYCKJ = MULD2H(ISYINT,ISYMB) 2217 ISYMI = MULD2H(ISYCKJ,ISWMAT) 2218 DO 200 ISYMJ = 1,NSYM 2219C 2220 ISYBJ = MULD2H(ISYMB,ISYMJ) 2221 ISYCK = MULD2H(ISYMJ,ISYCKJ) 2222C 2223 DO 210 J = 1,NRHF(ISYMJ) 2224C 2225 NBJ = IT1AM(ISYMB,ISYMJ) + NVIR(ISYMB)*(J - 1) + B 2226C 2227 DO 220 NCK = 1,NT1AM(ISYCK) 2228C 2229 NCKBJ = IT2AM(ISYCK,ISYBJ) + INDEX(NCK,NBJ) 2230 NCKJ = ICKI(ISYCK,ISYMJ) 2231 * + NT1AM(ISYCK)*(J - 1) + NCK 2232C 2233 WORK(NCKJ) = XIAJB(NCKBJ) 2234C 2235 220 CONTINUE 2236 210 CONTINUE 2237 200 CONTINUE 2238C 2239 NTOTA = MAX(NVIR(ISYMA),1) 2240 NTOCKJ = MAX(NCKI(ISYCKJ),1) 2241C 2242 KOFF1 = ISAIKJ(ISYCKJ,ISYMI) + 1 2243 KOFF2 = IT1AM(ISYMA,ISYMI) + A 2244C 2245 CALL DGEMV('T',NCKI(ISYCKJ),NRHF(ISYMI),-ONE,TMAT(KOFF1), 2246 * NTOCKJ,WORK,1,ONE,OMEGA1(KOFF2),NTOTA) 2247C 2248 CALL QEXIT('CC3_CY1') 2249C 2250 RETURN 2251 END 2252C 2253C /* Deck wbd_dia */ 2254 SUBROUTINE WBD_DIA(B,ISYMB,C,ISYMC,FREQY, 2255 * ISWMAT,WMAT,DIAG,FOCKD) 2256C 2257C Written by P. Jorgensen and F. Pawlowski, Spring 2002. 2258C 2259C Divide by orbital energies 2260C 2261C 2262 IMPLICIT NONE 2263C 2264#include "priunit.h" 2265#include "ccorb.h" 2266#include "ccsdsym.h" 2267C 2268 INTEGER ISYMB,ISYMC,ISWMAT,NB,NC 2269C 2270 REAL*8 WSMAT,DDOT,FREQY,EPSIBC 2271 REAL*8 WMAT(*), DIAG(*), FOCKD(*) 2272C 2273 CALL QENTER('WBD_DIA') 2274C 2275C----------------------------------------- 2276C Divide by the Fock matrix diagonals. 2277C----------------------------------------- 2278C 2279 NB = IORB(ISYMB) + NRHF(ISYMB) + B 2280 NC = IORB(ISYMC) + NRHF(ISYMC) + C 2281C 2282 EPSIBC = FOCKD(NB) + FOCKD(NC) - FREQY 2283C 2284 DO 500 L = 1,NCKIJ(ISWMAT) 2285C 2286 WMAT(L) = WMAT(L)/(DIAG(L) + EPSIBC) 2287C 2288 500 CONTINUE 2289C 2290c IF (IPRINT .GT. 55) THEN 2291c WSMAT = DDOT(NCKIJ(ISWMAT),WMAT,1,WMAT,1) 2292c WRITE(LUPRI,*) 'In wbd_dia: 5. Norm of WMAT ',WSMAT 2293c ENDIF 2294C 2295C--------------------- 2296C End. 2297C--------------------- 2298C 2299 CALL QEXIT('WBD_DIA') 2300C 2301 RETURN 2302 END 2303 2304C /* Deck wbdtest */ 2305 SUBROUTINE WBDtest(omega1,IB,IC, 2306 * TMAT,xiajb,iopt,work,lwrk) 2307C 2308C Written by P. Jorgensen and F. Pawlowski, Spring 2002. 2309C 2310C noddy code for the three cont to omega1 2311C 2312C 2313 IMPLICIT NONE 2314C 2315#include "priunit.h" 2316#include "ccorb.h" 2317#include "ccsdsym.h" 2318C 2319 INTEGER iopt,lwrk,ai,nck,nbj,kckbj,kaikj,kbjki,index,ib,ic 2320 integer kckji 2321C 2322 REAL*8 tMAT(*), xiajb(*), work(*), omega1(*) 2323C 2324 INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 2325C 2326c write(lupri,*)' entering wbdtest iop',iopt 2327 write(lupri,*)' tmat i test' 2328C call output(tmat,1,nvirt*nrhft,nrhft,nrhft,1,1, 2329C * nvirt*nrhft,nrhft,nrhft,1,lupri) 2330 do i = 1,nvirt*nrhft*nrhft*nrhft 2331 write(lupri,*) i , tmat(I) 2332 end do 2333 2334C 2335C----------------------------------------- 2336 call dcopy(nt1amx,omega1,1,work,1) 2337 write(lupri,*)' xi1eff wbdtest before cont is added ' 2338 do i = 1,nt1am(1) 2339 write(lupri,*) i , work(I) 2340 end do 2341C call dzero(work,nt1amx) 2342C 2343 if (iopt.eq.1) then 2344 b = ib 2345 c = ic 2346 write(lupri,*)' b,c third contribution',b,c 2347 DO AI = 1,NT1AMX 2348 DO K = 1,NRHFT 2349 DO J= 1,NRHFT 2350 nck = nvirt*(k-1) + c 2351 nbj = nvirt*(j-1) + b 2352 kckbj = index(nck,nbj) 2353 2354 kaikj = nvirt*nrhft*nrhft*(j-1) 2355 * + nvirt*nrhft*(k-1) + ai 2356c write(lupri,*)'nck,nbj,kckbj, xckbj' 2357c write(lupri,*)nck,nbj,kckbj,XIAJB(KCKBJ) 2358c write(lupri,*)'ai,j,k,kaikj,tmat(kaikj)' 2359c write(lupri,*)ai,j,k,kaikj,tmat(kaikj) 2360c write(lupri,*)'WORK(AI),XIAJB(KCKBJ)*tmat(kaikj' 2361c write(lupri,*)WORK(AI),XIAJB(KCKBJ)*tmat(kaikj) 2362 WORK(AI) = WORK(AI) + XIAJB(KCKBJ)*tmat(kaikj) 2363c write(lupri,*)'WORK(AI)',WORK(AI) 2364 end do 2365 end do 2366 end do 2367 2368 write(lupri,*)' xi1eff wbdtest first cont b,c',b,c 2369 do i = 1,nt1am(1) 2370 write(lupri,*) i , work(I) 2371 end do 2372 end if 2373cc 2374cc 2375 if (iopt.eq.2) then 2376 A = ic 2377 C = ib 2378 write(lupri,*)' a,c third contribution',a,c 2379 do i = 1,nrhft 2380 DO B = 1,nvirt 2381 ai = nvirt*(i-1)+a 2382 DO K = 1,NRHFT 2383 DO J= 1,NRHFT 2384 nck = nvirt*(k-1) + c 2385 nbj = nvirt*(j-1) + b 2386 kckbj = index(nck,nbj) 2387 2388 kbjki =nvirt*nrhft*nrhft*(i-1) 2389 * + nvirt*nrhft*(k-1) + nbj 2390c write(lupri,*)'nck,nbj,kckbj, xckbj' 2391c write(lupri,*)nck,nbj,kckbj,XIAJB(kckbj) 2392c WRITE(LUPRI,*)'AI,I,A' 2393c WRITE(LUPRI,*)AI,I,A 2394c write(lupri,*)'bj,j,k,kbjki,tmat(kbjki)' 2395c write(lupri,*)nbj,j,k,kbjki,tmat(kbjki) 2396c write(lupri,*)'WORK(AI),XIAJB(kckbj)*tmat(kbjki' 2397c write(lupri,*)WORK(AI),XIAJB(kckbj)*tmat(kbjki) 2398 WORK(AI) = WORK(AI) + XIAJB(kckbj)*tmat(kbjki) 2399c write(lupri,*)'WORK(AI)',WORK(AI) 2400 end do 2401 end do 2402 end do 2403 end do 2404 2405 write(lupri,*)' xi1eff wbdtest second cont a,c',a,c 2406 do i = 1,nt1am(1) 2407 write(lupri,*) i , work(I) 2408 end do 2409 end if 2410C 2411C 2412cc 2413 if (iopt.eq.3) then 2414 b = ic 2415 A = ib 2416 write(lupri,*)' a,b third contribution',a,b 2417 do i = 1,nrhft 2418 DO c = 1,nvirt 2419 ai = nvirt*(i-1)+a 2420 DO K = 1,NRHFT 2421 DO J= 1,NRHFT 2422 nck = nvirt*(k-1) + c 2423 nbj = nvirt*(j-1) + b 2424 kckbj = index(nck,nbj) 2425 2426 kckji =nvirt*nrhft*nrhft*(i-1) 2427 * + nvirt*nrhft*(j-1) + nck 2428c write(lupri,*)'nck,nbj,kckbj, xckbj' 2429c write(lupri,*)nck,nbj,kckbj,XIAJB(kckbj) 2430c WRITE(LUPRI,*)'AI,I,A' 2431c WRITE(LUPRI,*)AI,I,A 2432c write(lupri,*)'bj,j,k,kbjki,tmat(kbjki)' 2433c write(lupri,*)nbj,j,k,kbjki,tmat(kbjki) 2434c write(lupri,*)'WORK(AI),XIAJB(kckbj)*tmat(kbjki' 2435c write(lupri,*)WORK(AI),XIAJB(kckbj)*tmat(kbjki) 2436 WORK(AI) = WORK(AI) + XIAJB(kckbj)*tmat(kckji) 2437c write(lupri,*)'WORK(AI)',WORK(AI) 2438 end do 2439 end do 2440 end do 2441 end do 2442 2443 write(lupri,*)' xi1eff wbdtest third cont a,c',a,c 2444 do i = 1,nt1am(1) 2445 write(lupri,*) i , work(I) 2446 end do 2447 end if 2448 2449 RETURN 2450 END 2451C /* Deck CC3_CY2F */ 2452 SUBROUTINE CC3_CY2F(XI2,ISYRES,WMAT,ISWMAT,TMAT, 2453 * FOCK0,ISYINT,INDSQ,LENSQ,WORK,LWORK, 2454 * ISYMIB,IB,ISYMID,ID) 2455C 2456C P. Jorgensen and F. Pawlowski. Spring 2002 2457C 2458 IMPLICIT NONE 2459C 2460C Calculate Fock part to XI2 2461C 2462C 2463C General symmetry: ISWMAT is symmetry of WMAT and TMAT intermdiates. 2464C (exclusive isymb,isymd) 2465C ISYINT is symmetry of FOCK0 2466C ISYRES is symmetry of XI2 and XI2EFF 2467C 2468C Omega1(aibj) 2469C 2470C = sum_{ck} ( W^{abc}_{ijk} - W^{abc}_{kji} ) * F_{kc} 2471C (1): 2472C = sum_{k} ( W^{BC}_{aikj} - W^{BC}_{akij} ) * F_{kC} 2473C (2): 2474C + sum_{k} ( W^{CA}_{bjik} - W^{CA}_{bjki} ) * F_{kC} 2475C (3): 2476C + sum_{ck} ( W^{AB}_{ckji} - W^{AB}_{cijk} ) * F_{kc} 2477C 2478C 2479C 2480 INTEGER ISYRES, ISWMAT, ISYINT , LENSQ, LWORK, KSTART 2481 INTEGER ISYMIB, IB, ISYMID, ID, INDEX 2482 INTEGER INDSQ(LENSQ,6) 2483 INTEGER KFMAT, KABIJ, KEND, LEND 2484 INTEGER ISYMB, ISYMC, ISYMK, ISYAIJ, KOFFT, KOFFF, NTOAIJ 2485 INTEGER ISYMA, ISYBJI, NTOBJI, ISYAB, ISYIJ 2486 INTEGER KOFF1, KOFF2, ISYCK, NTOCK 2487 REAL*8 XI2(*), WMAT(*), TMAT(*), FOCK0(*), WORK(*) 2488 REAL*8 ZERO, ONE, TWO 2489 2490C 2491#include "priunit.h" 2492#include "ccsdsym.h" 2493#include "ccorb.h" 2494#include "ccsdinp.h" 2495C 2496 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 2497C 2498C INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 2499C 2500 CALL QENTER('CC3_CY2F') 2501 2502C 2503C work space allocations 2504C 2505 KSTART = 1 2506 KEND = KSTART + NVIRT*NRHFT*NRHFT 2507 LEND = LWORK - KEND 2508C 2509 IF (LEND .LT. 0) THEN 2510 WRITE(LUPRI,*) 'Memory available : ',LWORK 2511 WRITE(LUPRI,*) 'Memory needed : ',KEND 2512 CALL QUIT('Insufficient space in CC3_CY2F') 2513 END IF 2514C 2515C 2516C---------------------------------- 2517C (1): 2518C Omega2(aibj) = Omega2(aibj) 2519C + sum_{k} ( W^{BC}_{aikj} - W^{BC}_{akij} ) * F_{kC} 2520C---------------------------------- 2521C 2522C 2523 B = IB 2524 C = ID 2525 ISYMB = ISYMIB 2526 ISYMC = ISYMID 2527C 2528C T(aijk) = W^{BC}_{aikj} - W^{BC}_{akij} ) 2529C 2530 DO I = 1,LENSQ 2531C 2532 TMAT(I) = WMAT(INDSQ(I,3)) 2533 * - WMAT(INDSQ(I,4)) 2534C 2535 ENDDO 2536C 2537 ISYMK = MULD2H(ISYINT,ISYMC) 2538 ISYAIJ = MULD2H(ISYMK,ISWMAT) 2539 KOFFT = ISAIKJ(ISYAIJ,ISYMK) + 1 2540 KOFFF = IFCVIR(ISYMK,ISYMC) + NORB(ISYMK)*(C-1) + 1 2541 NTOAIJ = MAX(NCKI(ISYAIJ),1) 2542C 2543 IF (NRHF(ISYMK).GT.0) THEN 2544 CALL DGEMV('N',NCKI(ISYAIJ),NRHF(ISYMK),-ONE,TMAT(KOFFT), 2545 * NTOAIJ,FOCK0(KOFFF),1,ZERO,WORK,1) 2546C 2547 IF (NCKI(ISYAIJ).GT.0 ) THEN 2548 CALL CC3_RACC(XI2,WORK,ISYMB,B,ISYRES) 2549 END IF 2550 END IF 2551C 2552C---------------------------------- 2553C (2): 2554C Omega2(aibj) = Omega2(aibj) 2555C + sum_{k} ( W^{CA}_{bjik} - W^{CA}_{bjki} ) * F_{kC} 2556C---------------------------------- 2557C 2558 C = IB 2559 A = ID 2560C 2561 ISYMC = ISYMIB 2562 ISYMA = ISYMID 2563C 2564C 2565C T(bjik) = W^{CA}_{bjik} - W^{AC}_{bjki} ) 2566C 2567 DO I = 1,LENSQ 2568C 2569 TMAT(I) = WMAT(I) 2570 * - WMAT(INDSQ(I,3)) 2571C 2572 ENDDO 2573C 2574C 2575 ISYMK = MULD2H(ISYINT,ISYMC) 2576 ISYBJI = MULD2H(ISYMK,ISWMAT) 2577 KOFFT = ISAIKJ(ISYBJI,ISYMK) + 1 2578 KOFFF = IFCVIR(ISYMK,ISYMC) + NORB(ISYMK)*(C-1) + 1 2579 NTOBJI = MAX(NCKI(ISYBJI),1) 2580 IF (NRHF(ISYMK).GT.0) THEN 2581 CALL DGEMV('N',NCKI(ISYBJI),NRHF(ISYMK),-ONE,TMAT(KOFFT), 2582 * NTOBJI,FOCK0(KOFFF),1,ZERO,WORK,1) 2583C 2584 IF (NCKI(ISYBJI).GT.0) THEN 2585 CALL CC3_RACC(XI2,WORK,ISYMA,A,ISYRES) 2586 END IF 2587 END IF 2588C 2589C---------------------------------- 2590C (3): 2591C Omega2(aibj) = Omega2(aibj) 2592C + sum_{ck} ( W^{AB}_{ckji} - W^{AB}_{cijk} ) * F_{kc} 2593C---------------------------------- 2594C 2595 C = IB 2596C 2597 A = IB 2598 B = ID 2599C 2600 ISYMA = ISYMIB 2601 ISYMB = ISYMID 2602C 2603C 2604C T(ckij) = W^{AB}_{ckji} - W^{AB}_{cijk} 2605C 2606 DO I = 1,LENSQ 2607C 2608 TMAT(I) = WMAT(INDSQ(I,3)) 2609 * - WMAT(INDSQ(I,2)) 2610C 2611 ENDDO 2612C 2613 IF (NSYM .GT. 1) THEN 2614 CALL CC_GATHER(LENSQ,WORK,TMAT,INDSQ(1,6)) 2615 CALL DCOPY(LENSQ,WORK,1,TMAT,1) 2616 ENDIF 2617C 2618 2619C 2620C WORK SPACE ALLOCATION 2621C 2622 ISYAB = MULD2H(ISYMA,ISYMB) 2623 ISYIJ = MULD2H(ISYRES,ISYAB) 2624C 2625 KFMAT = 1 2626 KABIJ = KFMAT + NT1AM(ISYINT) 2627 KEND = KABIJ + NMATIJ(ISYIJ) 2628 LEND = LWORK - KEND 2629C 2630 IF (LEND .LT. 0) THEN 2631 WRITE(LUPRI,*) 'Memory available : ',LWORK 2632 WRITE(LUPRI,*) 'Memory needed : ',KEND 2633 CALL QUIT('Insufficient space in CC3_CY2F') 2634 END IF 2635C 2636C 2637C sort fock matrix 2638C 2639C FMAT(C,K) = FOCK0(K,C) 2640C 2641 DO ISYMC = 1,NSYM 2642 ISYMK = MULD2H(ISYMC,ISYINT) 2643 DO K = 1,NRHF(ISYMK) 2644 DO C = 1,NVIR(ISYMC) 2645 KOFF1 = IFCVIR(ISYMK,ISYMC) + NORB(ISYMK)*(C - 1) + K 2646 KOFF2 = KFMAT + IT1AM(ISYMC,ISYMK) 2647 * + NVIR(ISYMC)*(K - 1) + C -1 2648C 2649 WORK(KOFF2) = FOCK0(KOFF1) 2650C 2651 END DO 2652 END DO 2653 END DO 2654C 2655 ISYCK = ISYINT 2656 KOFFT = ISAIKL(ISYCK,ISYIJ) + 1 2657 NTOCK = MAX(NT1AM(ISYCK),1) 2658 IF (NMATIJ(ISYIJ).GT.0) THEN 2659C 2660 CALL DGEMV('T',NT1AM(ISYCK),NMATIJ(ISYIJ),-ONE,TMAT(KOFFT), 2661 * NTOCK,WORK(KFMAT),1,ZERO,WORK(KABIJ),1) 2662 CALL CC3_RABIJ(XI2,WORK(KABIJ),ISYMA,A,ISYMB,B,ISYRES) 2663C 2664 END IF 2665C 2666 CALL QEXIT('CC3_CY2F') 2667C 2668C 2669 RETURN 2670 END 2671C 2672C /* Deck CC3_RABIJ */ 2673 SUBROUTINE CC3_RABIJ(XI2,ABIJ,ISYMA,A,ISYMB,B,ISYRES) 2674C 2675C distribute matrix elements in ABIJ(I,J) to result vector XI2(AIBJ) 2676C 2677 IMPLICIT NONE 2678C 2679 INTEGER ISYMA, ISYMB, ISYRES, INDEX 2680 INTEGER ISYAB, ISYIJ, ISYMJ, ISYMI, ISYMBJ, ISYMAI 2681 INTEGER NBJ, NIJ, NAI, NAIBJ 2682 REAL*8 XI2(*), ABIJ(*) 2683 REAL*8 TWO 2684C 2685#include "priunit.h" 2686#include "ccsdsym.h" 2687#include "ccorb.h" 2688#include "ccsdinp.h" 2689C 2690 PARAMETER (TWO = 2.0D0) 2691C 2692 INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 2693C 2694 CALL QENTER('CC3_RABIJ') 2695C 2696 ISYAB = MULD2H(ISYMA,ISYMB) 2697 ISYIJ = MULD2H(ISYAB,ISYRES) 2698C 2699 DO 300 ISYMJ = 1,NSYM 2700 ISYMI = MULD2H(ISYIJ,ISYMJ) 2701 ISYMBJ = MULD2H(ISYMB,ISYMJ) 2702 ISYMAI = MULD2H(ISYMA,ISYMI) 2703C 2704 DO 310 J = 1,NRHF(ISYMJ) 2705 NBJ = IT1AM(ISYMB,ISYMJ) + NVIR(ISYMB)*(J - 1) + B 2706 IF (ISYMAI .EQ. ISYMBJ) THEN 2707C 2708 DO 320 I = 1,NRHF(ISYMI) 2709 NIJ = IMATIJ(ISYMI,ISYMJ) + NRHF(ISYMI)*(J - 1) + I 2710 NAI = IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I - 1) + A 2711 IF (NAI .EQ. NBJ) ABIJ(NIJ) = TWO*ABIJ(NIJ) 2712 NAIBJ = IT2AM(ISYMAI,ISYMBJ) + INDEX(NAI,NBJ) 2713 XI2(NAIBJ) = XI2(NAIBJ) + ABIJ(NIJ) 2714 320 CONTINUE 2715C 2716 ELSE IF (ISYMAI .LT. ISYMBJ) THEN 2717C 2718 DO 330 I = 1,NRHF(ISYMI) 2719 NIJ = IMATIJ(ISYMI,ISYMJ) + NRHF(ISYMI)*(J - 1) + I 2720 NAI = IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I - 1) + A 2721 NAIBJ = IT2AM(ISYMAI,ISYMBJ) 2722 * + NT1AM(ISYMAI)*(NBJ-1) + NAI 2723 XI2(NAIBJ) = XI2(NAIBJ) + ABIJ(NIJ) 2724 330 CONTINUE 2725C 2726 ELSE IF (ISYMBJ .LT. ISYMAI) THEN 2727C 2728 DO 340 I = 1,NRHF(ISYMI) 2729 NIJ = IMATIJ(ISYMI,ISYMJ) + NRHF(ISYMI)*(J - 1) + I 2730 NAI = IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I - 1) + A 2731 NAIBJ = IT2AM(ISYMBJ,ISYMAI) 2732 * + NT1AM(ISYMBJ)*(NAI-1) + NBJ 2733 XI2(NAIBJ) = XI2(NAIBJ) + ABIJ(NIJ) 2734 340 CONTINUE 2735C 2736 ENDIF 2737 310 CONTINUE 2738C 2739 300 CONTINUE 2740C 2741 CALL QEXIT('CC3_RABIJ') 2742C 2743 RETURN 2744 END 2745C /* Deck cc3_cy2o */ 2746 SUBROUTINE CC3_CY2O(XI2,ISYRES,WMAT,ISWMAT,TMAT,TROCC, 2747 * TROCC1,ISYINT,WORK,LWORK,INDSQ,LENSQ, 2748 * ISYMIB,IB,ISYMID,ID) 2749C 2750C 2751C General symmetry: ISWMAT is symmetry of WMAT and TMAT intermediates. 2752C (exclusive isymib,isymid) 2753C ISYINT is symmetry of integrals in TROCC and TROCC1. 2754C ISYRES = ISWMAT*ISYINT*ISYMIB*ISYMID 2755C 2756C XI2(aibj) = XI2(aibj) 2757C 2758C + sum_{ckl} (2W^{abc}_{ikl}-W^{abc}_{lki}-W^{abc}_{ilk}) * (lc|kj) 2759C (1): 2760C = sum_{kl} (2W^{BC}_{ailk}-W^{BC}_{alik}-W^{BC}_{aikl}) * (lc|kj) 2761C (2): 2762C + sum_{kl} (2W^{CA}_{bkil}-W^{CA}_{bkli}-W^{CA}_{blik}) * (lc|kj) 2763C (3): 2764C + sum_{ckl} (2W^{AB}_{clki}-W^{AB}_{cikl}-W^{AB}_{ckli}) * (lc|kj) 2765C 2766 IMPLICIT NONE 2767C 2768 INTEGER ISYRES, ISWMAT, ISYINT, LWORK, LENSQ 2769 INTEGER ISYMIB, IB, ISYMID, ID 2770 INTEGER INDSQ(LENSQ,6) 2771 INTEGER ISYMC, ISYMB, ISYBC, JSAIKL, LENGTH, ISYMJ, ISYBJ 2772 INTEGER ISYAI, ISYKL, NTOTAI 2773 INTEGER NTOTKL, KOFF1, KOFF2, KOFF3 2774 INTEGER ISYMA, ISYAB, JSCKLI, JSBIKL, ISYMI 2775 INTEGER ISYCKL, NTOCKL, NRHFI 2776 INTEGER KRMAT1, KRMAT2, KEND, LEND, INDEX 2777 INTEGER ISYCK, ISYCJ, ISYAC, ISYBI, ISYKLJ, NTOTBI 2778 INTEGER ISWB, ISWD, NCKIMX 2779 INTEGER ISYAIJ, ISYBJI, ISYIJ 2780 INTEGER KTMATX 2781C 2782 REAL*8 XI2(*), WMAT(*), TMAT(*), TROCC(*), TROCC1(*) 2783 REAL*8 WORK(*) 2784 REAL*8 TWO, ONE, ZERO, XTMAT, XRMAT, DDOT 2785C 2786#include "priunit.h" 2787#include "ccsdsym.h" 2788#include "ccorb.h" 2789#include "ccsdinp.h" 2790C 2791 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 2792C 2793C INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 2794C 2795 CALL QENTER('CC3_CY2O') 2796C 2797 2798 ISWB = MULD2H(ISWMAT,ISYMIB) 2799 ISWD = MULD2H(ISWMAT,ISYMID) 2800 ISYAIJ = MULD2H(ISYMIB,ISYRES) 2801 ISYBJI = MULD2H(ISYMID,ISYRES) 2802 NCKIMX = MAX(NCKI(ISWB),NCKI(ISWD),NCKI(ISYAIJ),NCKI(ISYBJI)) 2803C 2804 KRMAT1 = 1 2805 KRMAT2 = KRMAT1 + NCKIMX 2806 KTMATX = KRMAT2 + NCKIMX 2807 KEND = KTMATX + NCKIJ(ISWMAT) 2808 LEND = LWORK - KEND 2809 IF (LWORK .LT. KEND) THEN 2810 CALL QUIT('Insufficient space in CY2O') 2811 ENDIF 2812C 2813C------------------------- 2814C First occupied term. 2815C------------------------- 2816C (1): 2817C XI2(aibj) = XI2(aibj) 2818C 2819C + sum_{kl} (2W^{BC}_{ailk}-W^{BC}_{alik}-W^{BC}_{aikl}) * (lc|kj) 2820 2821C 2822 B = IB 2823 C = ID 2824C 2825 ISYMB = ISYMIB 2826 ISYMC = ISYMID 2827C 2828 ISYAIJ = MULD2H(ISYMB,ISYRES) 2829 IF (NCKI(ISYAIJ).GT.0 ) THEN 2830 CALL DZERO(WORK(KRMAT1),NCKI(ISYAIJ)) 2831C 2832 ISYBC = MULD2H(ISYMB,ISYMC) 2833 JSAIKL = ISWMAT 2834C 2835 LENGTH = NCKIJ(JSAIKL) 2836C 2837C--------------------------------------- 2838C 2839C TMAT(aikl) = (2W^{BC}_{ailk}-W^{BC}_{alik}-W^{BC}_{aikl}) 2840C 2841C--------------------------------------- 2842C 2843 DO I = 1,LENGTH 2844C 2845 TMAT(I) = TWO*WMAT(INDSQ(I,3)) 2846 * - WMAT(INDSQ(I,4)) 2847 * - WMAT(I) 2848C 2849 ENDDO 2850C 2851C---------------------------------- 2852C Symmetry sorting if symmetry. 2853C---------------------------------- 2854C 2855 2856 IF (NSYM .GT. 1) THEN 2857 CALL CC_GATHER(LENGTH,WORK(KTMATX),TMAT,INDSQ(1,6)) 2858 CALL DCOPY(LENGTH,WORK(KTMATX),1,TMAT,1) 2859 ENDIF 2860C 2861 IF (IPRINT .GT. 55) THEN 2862 XTMAT = DDOT(NCKIJ(JSAIKL),TMAT,1,TMAT,1) 2863 WRITE(LUPRI,*) 'In CC3_CY2O: 1. Norm of TMAT = ',XTMAT 2864 ENDIF 2865C 2866C----------------------- 2867C First contraction. 2868C----------------------- 2869C 2870C TROCC(kljC) = (lc|kj) 2871C 2872C XI2(aibj) = XI2(aibj) 2873C + sum(kl) TMAT(aikl)* TROCC(kljC) 2874C 2875 DO 200 ISYMJ = 1,NSYM 2876C 2877 ISYBJ = MULD2H(ISYMB,ISYMJ) 2878 ISYAI = MULD2H(ISYBJ,ISYRES) 2879 ISYKL = MULD2H(JSAIKL,ISYAI) 2880 ISYKLJ = MULD2H(ISYKL,ISYMJ) 2881C 2882 NTOTAI = MAX(NT1AM(ISYAI),1) 2883 NTOTKL = MAX(NMATIJ(ISYKL),1) 2884C 2885 KOFF1 = ISAIKL(ISYAI,ISYKL) + 1 2886 KOFF2 = ISJIKA(ISYKLJ,ISYMC) 2887 * + NMAJIK(ISYKLJ)*(C - 1) 2888 * + ISJIK(ISYKL,ISYMJ) + 1 2889 KOFF3 = KRMAT1 + ISAIK(ISYAI,ISYMJ) 2890C 2891 CALL DGEMM('N','N',NT1AM(ISYAI),NRHF(ISYMJ),NMATIJ(ISYKL), 2892 * ONE,TMAT(KOFF1),NTOTAI,TROCC(KOFF2),NTOTKL, 2893 * ONE,WORK(KOFF3),NTOTAI) 2894C 2895 200 CONTINUE 2896C 2897 IF (IPRINT .GT. 55) THEN 2898 XRMAT = DDOT(NCKI(ISYAIJ),WORK(KRMAT1),1,WORK(KRMAT1),1) 2899 WRITE(LUPRI,*) '1. In CC3_CY2O: Norm of RMAT1 = ',XRMAT 2900 ENDIF 2901C 2902 CALL CC3_RACC(XI2,WORK(KRMAT1),ISYMB,B,ISYRES) 2903C 2904 END IF 2905C 2906C-------------------------- 2907C Second occupied term. 2908C-------------------------- 2909C 2910C (2): 2911C 2912C XI2(aibj) = XI2(aibj) 2913C 2914C + sum_{kl} (2W^{CA}_{bkil}-W^{CA}_{bkli}-W^{CA}_{blik}) * (lc|kj) 2915C 2916 A = ID 2917 C = IB 2918C 2919 ISYMA = ISYMID 2920 ISYMC = ISYMIB 2921C 2922 ISYBJI = MULD2H(ISYMA,ISYRES) 2923 IF (NCKI(ISYBJI).GT.0) THEN 2924 CALL DZERO(WORK(KRMAT1),NCKI(ISYBJI)) 2925C 2926 ISYAC = MULD2H(ISYMA,ISYMC) 2927 JSBIKL = ISWMAT 2928C 2929C--------------------------------------- 2930C 2931C TMAT(bikl) = (2W^{CA}_{bkil}-W^{CA}_{bkli}-W^{CA}_{blik}) 2932C 2933C--------------------------------------- 2934C 2935 LENGTH = NCKIJ(JSBIKL) 2936C 2937 DO I = 1,LENGTH 2938C 2939 TMAT(I) = TWO*WMAT(INDSQ(I,1)) 2940 * - WMAT(INDSQ(I,2)) 2941 * - WMAT(INDSQ(I,4)) 2942C 2943 ENDDO 2944C 2945C---------------------------------- 2946C Symmetry sorting if symmetry. 2947C---------------------------------- 2948C 2949 IF (NSYM .GT. 1) THEN 2950 CALL CC_GATHER(LENGTH,WORK(KTMATX),TMAT,INDSQ(1,6)) 2951 CALL DCOPY(LENGTH,WORK(KTMATX),1,TMAT,1) 2952 ENDIF 2953C 2954 IF (IPRINT .GT. 55) THEN 2955 XTMAT = DDOT(NCKIJ(JSAIKL),TMAT,1,TMAT,1) 2956 WRITE(LUPRI,*) 'In CC3_CY2O: 2. Norm of TMAT = ',XTMAT 2957 ENDIF 2958C 2959C------------------------ 2960C Second contraction. 2961C------------------------ 2962C 2963C TROCC(kljC) = (lc|kj) 2964C 2965C M^A_(bij) = sum(kl) TMAT(bikl)* TTROCC(kljC) 2966C 2967 2968 DO 400 ISYMJ = 1,NSYM 2969C 2970 ISYCJ = MULD2H(ISYMC,ISYMJ) 2971 ISYKL = MULD2H(ISYINT,ISYCJ) 2972 ISYBI = MULD2H(ISYKL,ISWMAT) 2973 ISYKLJ = MULD2H(ISYKL,ISYMJ) 2974C 2975 NTOTBI = MAX(NT1AM(ISYBI),1) 2976 NTOTKL = MAX(NMATIJ(ISYKL),1) 2977C 2978 KOFF1 = ISAIKL(ISYBI,ISYKL) + 1 2979 KOFF2 = ISJIKA(ISYKLJ,ISYMC) 2980 * + NMAJIK(ISYKLJ)*(C - 1) 2981 * + ISJIK(ISYKL,ISYMJ) + 1 2982 KOFF3 = KRMAT1 + ISAIK(ISYBI,ISYMJ) 2983C 2984 2985 CALL DGEMM('N','N',NT1AM(ISYBI),NRHF(ISYMJ),NMATIJ(ISYKL), 2986 * ONE,TMAT(KOFF1),NTOTBI,TROCC(KOFF2),NTOTKL, 2987 * ONE,WORK(KOFF3),NTOTBI) 2988C 2989 400 CONTINUE 2990C 2991 2992 IF (IPRINT .GT. 55) THEN 2993 XRMAT = DDOT(NCKI(ISYBJI),WORK(KRMAT1),1,WORK(KRMAT1),1) 2994 WRITE(LUPRI,*) '2.In CC3_CY2O: Norm of RMAT1 = ',XRMAT 2995 ENDIF 2996C 2997C reorder result vector M^A_(bij) as M^A_(bji) 2998C 2999 CALL CC3_MAJI(WORK(KRMAT1),WORK(KRMAT2),ISYMA,A,ISYRES) 3000C 3001 IF (IPRINT .GT. 55) THEN 3002 XRMAT = DDOT(NCKI(ISYBJI),WORK(KRMAT2),1,WORK(KRMAT2),1) 3003 WRITE(LUPRI,*) 3004 * 'In CC3_CY2O: Reorder :Norm of RMAT2 = ',XRMAT 3005 ENDIF 3006C 3007C add vector to XI2 3008C 3009 CALL CC3_RACC(XI2,WORK(KRMAT2),ISYMA,A,ISYRES) 3010C 3011 END IF 3012C 3013C------------------------- 3014C Third occupied term. 3015C------------------------- 3016C (3): 3017C XI2(aibj) = XI2(aibj) 3018C 3019C + sum_{ckl} (2W^{AB}_{clki}-W^{AB}_{cikl}-W^{AB}_{ckli}) * (lc|kj) 3020C 3021 B = ID 3022 A = IB 3023C 3024 ISYMB = ISYMID 3025 ISYMA = ISYMIB 3026C 3027 ISYAB = MULD2H(ISYMA,ISYMB) 3028 ISYIJ = MULD2H(ISYAB,ISYRES) 3029 IF (NMATIJ(ISYIJ).GT.0) THEN 3030 CALL DZERO(WORK,NMATIJ(ISYIJ)) 3031C 3032 JSCKLI = ISWMAT 3033C 3034 LENGTH = NCKIJ(JSCKLI) 3035C 3036C---------------------------------- 3037C 3038C TMAT(ckli) = 2W^{AB}_{clki}-W^{AB}_{cikl}-W^{AB}_{ckli} 3039C 3040C---------------------------------- 3041C 3042 DO I = 1,LENGTH 3043C 3044 TMAT(I) = TWO*WMAT(INDSQ(I,1)) 3045 * - WMAT(INDSQ(I,4)) 3046 * - WMAT(I) 3047C 3048 ENDDO 3049C 3050C 3051 IF (IPRINT .GT. 55) THEN 3052 XTMAT = DDOT(NCKIJ(JSAIKL),TMAT,1,TMAT,1) 3053 WRITE(LUPRI,*) 'In CC3_CY2O: 3. Norm of TMAT = ',XTMAT 3054 ENDIF 3055C 3056C----------------------- 3057C Third contraction. 3058C----------------------- 3059C 3060C 3061C TROCC1(cklj) = (lc|kj) 3062C 3063C XI2(aibj) = XI2(aibj) 3064C + sum(kl) TMAT(ckli)* TROCC1(cklj) 3065C 3066C 3067 DO 600 ISYMJ = 1,NSYM 3068C 3069 ISYBJ = MULD2H(ISYMB,ISYMJ) 3070 ISYAI = MULD2H(ISYBJ,ISYRES) 3071 ISYMI = MULD2H(ISYAI,ISYMA) 3072 ISYCKL = MULD2H(ISYMI,JSCKLI) 3073C 3074 IF (LWORK .LT. NRHF(ISYMI)*NRHF(ISYMJ)) THEN 3075 CALL QUIT('Insufficient memory in CCSDT_CY2O') 3076 END IF 3077C 3078 NTOCKL = MAX(NCKI(ISYCKL),1) 3079 NRHFI = MAX(NRHF(ISYMI),1) 3080C 3081 KOFF1 = ISAIKJ(ISYCKL,ISYMI) + 1 3082 KOFF2 = ISAIKJ(ISYCKL,ISYMJ) + 1 3083 KOFF3 = IMATIJ(ISYMI,ISYMJ) + 1 3084C 3085 3086 CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),NCKI(ISYCKL), 3087 * ONE,TMAT(KOFF1),NTOCKL,TROCC1(KOFF2),NTOCKL, 3088 * ONE,WORK(KOFF3),NRHFI) 3089C 3090C 3091 3092 600 CONTINUE 3093C 3094 IF (IPRINT .GT. 55) THEN 3095 XRMAT = DDOT(NMATIJ(ISYIJ),WORK,1,WORK,1) 3096 WRITE(LUPRI,*) '3.In CC3_CY2O: Norm of RABIJ = ',XRMAT 3097 ENDIF 3098C 3099 CALL CC3_RABIJ(XI2,WORK,ISYMA,A,ISYMB,B,ISYRES) 3100C 3101 END IF 3102C 3103 CALL QEXIT('CC3_CY2O') 3104C 3105 RETURN 3106 END 3107C /* Deck CC3_MAJI */ 3108 SUBROUTINE CC3_MAJI(RMAT1,RMAT2,ISYMB,B,ISYRES) 3109C 3110C order matrix M^B_(aji) as M^B_(aij) 3111C 3112 IMPLICIT NONE 3113C 3114 INTEGER ISYMB,ISYRES 3115 INTEGER ISYAJI, ISYMI, ISYAJ, ISYMJ, ISYMA, ISYAI 3116 INTEGER KOFF1, KOFF2, KMAT, KEND0, LWRK0 3117 REAL*8 RMAT1(*), RMAT2(*) 3118C 3119#include "priunit.h" 3120#include "ccsdsym.h" 3121#include "ccorb.h" 3122#include "ccsdinp.h" 3123 3124C 3125 CALL QENTER('CC3_MAJI') 3126C 3127 ISYAJI = MULD2H(ISYRES,ISYMB) 3128 DO ISYMI = 1,NSYM 3129 ISYAJ = MULD2H(ISYMI,ISYAJI) 3130 DO I = 1,NRHF(ISYMI) 3131 DO ISYMJ = 1,NSYM 3132 ISYMA = MULD2H(ISYAJ,ISYMJ) 3133 ISYAI = MULD2H(ISYMA,ISYMI) 3134 DO J = 1,NRHF(ISYMJ) 3135 DO A = 1,NVIR(ISYMA) 3136 KOFF1 = ISAIK(ISYAJ,ISYMI) 3137 * + NT1AM(ISYAJ)*(I-1) 3138 * + IT1AM(ISYMA,ISYMJ) 3139 * + NVIR(ISYMA)*(J-1) + A 3140 KOFF2 = ISAIK(ISYAI,ISYMJ) 3141 * + NT1AM(ISYAI)*(J-1) 3142 * + IT1AM(ISYMA,ISYMI) 3143 * + NVIR(ISYMA)*(I-1) + A 3144 3145 RMAT2(KOFF2) = RMAT1(KOFF1) 3146 END DO 3147 END DO 3148 END DO 3149 END DO 3150 END DO 3151C 3152 CALL QEXIT('CC3_MAJI') 3153C 3154 3155 RETURN 3156 END 3157C /* Deck cc3_convir */ 3158 SUBROUTINE CC3_CY2V(XI2,ISYRES,RBJIA,WMAT,ISWMAT,TMAT,TRVIR, 3159 * TRVIR1,ISYINT,WORK,LWORK,INDSQ,LENSQ, 3160 * ISYMIB,IB,ISYMID,ID) 3161C 3162C 3163C General symmetry: ISWMAT is symmetry of WMAT and TMAT intermediates. 3164C (exclusive isymib,isymid) 3165C ISYINT is symmetry of integrals in TRVIR and TROVIR1. 3166C ISYRES = ISWMAT*ISYINT*ISYMIB*ISYMID 3167C RBJIA intermediate result vector 3168C 3169C XI2(aibj) = XI2(aibj) 3170C 3171C + sum_{cdk} (2W^{bcd}_{jik}-W^{bcd}_{kij}-W^{bcd}_{jki}) * (ac|kd) 3172C (1): 3173C = sum_{ck} (2W^{DB}_{cijk}-W^{DB}_{cikj}-W^{DB}_{ckji}) * (ac|kD) 3174C (2): 3175C + sum_{dk} (2W^{BC}_{dkij}-W^{BC}_{djik}-W^{BC}_{dikj}) * (aC|kd) 3176C (3): 3177C + sum_{k} (2W^{CD}_{bjki}-W^{CD}_{bkji}-W^{CD}_{bjik}) * (aC|kD) 3178C 3179 IMPLICIT NONE 3180C 3181 INTEGER ISYRES, ISWMAT, ISYINT , LENSQ, LWORK 3182 INTEGER INDSQ(LENSQ,6) 3183 INTEGER ISYMIB, IB, ISYMID, ID 3184 INTEGER INDEX ,LENGTH, ISCKIJ, ISYMB, ISYAIJ, KRMAT, KEND, LEND 3185 INTEGER ISYMJ ,ISYBJ, ISYAI, ISYCKI, ISYCK, ISYMA, NTOTCK, NVIRA 3186 INTEGER KOFF1 , KOFF2, KOFF3, ISYMD, ISDKIJ, ISYDKI 3187 INTEGER ISYDK, NTOTDK, ISBJIK, ISYCKA 3188 INTEGER ISYKA, KINT, ISYBJI, KOFFT, KOFFM, KOFFR 3189 INTEGER NTOK, NTOBJI, ISYMK, ISYMC, ISYMI 3190 REAL*8 XI2(*), RBJIA(*), WMAT(*), TMAT(*), TRVIR(*), TRVIR1(*) 3191 REAL*8 WORK(*) 3192 REAL*8 ZERO, ONE, TWO 3193 3194C 3195 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 3196C 3197#include "priunit.h" 3198#include "ccsdsym.h" 3199#include "ccorb.h" 3200#include "ccsdinp.h" 3201C 3202C 3203C INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 3204C 3205 CALL QENTER('CC3_CY2V') 3206 LENGTH = NCKIJ(ISWMAT) 3207C 3208C------------------------ 3209C (1): 3210C = sum_{ck} (2W^{DB}_{cijk}-W^{DB}_{cikj}-W^{DB}_{ckji}) * (ac|kD) 3211C 3212C------------------------ 3213C 3214C TMAT(ckij) = 2W^{DB}_{cijk}-W^{DB}_{cikj}-W^{DB}_{ckji} 3215C 3216C Use here the symmetry between BJ, DK !! 3217C 3218 ISCKIJ = ISWMAT 3219C 3220 ISYMB = ISYMIB 3221 ISYMD = ISYMID 3222 B = IB 3223 D = ID 3224C 3225 ISYAIJ = MULD2H(ISYMB,ISYRES) 3226 KRMAT = 1 3227 KEND = KRMAT + NCKI(ISYAIJ) 3228 LEND = LWORK - KEND 3229 IF (LEND .LT. 0)THEN 3230 CALL QUIT('1. Insufficient core in CC3_CY2V') 3231 ENDIF 3232 CALL DZERO(WORK(KRMAT), NCKI(ISYAIJ)) 3233 DO I = 1,LENGTH 3234 TMAT(I) = TWO*WMAT(INDSQ(I,1)) 3235 * - WMAT(INDSQ(I,2)) 3236 * - WMAT(I) 3237 ENDDO 3238C 3239C--------------------------- 3240C (ac|kD) = TRVIR^D_{ck,a} 3241C 3242C RMAT^B_(aij) = TRVIR^D_{ck,a} )^T * TMAT(ckij) 3243C--------------------------- 3244C 3245 DO 200 ISYMJ = 1,NSYM 3246C 3247 ISYBJ = MULD2H(ISYMB,ISYMJ) 3248 ISYAI = MULD2H(ISYBJ,ISYRES) 3249 ISYCKI = MULD2H(ISCKIJ,ISYMJ) 3250C 3251 DO 210 J = 1,NRHF(ISYMJ) 3252C 3253 DO 220 ISYMI = 1,NSYM 3254C 3255 ISYCK = MULD2H(ISYCKI,ISYMI) 3256 ISYMA = MULD2H(ISYAI,ISYMI) 3257C 3258 NTOTCK = MAX(NT1AM(ISYCK),1) 3259 NVIRA = MAX(NVIR(ISYMA),1) 3260C 3261 KOFF1 = ICKATR(ISYCK,ISYMA) + 1 3262 KOFF2 = ISAIKJ(ISYCKI,ISYMJ) 3263 * + NCKI(ISYCKI)*(J - 1) 3264 * + ISAIK(ISYCK,ISYMI) + 1 3265 KOFF3 = KRMAT 3266 * + ISAIK(ISYAI,ISYMJ) 3267 * + NT1AM(ISYAI)*(J - 1) 3268 * + IT1AM(ISYMA,ISYMI) 3269C 3270 CALL DGEMM('T','N',NVIR(ISYMA),NRHF(ISYMI),NT1AM(ISYCK), 3271 * -ONE,TRVIR(KOFF1),NTOTCK,TMAT(KOFF2),NTOTCK, 3272 * ONE,WORK(KOFF3),NVIRA) 3273C 3274 220 CONTINUE 3275C 3276 210 CONTINUE 3277C 3278 200 CONTINUE 3279C 3280 ISYAIJ = MULD2H(ISYMB,ISYRES) 3281 IF ( NCKI(ISYAIJ).GT.0 ) THEN 3282 CALL CC3_RACC(XI2,WORK(KRMAT),ISYMB,B,ISYRES) 3283 END IF 3284C 3285C 3286C (2): 3287C 3288C XI2(aibj) = XI2(aibj) 3289C 3290C + sum_{dk} (2W^{BC}_{dkij}-W^{BC}_{djik}-W^{BC}_{dikj}) * (aC|kd) 3291C 3292C TMAT(dkij) = 2W^{BC}_{dkij}-W^{BC}_{djik}-W^{BC}_{dikj} 3293C 3294 ISYMC = ISYMID 3295 ISYMB = ISYMIB 3296 B = IB 3297 C = ID 3298C 3299 ISDKIJ = ISWMAT 3300 ISYAIJ = MULD2H(ISYMB,ISYRES) 3301 KRMAT = 1 3302 KEND = KRMAT + NCKI(ISYAIJ) 3303 LEND = LWORK - KEND 3304 IF (LEND .LT. 0)THEN 3305 CALL QUIT('2. Insufficient core in CC3_CY2V') 3306 ENDIF 3307 CALL DZERO(WORK(KRMAT), NCKI(ISYAIJ)) 3308C 3309 DO I = 1,LENGTH 3310 TMAT(I) = TWO*WMAT(I) 3311 * - WMAT(INDSQ(I,5)) 3312 * - WMAT(INDSQ(I,1)) 3313 ENDDO 3314C 3315C (ad|kC) = TRVIR1^C(dka) 3316C 3317C RMAT^B_(aij) = ( TRVIR1^C(dka) )^T * TMAT(dkij) 3318C 3319 DO 300 ISYMJ = 1,NSYM 3320C 3321 ISYBJ = MULD2H(ISYMB,ISYMJ) 3322 ISYAI = MULD2H(ISYBJ,ISYRES) 3323 ISYDKI = MULD2H(ISDKIJ,ISYMJ) 3324C 3325 DO 310 J = 1,NRHF(ISYMJ) 3326C 3327 DO 320 ISYMI = 1,NSYM 3328C 3329 ISYDK = MULD2H(ISYDKI,ISYMI) 3330 ISYMA = MULD2H(ISYAI,ISYMI) 3331C 3332 NTOTDK = MAX(NT1AM(ISYDK),1) 3333 NVIRA = MAX(NVIR(ISYMA),1) 3334C 3335 KOFF1 = ICKATR(ISYDK,ISYMA) + 1 3336 KOFF2 = ISAIKJ(ISYDKI,ISYMJ) 3337 * + NCKI(ISYDKI)*(J - 1) 3338 * + ISAIK(ISYDK,ISYMI) + 1 3339 KOFF3 = KRMAT + ISAIK(ISYAI,ISYMJ) 3340 * + NT1AM(ISYAI)*(J - 1) 3341 * + IT1AM(ISYMA,ISYMI) 3342C 3343 CALL DGEMM('T','N',NVIR(ISYMA),NRHF(ISYMI),NT1AM(ISYDK), 3344 * -ONE,TRVIR1(KOFF1),NTOTDK,TMAT(KOFF2),NTOTDK, 3345 * ONE,WORK(KOFF3),NVIRA) 3346C 3347 320 CONTINUE 3348 310 CONTINUE 3349 300 CONTINUE 3350C 3351 ISYAIJ = MULD2H(ISYMB,ISYRES) 3352 IF ( NCKI(ISYAIJ).GT.0 ) THEN 3353 CALL CC3_RACC(XI2,WORK(KRMAT),ISYMB,B,ISYRES) 3354 END IF 3355C 3356C----------------------------------- 3357C (3): 3358C + sum_{k} (2W^{CD}_{bjki}-W^{CD}_{bkji}-W^{CD}_{bjik}) * (ac|kD) 3359C 3360C----------------------------------- 3361C 3362 ISYMC = ISYMIB 3363 ISYMD = ISYMID 3364 C = IB 3365 D = ID 3366C 3367C TMAT(bjik) = 2W^{CD}_{bjki}-W^{CD}_{bkji}-W^{CD}_{bjik} 3368C 3369 ISBJIK = ISWMAT 3370C 3371 DO I = 1,LENGTH 3372 TMAT(I) = TWO*WMAT(INDSQ(I,3)) 3373 * - WMAT(INDSQ(I,4)) 3374 * - WMAT(I) 3375 ENDDO 3376C 3377C (ac|kD) = TRVIR^D_{Cka} = M^DC_{ka} 3378C 3379c 3380 ISYCKA = MULD2H(ISYINT,ISYMD) 3381 ISYKA = MULD2H(ISYCKA,ISYMC) 3382C 3383 DO ISYMA = 1,NSYM 3384 ISYCK = MULD2H(ISYCKA,ISYMA) 3385 ISYMK = MULD2H(ISYMC,ISYCK) 3386 DO A = 1,NVIR(ISYMA) 3387 DO K = 1,NRHF(ISYMK) 3388 KOFFM = IT1AMT(ISYMK,ISYMA) 3389 * + NRHF(ISYMK)*(A-1) + K 3390 KINT = ICKATR(ISYCK,ISYMA) 3391 * + NT1AM(ISYCK)*(A-1) 3392 * + IT1AM(ISYMC,ISYMK) 3393 * + NVIR(ISYMC)*(K-1) + C 3394 WORK(KOFFM) = TRVIR(KINT) 3395 END DO 3396 END DO 3397 END DO 3398C 3399C RBJIA(bjia) = RBJIA(bjia) + TMAT(bjik) * M^DC_{ka} 3400C 3401 DO ISYMA = 1,NSYM 3402 ISYMK = MULD2H(ISYMA,ISYKA) 3403 ISYBJI = MULD2H(ISBJIK,ISYMK) 3404 KOFFT = ISAIKJ(ISYBJI,ISYMK) + 1 3405 KOFFM = IT1AMT(ISYMK,ISYMA) + 1 3406 KOFFR = IT2SP(ISYBJI,ISYMA) + 1 3407 NTOK = MAX(NRHF(ISYMK),1) 3408 NTOBJI = MAX(NCKI(ISYBJI),1) 3409C 3410 CALL DGEMM('N','N',NCKI(ISYBJI),NVIR(ISYMA), 3411 * NRHF(ISYMK),-ONE,TMAT(KOFFT),NTOBJI, 3412 * WORK(KOFFM),NTOK,ONE,RBJIA(KOFFR),NTOBJI) 3413C 3414 END DO 3415C 3416 CALL QEXIT('CC3_CY2V') 3417 RETURN 3418 END 3419C /* Deck CC3_RBJIA */ 3420 SUBROUTINE CC3_RBJIA(XI2,ISYRES,RBJIA) 3421C 3422C distribute RBJIA in result vector XI2 3423C XI2(ai,bj) = XI2(ai,bj) + RBJIA(bjai) + RBJIA(aibj) 3424C 3425C XI2 triangular packed 3426C RBJIA squared packed 3427C 3428 IMPLICIT NONE 3429C 3430 INTEGER ISYMA,ISYRES,ISYBJI,KOFFR 3431C 3432 REAL*8 XI2(*), RBJIA(*) 3433C 3434#include "priunit.h" 3435#include "ccsdsym.h" 3436#include "ccorb.h" 3437#include "ccsdinp.h" 3438 3439C 3440 CALL QENTER('CC3_RBJIA') 3441C 3442 DO ISYMA = 1,NSYM 3443 ISYBJI = MULD2H(ISYRES,ISYMA) 3444 DO A = 1,NVIR(ISYMA) 3445 KOFFR = IT2SP(ISYBJI,ISYMA) 3446 * + NCKI(ISYBJI)*(A-1) + 1 3447 CALL CC3_RACC(XI2,RBJIA(KOFFR),ISYMA,A,ISYRES) 3448 END DO 3449 END DO 3450C 3451 CALL QEXIT('CC3_RBJIA') 3452C 3453 RETURN 3454 END 3455C 3456C /* Deck get_filename */ 3457 SUBROUTINE GET_FILENAME(FNCKJDR,FNDKBCR,MCAU) 3458* 3459************************************************************ 3460* Based on MCAU value get the name of the files for C1-transformed 3461* integrals in CC3 Cauchy vector calculations 3462* 3463* Currently MCAU may not be larger than 999 3464* 3465* Filip Pawlowski, 11-May-2004, Aarhus 3466************************************************************ 3467* 3468 IMPLICIT NONE 3469C 3470#include "priunit.h" 3471C 3472 LOGICAL LCDBG 3473 PARAMETER (LCDBG = .FALSE.) 3474C 3475 INTEGER MCAU 3476 INTEGER UNITS,TENTH,HUNDREDS 3477 CHARACTER*(*) FNCKJDR,FNDKBCR 3478 CHARACTER*3 NFIL 3479C 3480 CALL QENTER('GETFN') 3481C 3482 !Convert MCAU to character string NFIL 3483 IF ((MCAU.GE.0) .AND. (MCAU.LE.9)) THEN 3484C 3485 NFIL = '_'//'_'//CHAR(MCAU+48) 3486C 3487 ELSE IF ((MCAU.GE.10) .AND. (MCAU.LE.99)) THEN 3488C 3489 TENTH = DINT(AINT(DFLOAT(MCAU)/10.0D0)) 3490 UNITS = MCAU - TENTH*10 3491 NFIL = '_'//CHAR(TENTH+48)//CHAR(UNITS+48) 3492C 3493 ELSE IF ((MCAU.GE.100) .AND. (MCAU.LE.999)) THEN 3494C 3495 HUNDREDS = DINT(AINT(DFLOAT(MCAU)/100.0D0)) 3496 TENTH = DINT(AINT(DFLOAT(MCAU-HUNDREDS*100)/10.0D0)) 3497 UNITS = MCAU - HUNDREDS*100 - TENTH*10 3498 NFIL = CHAR(HUNDREDS+48)//CHAR(TENTH+48)//CHAR(UNITS+48) 3499C 3500 ELSE 3501 WRITE(LUPRI,*)'MCAU = ', MCAU 3502 CALL QUIT('Case not implemented in GET_FILENAME') 3503 END IF 3504C 3505 !Generate file name for C1-transformed integrals 3506 FNCKJDR = 'CC3_CAUINT_O_'//NFIL 3507 FNDKBCR = 'CC3_CAUINT_V_'//NFIL 3508C 3509 IF (LCDBG) THEN 3510 WRITE(LUPRI,*)'MCAU = ',MCAU 3511 WRITE(LUPRI,*)'FNCKJDR = ',FNCKJDR 3512 WRITE(LUPRI,*)'FNDKBCR = ',FNDKBCR 3513 WRITE(LUPRI,*)' ' 3514 END IF 3515C 3516 !End 3517 CALL QEXIT('GETFN') 3518C 3519 RETURN 3520 END 3521C /* Deck delete_file */ 3522 SUBROUTINE DELETE_FILES(FILENAME) 3523* 3524********************************************************************* 3525* Delete file(s) with name FILENAME 3526* 3527* NOTE: the file name FILENAME can contain wildcard 3528* characters, for example CALL DELETE_FILES('CC3_CAUINT_V*') 3529* 3530* Filip Pawlowski, Aarhus, 16-Apr-2004 3531********************************************************************* 3532* 3533 IMPLICIT NONE 3534C 3535 INTEGER ILEN,INFO 3536 CHARACTER*(*) FILENAME 3537 CHARACTER*(80) COMMAND 3538C 3539 !Determine length of the filename 3540 !(cannot be larger than 74 due to the way COMMAND is declared) 3541 ILEN = MIN(LEN(FILENAME),74) 3542C 3543 !Generate the command to delete the file 3544 WRITE(COMMAND,'(2A)') 'rm -f ', FILENAME(1:ILEN) 3545C 3546 !Execute the command 3547#if defined (SYS_CRAY) 3548 CALL INFO = ISHELL(COMMAND) 3549#else 3550 CALL SYSTEM(COMMAND) 3551#endif 3552 3553 RETURN 3554 END 3555 3556 3557