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_ft3b */ 20 SUBROUTINE CC3_FT3B(LISTL,IDLSTL,LISTB,IDLSTB,IOPTRES, 21 * OMEGA1,ISYRES,WORK,LWORK, 22 * FNCKJD,FNDKBC,FNDELD,FNTOC) 23 24C 25C Calculate <L2|[[H,E_AI],R_3]|HF> 26C 27C Initially construct 28C 29C t3(ai,Bj,Ck) = S^BC(aikj) + U^BC(aikj) + S^CB(aijk) + U^CB(aijk) 30C 31C W^BC(aikj) = sum_d X(ad)*t3(di,Bj,Ck) - sum_l X(li)*t3(al,Bj,Ck) 32C - sum_ld X(ld)*t2(ai,Cl)*t2(Bj,dk) 33C - sum_ld X(ld)*t2(ai,Bl)*t2(Ck,dj) 34C + ( -<mu3 | [[H,T1^Y],T2] | HF> -<mu3| [H,T2^Y] |HF> ) 35C 36C Written by Poul Jorgensen and Filip Pawlowski, Aarhus, Fall 2002 37C 38 IMPLICIT NONE 39C 40 INTEGER LWORK, IOPTRES 41 INTEGER ISINT1, ISINT2, ISYMT2, ISYMT3, ISYMW 42 INTEGER LU3SRTR, LUCKJDR, LUDELDR, LUDKBCR, LUDKBCR4 43 INTEGER KT1AM, KT2TP, KFOCKD, KEND0, LWRK0, IOPT 44 INTEGER IDLSTB, ISYMB, KT1B, KT2B 45 INTEGER ISINTR1, ISINTR2, KLAMDP, KLAMDH, KEND1, LWRK1 46 INTEGER KTROC, KTROC1, KW3XOG1, KT3OG2, KW3XOGX1 47 INTEGER KINTOC, KEND2, LWRK2 48 INTEGER IOFF, ISYMD, ISCKB1, ISCKB2, ISCKB2Y 49 INTEGER KTRVI, KTRVI1, KT3VDG2, KT3VDG3, KW3XVDG1, KW3XVDGX1 50 INTEGER KEND3, LWRK3, KEND4, LWRK4, KINTVI, LENSQ, LENSQW 51 INTEGER ISYALJ, ISYALJ2, ISYMBD, ISCKIJ, ISCKD2, ISCKD2Y 52 INTEGER ISWMAT, KSMAT, KSMAT3, KUMAT, KUMAT3, KDIAG, KDIAGW 53 INTEGER KWMAT, KINDSQ, KINDSQW, KINDEX, KINDEX2, KTMAT 54 INTEGER KW3XVBG1, KW3XVDGX2, KT3VBG2, KT3VBG3, KEND5, LWRK5 55 INTEGER LUCKJD, LUTOC, LUDKBC, LUDELD 56 INTEGER KSAVE, KEND00, LWRK00, KEND01, LWRK01 57 INTEGER ISYMRB, ISYRES 58 INTEGER KLAMPB, KLAMHB, KFCKB, LUFCK, IVEC, IDLSTL, ISYML 59 INTEGER ISYFCKB 60 INTEGER IRREP, IERR, KXIAJB, KYIAJB, IOPTTCME, ISYMK, KOFF1, KOFF2 61 INTEGER IRREP2, ISYCKB 62 INTEGER KXAIJB, KYAIJB, KYAIBJ 63 INTEGER ISINT3, KL2TP, ISYML2 64C 65C Functions : 66C 67 INTEGER ILSTSYM 68C 69 integer kx3am 70C 71#if defined (SYS_CRAY) 72 REAL WORK(LWORK) 73 REAL OMEGA1(*) 74 REAL FREQB, TCON 75 REAL XNORMVAL, DDOT, HALF, TWO, XT2TP, XUMAT 76#else 77 DOUBLE PRECISION WORK(LWORK) 78 DOUBLE PRECISION OMEGA1(*) 79 DOUBLE PRECISION FREQB, TCON 80 DOUBLE PRECISION XNORMVAL, DDOT, HALF, TWO, XT2TP, XUMAT 81#endif 82 CHARACTER*(*) FNCKJD, FNDKBC, FNDELD, FNTOC 83 CHARACTER*12 FN3SRTR, FNCKJDR, FNDELDR, FNDKBCR, FNDKBCR4 84 CHARACTER*10 MODEL 85 CHARACTER*8 LABELB 86 CHARACTER*3 LISTL, LISTB 87 CHARACTER CDUMMY*1 88C 89C 90#include "priunit.h" 91#include "dummy.h" 92#include "iratdef.h" 93#include "ccsdsym.h" 94#include "inftap.h" 95#include "ccinftap.h" 96#include "ccorb.h" 97#include "ccsdinp.h" 98#include "ccr1rsp.h" 99#include "ccexci.h" 100C 101 PARAMETER(HALF = 0.5D0, TWO = 2.0D0) 102 PARAMETER(FN3SRTR = 'CCSDT_FBMAT1', FNCKJDR = 'CCSDT_FBMAT2', 103 * FNDELDR = 'CCSDT_FBMAT3', FNDKBCR = 'CCSDT_FBMAT4', 104 * FNDKBCR4 = 'CCSDT_FBMAT5') 105C 106 CALL QENTER('CC3_FT3B') 107C 108C-------------------------- 109C Save and set flags. 110C-------------------------- 111C 112C Set symmetry flags. 113C 114C 115C ISYMT2 is symmetry of T2TP 116C ISINT2 is symmetry of integrals in triples equation (ISINT2=1) 117C ISINT1 is symmetry of integrals in contraction (ISINT1=1) 118C ISYMT3 is symmetry of S and U intermediate 119C ISYMW is symmetry of W intermediate 120C ISINT3 is symmetry of integrals used in contraction with W 121C ISYRES is symmetry of result vector 122C 123C------------------------------------------------------------- 124C 125COMMENT COMMENT COMMENT 126c 127c summing up contributions in omega1eff 128c must be changed when introduced in real implementation 129c 130c CALL DZERO(OMEGA1EFF,NT1AM(1)) 131c 132COMMENT COMMENT COMMENT 133 IPRCC = IPRINT 134 ISINT1 = 1 135 ISINT2 = 1 136 ISINT3 = 1 137 ISYMT2 = 1 138 ISYMT3 = MULD2H(ISINT2,ISYMT2) 139C 140C-------------------------------- 141C Open files 142C-------------------------------- 143C 144 LUCKJD = -1 145 LUTOC = -1 146 LUDKBC = -1 147 LUDELD = -1 148C 149 CALL WOPEN2(LUCKJD,FNCKJD,64,0) 150 CALL WOPEN2(LUTOC,FNTOC,64,0) 151 CALL WOPEN2(LUDKBC,FNDKBC,64,0) 152 CALL WOPEN2(LUDELD,FNDELD,64,0) 153C 154C-------------------------------- 155C Open temporary files 156C-------------------------------- 157C 158 LU3SRTR = -1 159 LUCKJDR = -1 160 LUDELDR = -1 161 LUDKBCR = -1 162 LUDKBCR4 = -1 163C 164 CALL WOPEN2(LU3SRTR,FN3SRTR,64,0) 165 CALL WOPEN2(LUCKJDR,FNCKJDR,64,0) 166 CALL WOPEN2(LUDELDR,FNDELDR,64,0) 167 CALL WOPEN2(LUDKBCR,FNDKBCR,64,0) 168 CALL WOPEN2(LUDKBCR4,FNDKBCR4,64,0) 169C 170C---------------------------------------------- 171C Calculate the zeroth order stuff once 172C---------------------------------------------- 173C 174 KT2TP = 1 175 KL2TP = KT2TP + NT2SQ(ISYMT2) 176 KFOCKD = KL2TP + NT2SQ(ISYMT2) 177 KLAMDP = KFOCKD + NORBTS 178 KLAMDH = KLAMDP + NLAMDT 179 KXIAJB = KLAMDH + NLAMDT 180 KYIAJB = KXIAJB + NT2AM(ISINT1) 181 KXAIJB = KYIAJB + NT2AM(ISINT1) 182 KYAIJB = KXAIJB + NT2SQ(ISINT1) 183 KYAIBJ = KYAIJB + NT2SQ(ISINT1) 184 KEND00 = KYAIBJ + NT2SQ(ISINT1) 185 LWRK00 = LWORK - KEND00 186C 187 KT1AM = KEND00 188 KEND01 = KT1AM + NT1AM(ISYMT2) 189 LWRK01 = LWORK - KEND01 190C 191 IF (LWRK01 .LT. 0) THEN 192 CALL QUIT('Out of memory in CC3_FT3B (zeroth allo.') 193 ENDIF 194C 195C----------------------- 196C read in doubles component of L0 197C----------------------- 198C 199 ISYML2 = ILSTSYM(LISTL,IDLSTL) 200 IF (LWRK01 .LT. NT2SQ(ISYML2)) THEN 201 CALL QUIT('Not enough memory to construct L2TP in CC3_FT3B') 202 ENDIF 203 204 IOPT = 2 205 CALL GET_T1_T2(IOPT,.FALSE.,DUMMY,WORK(KL2TP),LISTL,IDLSTL,ISYML2, 206 * WORK(KEND01),LWRK01) 207C 208C------------------------ 209C Construct L(ia,jb). 210C------------------------ 211C 212 REWIND(LUIAJB) 213 CALL READI(LUIAJB,IRAT*NT2AM(ISINT1),WORK(KXIAJB)) 214 CALL DCOPY(NT2AM(ISINT1),WORK(KXIAJB),1,WORK(KYIAJB),1) 215C 216 IOPTTCME = 1 217 CALL CCSD_TCMEPK(WORK(KYIAJB),1.0D0,ISINT1,IOPTTCME) 218C 219 IF ( IPRINT .GT. 55) THEN 220 XNORMVAL = DDOT(NT2AM(ISINT1),WORK(KYIAJB),1, 221 * WORK(KYIAJB),1) 222 WRITE(LUPRI,*) 'Norm of L(IAJB) ',XNORMVAL 223 ENDIF 224 CALL CC_T2SQ(WORK(KXIAJB),WORK(KEND01),ISINT1) 225 CALL CC_T2SQ(WORK(KYIAJB),WORK(KYAIBJ),ISINT1) 226C 227 CALL CC3_T2TP(WORK(KXAIJB),WORK(KEND01),ISINT1) 228 CALL CC3_T2TP(WORK(KYAIJB),WORK(KYAIBJ),ISINT1) 229C 230 IF (IPRINT .GT. 55) THEN 231 XT2TP = DDOT(NT2SQ(ISINT1),WORK(KYAIJB),1,WORK(KYAIJB),1) 232 WRITE(LUPRI,*) 'Norm of L(AIJB) ', XT2TP 233 ENDIF 234C 235C---------------------------------------------------------------- 236C Read t1 and calculate the zero'th order Lambda matrices 237C---------------------------------------------------------------- 238C 239 CALL GET_LAMBDA0(WORK(KLAMDP),WORK(KLAMDH),WORK(KEND01),LWRK01) 240C 241C------------------------------------------- 242C Read in t2 , square it and reorder 243C------------------------------------------- 244C 245 IOPT = 2 246 CALL GET_T1_T2(IOPT,.FALSE.,DUMMY,WORK(KT2TP),'R0',0,ISYMT2, 247 * WORK(KEND01),LWRK01) 248 249 IF (IPRINT .GT. 55) THEN 250 XNORMVAL = DDOT(NT2SQ(ISYMT2),WORK(KT2TP),1,WORK(KT2TP),1) 251 WRITE(LUPRI,*) 'Norm of T2TP ',XNORMVAL 252 ENDIF 253C 254C-------------------------------------- 255C Read in orbital energies 256C-------------------------------------- 257C 258C 259 CALL GET_ORBEN(WORK(KFOCKD),WORK(KEND01),LWRK01) 260C If we want to sum the T3 amplitudes 261 if (.false.) then 262 kx3am = kend00 263 kend00 = kx3am + nt1amx*nt1amx*nt1amx 264 call dzero(work(kx3am),nt1amx*nt1amx*nt1amx) 265 lwrk00 = lwork - kend00 266 if (lwrk00 .lt. 0) then 267 write(lupri,*) 'Memory available : ',lwork 268 write(lupri,*) 'Memory needed : ',kend00 269 call quit('Insufficient space (T3) in CC3_FT3B') 270 END IF 271 endif 272C 273C determining R1 labels 274C 275 IF (LISTB(1:3).EQ.'R1 ') THEN 276 ISYMRB = ISYLRT(IDLSTB) 277 FREQB = FRQLRT(IDLSTB) 278 LABELB = LRTLBL(IDLSTB) 279 ELSE IF (LISTB(1:3).EQ.'RE ') THEN 280 ISYMRB = ILSTSYM(LISTB,IDLSTB) 281 FREQB = EIGVAL(IDLSTB) 282 LABELB = '- none -' 283C 284 ELSE 285 CALL QUIT('Unknown list in CC3_FT3B.') 286 END IF 287C 288 !symmetry check 289c ISYMW = MULD2H(ISYMT3,ISYMRB) 290c IF (ISYRES.NE.MULD2H(ISYMW,ISINT3)) THEN 291c WRITE(LUPRI,*) 'INCONSISTENT SYMMETRY SPECIFICATION' 292c WRITE(LUPRI,*) 293c * 'ISYRES =',ISYRES,'ISYMW =',ISYMW,'ISINT3 =',ISINT3 294c STOP' CC3_FT3B inconsistent symmetry specification' 295c END IF 296C 297 ISYFCKB = MULD2H(ISYMOP,ISYMRB) 298C 299C------------------------------------------------- 300C Read T1^B and T2^B 301C Calculate (ck|de)-tilde(B) and (ck|lm)-tilde(B) 302C Used to construct T3^B 303C------------------------------------------------- 304C 305 KT2B = KEND00 306 KFCKB = KT2B + NT2SQ(ISYMRB) 307 KEND0 = KFCKB + N2BST(ISYFCKB) 308 LWRK0 = LWORK - KEND0 309C 310 KT1B = KEND0 311 KLAMPB = KT1B + NT1AM(ISYMRB) 312 KLAMHB = KLAMPB + NLAMDT 313 KEND1 = KLAMHB + NLAMDT 314 LWRK1 = LWORK - KEND1 315C 316 IF (LWRK1 .LT. NT2SQ(ISYMRB)) THEN 317 CALL QUIT('Out of memory in CC3_XI (TOT_T3Y) ') 318 ENDIF 319C 320C Readin 321C 322 IOPT = 3 323 CALL GET_T1_T2(IOPT,.TRUE.,WORK(KT1B),WORK(KT2B),LISTB,IDLSTB, 324 * ISYMRB,WORK(KEND1),LWRK1) 325C 326 IF (IPRINT .GT. 55) THEN 327 XNORMVAL = DDOT(NT2SQ(ISYMRB),WORK(KT2B),1,WORK(KT2B),1) 328 WRITE(LUPRI,*) 'Norm of T2B ',XNORMVAL 329 XNORMVAL = DDOT(NT1AM(ISYMRB),WORK(KT1B),1,WORK(KT1B),1) 330 WRITE(LUPRI,*) 'Norm of T1B ',XNORMVAL 331 ENDIF 332C 333C Integrals 334 335 ISINTR1 = MULD2H(ISINT1,ISYMRB) 336 ISINTR2 = MULD2H(ISINT2,ISYMRB) 337C 338 CALL CC3_BARINT(WORK(KT1B),ISYMRB,WORK(KLAMDP), 339 * WORK(KLAMDH),WORK(KEND1),LWRK1, 340 * LU3SRTR,FN3SRTR,LUCKJDR,FNCKJDR) 341C 342 CALL CC3_SORT1(WORK(KEND1),LWRK1,2,ISINTR1,LU3SRTR,FN3SRTR, 343 * LUDELDR,FNDELDR,IDUMMY,CDUMMY,IDUMMY,CDUMMY, 344 * IDUMMY,CDUMMY) 345C 346 CALL CC3_SINT(WORK(KLAMDH),WORK(KEND1),LWRK1,ISINTR1, 347 * LUDELDR,FNDELDR,LUDKBCR,FNDKBCR) 348 349 CALL CC3_TCME(WORK(KLAMDP),ISINTR1,WORK(KEND1),LWRK1, 350 * IDUMMY,CDUMMY,LUDKBCR,FNDKBCR, 351 * IDUMMY,CDUMMY,IDUMMY,CDUMMY, 352 * IDUMMY,CDUMMY,LUDKBCR4,FNDKBCR4,2) 353 354C--------------------------------------- 355C Calculate the Lambda^B 356C--------------------------------------- 357C 358 CALL CCLR_LAMTRA(WORK(KLAMDP),WORK(KLAMPB),WORK(KLAMDH), 359 * WORK(KLAMHB),WORK(KT1B),ISYMRB) 360C 361C------------------------------------------ 362C Calculate the F^B matrix 363C------------------------------------------ 364C 365 IF ( .NOT.(LISTB(1:3).EQ.'RE ') ) THEN 366 CALL GET_FOCKX(WORK(KFCKB),LABELB,IDLSTB,ISYMRB,WORK(KLAMDP),1, 367 * WORK(KLAMDH),1,WORK(KEND1),LWRK1) 368 END IF 369C 370C----------------------------- 371C Read occupied integrals. 372C----------------------------- 373C 374C Memory allocation. 375C 376C--------------------------------------------------------- 377C Work allocation 378C--------------------------------------------------------- 379C 380 KW3XOG1 = KEND0 381 KT3OG2= KW3XOG1 + NTRAOC(ISINT2) 382 KEND1 = KT3OG2+ NTRAOC(ISINT2) 383 LWRK1 = LWORK - KEND1 384C 385 KW3XOGX1 = KEND1 386 KEND1 = KW3XOGX1 + NTRAOC(ISINTR2) 387C 388 KINTOC = KEND1 389 KEND2 = KINTOC 390 * + MAX(NTOTOC(ISINT2),NTOTOC(ISINT1),NTOTOC(ISYMOP)) 391 LWRK2 = LWORK - KEND2 392C 393 IF (LWRK2 .LT. 0) THEN 394 WRITE(LUPRI,*) 'Memory available : ',LWORK 395 WRITE(LUPRI,*) 'Memory needed : ',KEND2 396 CALL QUIT('Insufficient space in CC3_XI') 397 END IF 398C 399C Occupied integrals. 400C 401C----------------------- 402C Read in integrals for SMAT etc. 403C----------------------- 404C 405 CALL INTOCC_T30(LUCKJD,FNCKJD,WORK(KLAMDP),ISINT2,WORK(KW3XOG1), 406 * WORK(KT3OG2),WORK(KEND2),LWRK2) 407C 408C------------------------------------------ 409C B transformed Occupied integrals. 410C----------------------------------------- 411C 412 CALL INTOCC_T3X(LUCKJDR,FNCKJDR,WORK(KLAMDP),ISINTR2, 413 * WORK(KW3XOGX1),WORK(KEND2),LWRK2) 414C 415C---------------------------- 416C General loop structure. 417C---------------------------- 418C 419 DO ISYMD = 1,NSYM 420 421 ISYCKB = MULD2H(ISYMOP,ISYMD) 422 ISCKB2 = MULD2H(ISINT2,ISYMD) 423 ISCKB2Y = MULD2H(ISINTR2,ISYMD) 424 ISCKB1 = MULD2H(ISINT1,ISYMD) 425C 426 IF (IPRINT .GT. 55) THEN 427C 428 WRITE(LUPRI,*) 'In CC3_FT3B: ISYCKB :',ISYCKB 429 WRITE(LUPRI,*) 'In CC3_FT3B: ISCKB2 :',ISCKB2 430 WRITE(LUPRI,*) 'In CC3_FT3B: ISCKB2Y:',ISCKB2Y 431 WRITE(LUPRI,*) 'In CC3_FT3B: ISCKB1 :',ISCKB1 432 433 ENDIF 434 435C 436C-------------------------- 437C Memory allocation. 438C-------------------------- 439C 440 KT3VDG3 = KEND1 441 KT3VDG2 = KT3VDG3 + NCKATR(ISCKB2) 442 KEND2 = KT3VDG2 + NCKATR(ISCKB2) 443 LWRK2 = LWORK - KEND2 444 445 KW3XVDG1 = KEND2 446 KEND3 = KW3XVDG1 + NCKATR(ISCKB2) 447 LWRK3 = LWORK - KEND3 448C 449 KW3XVDGX1 = KEND3 450 KEND3 = KW3XVDGX1 + NCKATR(ISCKB2Y) 451 LWRK3 = LWORK - KEND3 452 453 KINTVI = KEND3 454 KEND4 = KINTVI 455 * + MAX(NCKA(ISCKB2),NCKA(ISCKB1),NCKA(ISYCKB)) 456 LWRK4 = LWORK - KEND4 457C 458 IF (LWRK4 .LT. 0) THEN 459 WRITE(LUPRI,*) 'Memory available : ',LWORK 460 WRITE(LUPRI,*) 'Memory needed : ',KEND4 461 CALL QUIT('Insufficient space in CC3_XI') 462 END IF 463C 464C--------------------- 465C Sum over D 466C--------------------- 467C 468 DO D = 1,NVIR(ISYMD) 469C 470C----------------------------------------------- 471C Read virtual integrals used in first s3am. 472C----------------------------------------------- 473C 474 CALL INTVIR_T30_D(LUDKBC,FNDKBC,LUDELD,FNDELD,ISINT2, 475 * WORK(KW3XVDG1),WORK(KT3VDG2), 476 * WORK(KT3VDG3),WORK(KLAMDH),ISYMD,D, 477 * WORK(KEND4),LWRK4) 478C 479C----------------------------------------------------------------------- 480C Read B transformed virtual integrals used for W for TOT_T3Y 481C----------------------------------------------------------------------- 482C 483 IOFF = ICKBD(ISCKB2Y,ISYMD) + NCKATR(ISCKB2Y)*(D - 1) + 1 484 IF (NCKATR(ISCKB2Y) .GT. 0) THEN 485 CALL GETWA2(LUDKBCR,FNDKBCR,WORK(KW3XVDGX1),IOFF, 486 & NCKATR(ISCKB2Y)) 487 ENDIF 488C 489C-------------------------------------------------------- 490C Sum over ISYMB 491C-------------------------------------------------------- 492C 493 DO ISYMB = 1,NSYM 494 495 ISYALJ = MULD2H(ISYMB,ISYMT2) 496 ISYALJ2 = MULD2H(ISYMD,ISYMT2) 497 ISYMBD = MULD2H(ISYMB,ISYMD) 498 ISCKIJ = MULD2H(ISYMBD,ISYMT3) 499 ISCKD2 = MULD2H(ISINT2,ISYMB) 500 ISCKD2Y = MULD2H(ISINTR2,ISYMB) 501 ISWMAT = MULD2H(ISYMRB,ISCKIJ) 502 503 IF (IPRINT .GT. 55) THEN 504 505 WRITE(LUPRI,*) 'In CC3_FT3B: ISYMD :',ISYMD 506 WRITE(LUPRI,*) 'In CC3_FT3B: ISYMB :',ISYMB 507 WRITE(LUPRI,*) 'In CC3_FT3B: ISYALJ :',ISYALJ 508 WRITE(LUPRI,*) 'In CC3_FT3B: ISYALJ2:',ISYALJ2 509 WRITE(LUPRI,*) 'In CC3_FT3B: ISYMBD :',ISYMBD 510 WRITE(LUPRI,*) 'In CC3_FT3B: ISCKIJ :',ISCKIJ 511 WRITE(LUPRI,*) 'In CC3_FT3B: ISWMAT :',ISWMAT 512 513 ENDIF 514 515C Can use kend3 since we do not need the integrals anymore. 516 KSMAT = KEND3 517 KSMAT3 = KSMAT + NCKIJ(ISCKIJ) 518 KUMAT = KSMAT3 + NCKIJ(ISCKIJ) 519 KUMAT3 = KUMAT + NCKIJ(ISCKIJ) 520 KDIAG = KUMAT3 + NCKIJ(ISCKIJ) 521 KDIAGW = KDIAG + NCKIJ(ISCKIJ) 522 KWMAT = KDIAGW + NCKIJ(ISWMAT) 523 KINDSQW = KWMAT + NCKIJ(ISWMAT) 524 KINDSQ = KINDSQW + (6*NCKIJ(ISWMAT) - 1)/IRAT + 1 525 KINDEX = KINDSQ + (6*NCKIJ(ISCKIJ) - 1)/IRAT + 1 526 KINDEX2 = KINDEX + (NCKI(ISYALJ) - 1)/IRAT + 1 527 KTMAT = KINDEX2 + (NCKI(ISYALJ2) - 1)/IRAT + 1 528 KW3XVBG1 = KTMAT + MAX(NCKIJ(ISCKIJ),NCKIJ(ISWMAT)) 529 KT3VBG2 = KW3XVBG1 + NCKATR(ISCKD2) 530 KT3VBG3 = KT3VBG2 + NCKATR(ISCKD2) 531 KEND4 = KT3VBG3 + NCKATR(ISCKD2) 532 LWRK4 = LWORK - KEND4 533C 534 KW3XVDGX2 = KEND4 535 KEND4 = KW3XVDGX2 + NCKATR(ISCKD2Y) 536 537C 538 KINTVI = KEND4 539 KEND5 = KINTVI + NCKA(ISCKD2) 540 LWRK5 = LWORK - KEND5 541C 542 IF (LWRK5 .LT. 0) THEN 543 WRITE(LUPRI,*) 'Memory available : ',LWORK 544 WRITE(LUPRI,*) 'Memory needed : ',KEND5 545 CALL QUIT('Insufficient space in CC3_XI') 546 END IF 547C 548C--------------------------------------------- 549C Construct part of the diagonal. 550C--------------------------------------------- 551C 552 CALL CC3_DIAG(WORK(KDIAG),WORK(KFOCKD),ISCKIJ) 553 CALL CC3_DIAG(WORK(KDIAGW),WORK(KFOCKD),ISWMAT) 554C 555 IF (IPRINT .GT. 55) THEN 556 XNORMVAL = DDOT(NCKIJ(ISCKIJ),WORK(KDIAG),1, 557 * WORK(KDIAG),1) 558 WRITE(LUPRI,*) 'Norm of DIA ',XNORMVAL 559 ENDIF 560C 561C------------------------------------- 562C Construct index arrays. 563C---------------------------------- 564C 565 LENSQ = NCKIJ(ISCKIJ) 566 CALL CC3_INDSQ(WORK(KINDSQ),LENSQ,ISCKIJ) 567 CALL CC3_INDEX(WORK(KINDEX),ISYALJ) 568 CALL CC3_INDEX(WORK(KINDEX2),ISYALJ2) 569 LENSQW = NCKIJ(ISWMAT) 570 CALL CC3_INDSQ(WORK(KINDSQW),LENSQW,ISWMAT) 571 572 DO B = 1,NVIR(ISYMB) 573C 574C------------------------------------ 575C Initialise 576C--------------------------------------- 577C 578 CALL DZERO(WORK(KWMAT),NCKIJ(ISWMAT)) 579C 580C----------------------------------------------- 581C Read virtual integrals used in second s3am. 582C----------------------------------------------- 583C 584 CALL INTVIR_T30_D(LUDKBC,FNDKBC,LUDELD,FNDELD,ISINT2, 585 * WORK(KW3XVBG1),WORK(KT3VBG2), 586 * WORK(KT3VBG3),WORK(KLAMDH),ISYMB,B, 587 * WORK(KEND5),LWRK5) 588C 589C-------------------------------------------------------------------- 590C Read B transformed virtual integrals used in W 591C-------------------------------------------------------------------- 592C 593 IOFF = ICKBD(ISCKD2Y,ISYMB) + 594 & NCKATR(ISCKD2Y)*(B - 1) + 1 595 IF (NCKATR(ISCKD2Y) .GT. 0) THEN 596 CALL GETWA2(LUDKBCR,FNDKBCR,WORK(KW3XVDGX2),IOFF, 597 & NCKATR(ISCKD2Y)) 598 ENDIF 599C 600C------------------------------------------------------------------------ 601C Calculate the S(ci,bk,dj) matrix for T3 for B,D. 602C------------------------------------------------------------------- 603C 604 CALL GET_T30_BD(ISYMT3,ISINT2,WORK(KT2TP),ISYMT2, 605 * WORK(KTMAT),WORK(KFOCKD),WORK(KDIAG), 606 * WORK(KINDSQ),LENSQ,WORK(KSMAT), 607 * WORK(KW3XVDG1),WORK(KT3VDG2), 608 * WORK(KW3XOG1),WORK(KINDEX), 609 * WORK(KSMAT3),WORK(KW3XVBG1), 610 * WORK(KT3VBG2),WORK(KINDEX2), 611 * WORK(KUMAT),WORK(KT3VDG3), 612 * WORK(KT3OG2),WORK(KUMAT3), 613 * WORK(KT3VBG3),ISYMB,B,ISYMD,D, 614 * ISCKIJ,WORK(KEND4),LWRK4) 615C 616C------------------------------------------------------ 617C Based on T3 BD amplitudes calculate W BD intermediates 618C------------------------------------------------------ 619C 620 621 IF (LISTB(1:3).EQ.'R1 ') THEN 622 623 CALL GET_T3X_BD(ISYMRB,WORK(KWMAT),ISWMAT, 624 * WORK(KTMAT),ISCKIJ,WORK(KFCKB), 625 * ISYMRB,WORK(KT2TP),ISYMT2, 626 * WORK(KT2B),ISYMRB,WORK(KW3XVDG1), 627 * WORK(KW3XVBG1),WORK(KW3XOG1), 628 * ISINT2,WORK(KW3XVDGX1), 629 * WORK(KW3XVDGX2),WORK(KW3XOGX1), 630 * ISINTR2,FREQB,WORK(KDIAGW), 631 * WORK(KFOCKD),WORK(KINDSQW),LENSQW, 632 * B,ISYMB,D,ISYMD,WORK(KEND4),LWRK4) 633 634 END IF 635C 636 IF (LISTB(1:3).EQ.'RE ') THEN 637 CALL WBD_GROUND(WORK(KT2B),ISYMRB,WORK(KTMAT), 638 * WORK(KW3XVDG1),WORK(KW3XVBG1), 639 * WORK(KW3XOG1),ISINT2,WORK(KWMAT), 640 * WORK(KEND4),LWRK4, 641 * WORK(KINDSQW),LENSQW, 642 * ISYMB,B,ISYMD,D) 643C 644 CALL WBD_GROUND(WORK(KT2TP),ISYMT2,WORK(KTMAT), 645 * WORK(KW3XVDGX1),WORK(KW3XVDGX2), 646 * WORK(KW3XOGX1),ISINTR2,WORK(KWMAT), 647 * WORK(KEND4),LWRK4, 648 * WORK(KINDSQW),LENSQW, 649 * ISYMB,B,ISYMD,D) 650 651C 652C------------------------------------------------ 653C Divide by the energy difference and 654C remove the forbidden elements 655C------------------------------------------------ 656C 657 CALL WBD_DIA(B,ISYMB,D,ISYMD,FREQB, 658 * ISWMAT,WORK(KWMAT), 659 * WORK(KDIAGW),WORK(KFOCKD)) 660 661 CALL T3_FORBIDDEN(WORK(KWMAT),ISYMRB, 662 * ISYMB,B,ISYMD,D) 663 END IF 664C 665C---------------------------------------------------------------------- 666C Calculate the term <L2|[H,E_ai],W^BD(3)]]|HF> 667C---------------------------------------------------------------------- 668C 669* call sum_pt3(work(kwmat),isymb,b,isymd,d, 670* * iswmat,work(kx3am),4) 671C 672 CALL T3_ONEL1W(OMEGA1,WORK(KWMAT),WORK(KTMAT), 673 * ISWMAT,WORK(KXIAJB),WORK(KXAIJB), 674 * WORK(KYIAJB),WORK(KYAIJB),ISINT3, 675 * WORK(KL2TP),ISYML2,WORK(KINDSQW), 676 * LENSQW,WORK(KEND4), 677 * LWRK4,ISYMB,B,ISYMD,D) 678C 679 CALL T3_ONEL2W(OMEGA1,WORK(KWMAT),WORK(KTMAT), 680 * ISWMAT,WORK(KXIAJB),WORK(KYIAJB), 681 * ISINT3, 682 * WORK(KL2TP),ISYML2,WORK(KINDSQW), 683 * LENSQW,WORK(KEND4), 684 * LWRK4,ISYMB,B,ISYMD,D) 685C 686 CALL T3_ONEL3W(OMEGA1,WORK(KWMAT),WORK(KTMAT), 687 * ISWMAT,WORK(KYIAJB),WORK(KYAIBJ), 688 * ISINT3, 689 * WORK(KL2TP),ISYML2,WORK(KINDSQW), 690 * LENSQW,WORK(KEND4), 691 * LWRK4,ISYMB,B,ISYMD,D) 692C 693 END DO ! B 694 END DO ! ISYMB 695 END DO ! D 696 END DO ! ISYMD 697C 698* write(lupri,*)'omega1 ' 699* do i = 1,nt1am(isyres) 700* write(lupri,*) i , OMEGA1(i) 701* end do 702* 703* write(lupri,*) 'cont to t3x in CC3_FT3B' 704* call print_pt3(work(kx3am),1,4) 705* write(lupri,*) 'The summed S terms : ' 706* call print_pt3(work(kx3am),1,1) 707C 708C-------------------------------- 709C Close files 710C-------------------------------- 711C 712 CALL WCLOSE2(LUCKJD,FNCKJD,'KEEP') 713 CALL WCLOSE2(LUTOC,FNTOC,'KEEP') 714 CALL WCLOSE2(LUDKBC,FNDKBC,'KEEP') 715 CALL WCLOSE2(LUDELD,FNDELD,'KEEP') 716C 717C-------------------------------- 718C Close files for "response" 719C-------------------------------- 720C 721 CALL WCLOSE2(LU3SRTR,FN3SRTR,'DELETE') 722 CALL WCLOSE2(LUCKJDR,FNCKJDR,'DELETE') 723 CALL WCLOSE2(LUDELDR,FNDELDR,'DELETE') 724 CALL WCLOSE2(LUDKBCR,FNDKBCR,'DELETE') 725 CALL WCLOSE2(LUDKBCR4,FNDKBCR4,'DELETE') 726C 727C------------- 728C End 729C------------- 730C 731 CALL QEXIT('CC3_FT3B') 732C 733 RETURN 734 END 735C /* Deck t3_onel1w */ 736 SUBROUTINE T3_ONEL1W(OMEGA1,WMAT,TMAT,ISYMIM,XIAJB,XAIJB,YIAJB, 737 * YAIJB,ISYINT,C2TP,ISYMC2,INDSQ,LENSQ,WORK, 738 * LWORK,ISYMB,B,ISYMD,D) 739C 740C Written by Filip Pawlowski and Poul Jorgensen based on 741C SUBROUTINE T3_ONEL1 by K. Hald, Spring 2002. 742C 743C Calculate the term t^{def}_{lmn} T^{fd}_{mi} g_{nela} 744C - t^{def}_{lnm} T^{fd}_{mi} L_{nela} 745C using W intermediates 746C 747C g contributions 748C 749C 750C 1) 751C W^{ed}(flnm) T^{df}_{mi} g_{nela} 752C 2) 753C W^{ed}(fnlm) T^{fd}_{mi} g_{nela} 754C 3) 755C W^{fd}(emln) T^{fd}_{mi} g_{nela} 756C 757C L contributions 758C 759C 4) 760C W^{ed}(flmn) T^{df}_{mi} L_{nela} 761C 5) 762C W^{ed}(fmln) T^{fd}_{mi} L_{nela} 763C 6) 764C W^{fd}(enlm) T^{fd}_{mi} L_{nela} 765C 766C 767C XIAJB contains g and YIAJB contains L 768C 769 IMPLICIT NONE 770C 771#include "priunit.h" 772#include "ccsdsym.h" 773#include "ccorb.h" 774C 775 INTEGER ISYMIM, ISYINT, ISYMC2, LENSQ, LWORK, ISYMB, ISYMD 776 INTEGER INDSQ(LENSQ,6), INDEX 777 INTEGER ISYRE1, ISYRES, ISYMBD, ISFLMN, ISYANL, LENGTH 778 INTEGER ISYFIM, KTMAT, KC2TEMP, KINT, KEND1, LWRK1 779 INTEGER ISYMM, ISYMFI, ISYMF, KOFF1, KOFF2 780 INTEGER ISYML, ISYMAN, ISYMA, ISYMN, ISYMLN, NBN, NAN, NAL 781 INTEGER ISYMI, ISYMAB, ISYMFM, KOFF3, NUMBFM, NUMBLN, NUMBA 782 INTEGER ISYMBN, ISYMAL, KINT2 783 INTEGER ISYMMI, KMIMAT, KEND2, LWRK2, ISYBMI, ISYMBM 784 INTEGER ISENLM, ISYENL, NUMENL, NUMM, NUMA 785 INTEGER ISYRE2 786C 787#if defined (SYS_CRAY) 788 REAL OMEGA1(*), WMAT(*), TMAT(*), XIAJB(*), XAIJB(*), YAIJB(*) 789 REAL YIAJB(*), C2TP(*), WORK(LWORK), ZERO, ONE 790#else 791 DOUBLE PRECISION OMEGA1(*), WMAT(*), TMAT(*), XIAJB(*), XAIJB(*) 792 DOUBLE PRECISION YAIJB(*) 793 DOUBLE PRECISION YIAJB(*), C2TP(*), WORK(LWORK), ZERO, ONE 794#endif 795C 796 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 797C 798 INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 799C 800 CALL QENTER('T3_ONEL1W') 801C 802 ISYMBD = MULD2H(ISYMB,ISYMD) 803C 804 ISYRE1 = MULD2H(ISYMIM,ISYMC2) 805 ISYRE2 = MULD2H(ISYRE1,ISYMBD) 806 ISYRES = MULD2H(ISYRE2,ISYINT) 807C 808 ISFLMN = ISYMIM 809 ISYANL = MULD2H(ISYMB,ISYINT) 810C 811 LENGTH = NCKIJ(ISFLMN) 812C 813C calculate the terms 814C 815C 1) 816C W^{ed}(flnm) T^{df}_{mi} g_{nela} 817C 2) 818C W^{ed}(fnlm) T^{fd}_{mi} g_{nela} 819C 820C----------------------------- 821C Sort C2 822C----------------------------- 823C 824 ISYFIM = MULD2H(ISYMC2,ISYMD) 825C 826 KTMAT = 1 827 KC2TEMP = KTMAT + NCKIJ(ISFLMN) 828 KINT = KC2TEMP + NMAIJA(ISYFIM) 829 KINT2 = KINT + NCKI(ISYANL) 830 KEND1 = KINT2 + NCKI(ISYANL) 831 LWRK1 = LWORK - KEND1 832C 833 IF (LWRK1 .LT. 0) THEN 834 CALL QUIT('Out of memory in T3_ONEL1W (sort)') 835 ENDIF 836C 837 DO ISYMM = 1, NSYM 838 ISYMFI = MULD2H(ISYFIM,ISYMM) 839 DO ISYMF = 1, NSYM 840 ISYMI = MULD2H(ISYMFI,ISYMF) 841 ISYMFM = MULD2H(ISYMF,ISYMM) 842C 843 DO M = 1, NRHF(ISYMM) 844 DO I = 1, NRHF(ISYMI) 845C 846 KOFF1 = IT2SP(ISYFIM,ISYMD) 847 * + NCKI(ISYFIM)*(D-1) 848 * + ICKI(ISYMFI,ISYMM) 849 * + NT1AM(ISYMFI)*(M-1) 850 * + IT1AM(ISYMF,ISYMI) 851 * + NVIR(ISYMF)*(I-1) 852 * + 1 853C 854 KOFF2 = KC2TEMP 855 * + ICKI(ISYMFM,ISYMI) 856 * + NT1AM(ISYMFM)*(I-1) 857 * + IT1AM(ISYMF,ISYMM) 858 * + NVIR(ISYMF)*(M-1) 859 860C 861 CALL DCOPY(NVIR(ISYMF),C2TP(KOFF1),1,WORK(KOFF2),1) 862C 863 ENDDO 864 ENDDO 865 ENDDO 866 ENDDO 867C 868C--------------------------- 869C Sort integrals. 870C--------------------------- 871C 872 DO ISYML = 1, NSYM 873 ISYMAN = MULD2H(ISYANL,ISYML) 874 DO ISYMA = 1, NSYM 875 ISYMN = MULD2H(ISYMAN,ISYMA) 876 ISYMLN = MULD2H(ISYMN,ISYML) 877 ISYMBN = MULD2H(ISYMB,ISYMN) 878 ISYMAL = MULD2H(ISYMA,ISYML) 879C 880 DO N = 1, NRHF(ISYMN) 881 NBN = IT1AM(ISYMB,ISYMN) + NVIR(ISYMB)*(N-1) + B 882 DO A = 1, NVIR(ISYMA) 883 NAN = IT1AM(ISYMA,ISYMN) + NVIR(ISYMA)*(N-1) + A 884 DO L = 1, NRHF(ISYML) 885 NAL = IT1AM(ISYMA,ISYML) + NVIR(ISYMA)*(L-1) + A 886C 887 KOFF1 = IT2AM(ISYMBN,ISYMAL) + INDEX(NBN,NAL) 888 KOFF2 = KINT - 1 889 * + IMAIJA(ISYMLN,ISYMA) 890 * + NMATIJ(ISYMLN)*(A-1) 891 * + IMATIJ(ISYML,ISYMN) 892 * + NRHF(ISYML)*(N-1) 893 * + L 894 KOFF3 = KINT2 - 1 895 * + IMAIJA(ISYMLN,ISYMA) 896 * + NMATIJ(ISYMLN)*(A-1) 897 * + IMATIJ(ISYML,ISYMN) 898 * + NRHF(ISYML)*(N-1) 899 * + L 900C 901 WORK(KOFF2) = XIAJB(KOFF1) 902 WORK(KOFF3) = YIAJB(KOFF1) 903C 904 ENDDO 905 ENDDO 906 ENDDO 907C 908 ENDDO 909 ENDDO 910C 911C---------------------- 912C Construct TMAT for the g term 913C---------------------- 914C 915 DO I = 1, LENGTH 916 TMAT(I) = WMAT(INDSQ(I,2)) 917 WORK(KTMAT-1+I) = WMAT(INDSQ(I,5)) 918 ENDDO 919C 920C--------------------------------------------- 921C Symmetry sorting if symmetry 922C--------------------------------------------- 923C 924 IF (NSYM .GT. 1) THEN 925 CALL CC_GATHER(LENGTH,WORK(KEND1),TMAT,INDSQ(1,6)) 926 CALL DCOPY(LENGTH,WORK(KEND1),1,TMAT,1) 927C 928 CALL CC_GATHER(LENGTH,WORK(KEND1),WORK(KTMAT),INDSQ(1,6)) 929 CALL DCOPY(LENGTH,WORK(KEND1),1,WORK(KTMAT),1) 930 ENDIF 931C 932C------------------------------------- 933C Contract 934C------------------------------------- 935C 936 DO ISYMA = 1, NSYM 937 ISYMI = MULD2H(ISYRES,ISYMA) 938 ISYMAB = MULD2H(ISYMA,ISYMB) 939 ISYMLN = MULD2H(ISYINT,ISYMAB) 940 ISYMFM = MULD2H(ISYFIM,ISYMI) 941C 942 CALL DZERO(WORK(KEND1),NMATIJ(ISYMLN)*NRHF(ISYMI)) 943C 944 KOFF1 = ISAIKL(ISYMFM,ISYMLN) + 1 945 KOFF2 = KC2TEMP 946 * + ICKI(ISYMFM,ISYMI) 947 KOFF3 = KEND1 948C 949 NUMBFM = MAX(1,NT1AM(ISYMFM)) 950 NUMBLN = MAX(1,NMATIJ(ISYMLN)) 951C 952 CALL DGEMM('T','N',NMATIJ(ISYMLN),NRHF(ISYMI),NT1AM(ISYMFM), 953 * ONE,TMAT(KOFF1),NUMBFM,WORK(KOFF2),NUMBFM, 954 * ONE,WORK(KOFF3),NUMBLN) 955C 956 KOFF1 = KTMAT 957 * + ISAIKL(ISYMFM,ISYMLN) 958 KOFF2 = IT2SP(ISYFIM,ISYMD) 959 * + NCKI(ISYFIM)*(D-1) 960 * + ICKI(ISYMFM,ISYMI) + 1 961 KOFF3 = KEND1 962C 963 NUMBFM = MAX(1,NT1AM(ISYMFM)) 964 NUMBLN = MAX(1,NMATIJ(ISYMLN)) 965C 966 CALL DGEMM('T','N',NMATIJ(ISYMLN),NRHF(ISYMI),NT1AM(ISYMFM), 967 * ONE,WORK(KOFF1),NUMBFM,C2TP(KOFF2),NUMBFM, 968 * ONE,WORK(KOFF3),NUMBLN) 969C 970 KOFF1 = KINT 971 * + IMAIJA(ISYMLN,ISYMA) 972 KOFF2 = KEND1 973 KOFF3 = IT1AM(ISYMA,ISYMI) + 1 974C 975 NUMBLN = MAX(1,NMATIJ(ISYMLN)) 976 NUMBA = MAX(1,NVIR(ISYMA)) 977C 978 CALL DGEMM('T','N',NVIR(ISYMA),NRHF(ISYMI),NMATIJ(ISYMLN), 979 * -ONE,WORK(KOFF1),NUMBLN,WORK(KOFF2),NUMBLN, 980 * ONE,OMEGA1(KOFF3),NUMBA) 981C 982 ENDDO 983C 984C------------------------------------- 985C 3) 986C W^{fd}(emln) T^{fd}_{mi} g_{nela} 987C-------------------------------------- 988C 989C ------------------------------- 990C Sort T^{BD}_{mi} as T_{mi} 991C ------------------------------- 992C 993 ISYMMI = MULD2H(ISYMBD,ISYMC2) 994 ISYBMI = MULD2H(ISYMMI,ISYMB) 995C 996 KMIMAT = KEND1 997 KEND2 = KMIMAT + NMATIJ(ISYMMI) 998 LWRK2 = LWORK - KEND2 999C 1000 IF (LWRK2 .LT. 0) THEN 1001 CALL QUIT('Out of memory in T3_ONEL1W (sort)') 1002 ENDIF 1003C 1004 DO ISYMI = 1,NSYM 1005 ISYMM = MULD2H(ISYMMI,ISYMI) 1006 ISYMBM = MULD2H(ISYBMI,ISYMI) 1007 DO I = 1,NRHF(ISYMI) 1008 DO M = 1,NRHF(ISYMM) 1009 KOFF1 = IT2SP(ISYBMI,ISYMD) 1010 * + NCKI(ISYBMI)*(D-1) 1011 * + ICKI(ISYMBM,ISYMI) 1012 * + NT1AM(ISYMBM)*(I-1) 1013 * + IT1AM(ISYMB,ISYMM) 1014 * + NVIR(ISYMB)*(M-1) 1015 * + B 1016 KOFF2 = IMATIJ(ISYMM,ISYMI) 1017 * + NRHF(ISYMM)*(I-1) 1018 * + M 1019 WORK(KMIMAT-1+KOFF2) = C2TP(KOFF1) 1020 ENDDO 1021 ENDDO 1022 ENDDO 1023C 1024C 1025C---------------------------------- 1026C Construct TMAT 1027C---------------------------------- 1028C 1029 DO I = 1, LENGTH 1030 TMAT(I) = WMAT(INDSQ(I,5)) 1031 ENDDO 1032C 1033C------------------------------------- 1034C Omega(ai) = Omega(ai) + 1035C W^{BD}(emln) T^{BD}_{mi} g_{nela} 1036C 1037C Tmat(enl,m) * T_{mi} * Xaijb(enla) 1038C-------------------------------------- 1039C 1040 ISENLM = ISYMIM 1041 DO ISYMA = 1,NSYM 1042 ISYMI = MULD2H(ISYRES,ISYMA) 1043 ISYMMI = MULD2H(ISYMBD,ISYMC2) 1044 ISYMM = MULD2H(ISYMMI,ISYMI) 1045 ISYENL = MULD2H(ISENLM,ISYMM) 1046C 1047 KOFF1 = ISAIKJ(ISYENL,ISYMM) + 1 1048 KOFF2 = KMIMAT + IMATIJ(ISYMM,ISYMI) 1049 KOFF3 = KEND2 + ISAIKJ(ISYENL,ISYMI) 1050 NUMENL = MAX(1,NCKI(ISYENL)) 1051 NUMM = MAX(1,NRHF(ISYMM)) 1052C 1053 CALL DGEMM('N','N',NCKI(ISYENL),NRHF(ISYMI),NRHF(ISYMM), 1054 * ONE,TMAT(KOFF1),NUMENL,WORK(KOFF2),NUMM, 1055 * ZERO,WORK(KOFF3),NUMENL) 1056C 1057 KOFF1 = IT2SP(ISYENL,ISYMA) + 1 1058 KOFF2 = KEND2 + ISAIKJ(ISYENL,ISYMI) 1059 KOFF3 = IT1AM(ISYMA,ISYMI) + 1 1060 NUMENL = MAX(1,NCKI(ISYENL)) 1061 NUMA = MAX(1,NVIR(ISYMA)) 1062C 1063 CALL DGEMM('T','N',NVIR(ISYMA),NRHF(ISYMI),NCKI(ISYENL), 1064 * -ONE,XAIJB(KOFF1),NUMENL,WORK(KOFF2),NUMENL, 1065 * ONE,OMEGA1(KOFF3),NUMA) 1066 END DO 1067C 1068C------------------------------------- 1069C 4) 1070C W^{ed}(flmn) T^{df}_{mi} L_{nela} 1071C 5) 1072C W^{ed}(fmln) T^{fd}_{mi} L_{nela} 1073C------------------------------------- 1074C 1075C---------------------------------- 1076C Construct TMAT for L term 1077C---------------------------------- 1078C 1079 DO I = 1, LENGTH 1080 TMAT(I) = WMAT(INDSQ(I,1)) 1081C 1082 WORK(KTMAT-1+I) = WMAT(I) 1083 ENDDO 1084C 1085C--------------------------------------------- 1086C Symmetry sorting if symmetry 1087C--------------------------------------------- 1088C 1089 IF (NSYM .GT. 1) THEN 1090 CALL CC_GATHER(LENGTH,WORK(KEND2),TMAT,INDSQ(1,6)) 1091 CALL DCOPY(LENGTH,WORK(KEND2),1,TMAT,1) 1092C 1093 CALL CC_GATHER(LENGTH,WORK(KEND2),WORK(KTMAT),INDSQ(1,6)) 1094 CALL DCOPY(LENGTH,WORK(KEND2),1,WORK(KTMAT),1) 1095 ENDIF 1096C 1097C------------------------------------- 1098C Contract 1099C------------------------------------- 1100C 1101 DO ISYMA = 1, NSYM 1102 ISYMI = MULD2H(ISYRES,ISYMA) 1103 ISYMAB = MULD2H(ISYMA,ISYMB) 1104 ISYMLN = MULD2H(ISYINT,ISYMAB) 1105 ISYMFM = MULD2H(ISYFIM,ISYMI) 1106C 1107 CALL DZERO(WORK(KEND2),NMATIJ(ISYMLN)*NRHF(ISYMI)) 1108C 1109 KOFF1 = ISAIKL(ISYMFM,ISYMLN) + 1 1110 KOFF2 = KC2TEMP 1111 * + ICKI(ISYMFM,ISYMI) 1112 KOFF3 = KEND2 1113C 1114 NUMBFM = MAX(1,NT1AM(ISYMFM)) 1115 NUMBLN = MAX(1,NMATIJ(ISYMLN)) 1116C 1117 CALL DGEMM('T','N',NMATIJ(ISYMLN),NRHF(ISYMI),NT1AM(ISYMFM), 1118 * ONE,TMAT(KOFF1),NUMBFM,WORK(KOFF2),NUMBFM, 1119 * ONE,WORK(KOFF3),NUMBLN) 1120C 1121 KOFF1 = KTMAT 1122 * + ISAIKL(ISYMFM,ISYMLN) 1123 KOFF2 = IT2SP(ISYFIM,ISYMD) 1124 * + NCKI(ISYFIM)*(D-1) 1125 * + ICKI(ISYMFM,ISYMI) + 1 1126 KOFF3 = KEND2 1127C 1128 NUMBFM = MAX(1,NT1AM(ISYMFM)) 1129 NUMBLN = MAX(1,NMATIJ(ISYMLN)) 1130C 1131 CALL DGEMM('T','N',NMATIJ(ISYMLN),NRHF(ISYMI),NT1AM(ISYMFM), 1132 * ONE,WORK(KOFF1),NUMBFM,C2TP(KOFF2),NUMBFM, 1133 * ONE,WORK(KOFF3),NUMBLN) 1134C 1135 KOFF1 = KINT2 1136 * + IMAIJA(ISYMLN,ISYMA) 1137 KOFF2 = KEND2 1138 KOFF3 = IT1AM(ISYMA,ISYMI) + 1 1139C 1140 NUMBLN = MAX(1,NMATIJ(ISYMLN)) 1141 NUMBA = MAX(1,NVIR(ISYMA)) 1142C 1143 CALL DGEMM('T','N',NVIR(ISYMA),NRHF(ISYMI),NMATIJ(ISYMLN), 1144 * ONE,WORK(KOFF1),NUMBLN,WORK(KOFF2),NUMBLN, 1145 * ONE,OMEGA1(KOFF3),NUMBA) 1146 ENDDO 1147C 1148C-------------------------------------- 1149C 6) 1150C W^{fd}(enlm) T^{fd}_{mi} L_{nela} 1151C-------------------------------------- 1152C 1153C---------------------------------- 1154C Construct TMAT 1155C---------------------------------- 1156C 1157 DO I = 1, LENGTH 1158 TMAT(I) = WMAT(I) 1159 ENDDO 1160C 1161C------------------------------------- 1162C Omega(ai) = Omega(ai) + 1163C W^{BD}(emln) T^{BD}_{mi} L_{nela} 1164C 1165C Tmat(enl,m) * T_{mi} * Yaijb(enla) 1166C-------------------------------------- 1167C 1168 1169 ISENLM = ISYMIM 1170 DO ISYMA = 1,NSYM 1171 ISYMI = MULD2H(ISYRES,ISYMA) 1172 ISYMMI = MULD2H(ISYMBD,ISYMC2) 1173 ISYMM = MULD2H(ISYMMI,ISYMI) 1174 ISYENL = MULD2H(ISENLM,ISYMM) 1175C 1176 KOFF1 = ISAIKJ(ISYENL,ISYMM) + 1 1177 KOFF2 = KMIMAT + IMATIJ(ISYMM,ISYMI) 1178 KOFF3 = KEND2 + ISAIKJ(ISYENL,ISYMI) 1179 NUMENL = MAX(1,NCKI(ISYENL)) 1180 NUMM = MAX(1,NRHF(ISYMM)) 1181C 1182 CALL DGEMM('N','N',NCKI(ISYENL),NRHF(ISYMI),NRHF(ISYMM), 1183 * ONE,TMAT(KOFF1),NUMENL,WORK(KOFF2),NUMM, 1184 * ZERO,WORK(KOFF3),NUMENL) 1185C 1186 KOFF1 = IT2SP(ISYENL,ISYMA) + 1 1187 KOFF2 = KEND2 + ISAIKJ(ISYENL,ISYMI) 1188 KOFF3 = IT1AM(ISYMA,ISYMI) + 1 1189 NUMENL = MAX(1,NCKI(ISYENL)) 1190 NUMA = MAX(1,NVIR(ISYMA)) 1191C 1192 CALL DGEMM('T','N',NVIR(ISYMA),NRHF(ISYMI),NCKI(ISYENL), 1193 * ONE,YAIJB(KOFF1),NUMENL,WORK(KOFF2),NUMENL, 1194 * ONE,OMEGA1(KOFF3),NUMA) 1195C 1196 END DO 1197 1198C 1199 1200C 1201C---------------------------- 1202C End. 1203C---------------------------- 1204C 1205 CALL QEXIT('T3_ONEL1W') 1206C 1207 RETURN 1208 END 1209C /* Deck t3_onel2w */ 1210 SUBROUTINE T3_ONEL2W(OMEGA1,WMAT,TMAT,ISYMIM,XIAJB,YIAJB, 1211 * ISYINT,C2TP,ISYMC2,INDSQ,LENSQ,WORK,LWORK, 1212 * ISYMB,B,ISYMD,D) 1213C 1214C Written Filip Pawlowski and Poul Jorgensen based on 1215C SUBROUTINE T3_ONEL2 by K. Hald, Spring 2002. 1216C 1217C Calculate the term t^{def}_{lmn} L^{ad}_{mn} g_{ielf} 1218C -t^{def}_{nml} L^{ad}_{mn} L_{ielf} 1219C using W intermediates 1220C 1221C g contributions 1222C 1) 1223C W^{df}(enml) T^{ad}_{mn} g_{ifle} 1224C 2) 1225C W^{df}(emnl) T^{ad}_{mn} g_{ielf} 1226C 3) 1227C W^{ef}(dlnm) T^{ad}_{mn} g_{ielf} 1228C 1229C L contributions 1230C 4) 1231C W^{df}(elmn) T^{ad}_{mn} L_{ifle} 1232C 5) 1233C W^{df}(emln) T^{ad}_{mn} L_{ielf} 1234C 6) 1235C W^{ef}(dnlm) T^{ad}_{mn} L_{ielf} 1236 1237 1238C XIAJB contains g and YIAJB contains L 1239C 1240 IMPLICIT NONE 1241C 1242#include "priunit.h" 1243#include "ccsdsym.h" 1244#include "ccorb.h" 1245C 1246 INTEGER ISYMIM, ISYINT, ISYMC2, LENSQ, LWORK, ISYMB, ISYMD 1247 INTEGER INDSQ(LENSQ,6), INDEX 1248 INTEGER ISYRE1, ISYRES, ISYMBD, ISELMN, ISYAMN, ISYELI 1249 INTEGER LENGTH, KTMAT, KINT1, KINT2, KC2TEMP, KEND1, LWRK1 1250 INTEGER ISYMN, ISYMAM, ISYMA, ISYMM, ISYMMN, KOFF1, KOFF2, KOFF3 1251 INTEGER ISYML, ISYMEI, ISYMDL, ISYME, ISYMI, ISYMEL, NDL 1252 INTEGER NEI, NUMBEL, NUMBMN, NUMBA 1253 INTEGER KINTIL, ISYMIL, ISYMBI, NBI, ISYENM, NUMENM, NUML, NUMA 1254 INTEGER ISYRE2 1255C 1256#if defined (SYS_CRAY) 1257 REAL OMEGA1(*), WMAT(*), TMAT(*), XIAJB(*) 1258 REAL YIAJB(*), C2TP(*), WORK(LWORK), ZERO, ONE 1259#else 1260 DOUBLE PRECISION OMEGA1(*), WMAT(*), TMAT(*), XIAJB(*) 1261 DOUBLE PRECISION YIAJB(*), C2TP(*), WORK(LWORK), ZERO, ONE 1262#endif 1263C 1264 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 1265C 1266 INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 1267C 1268 CALL QENTER('T3_ONEL2W') 1269C 1270 ISYMBD = MULD2H(ISYMB,ISYMD) 1271C 1272 ISYRE1 = MULD2H(ISYMIM,ISYINT) 1273 ISYRE2 = MULD2H(ISYMBD,ISYRE1) 1274 ISYRES = MULD2H(ISYRE2,ISYMC2) 1275C 1276 ISELMN = ISYMIM 1277 ISYMIL = MULD2H(ISYMBD,ISYINT) 1278C 1279 ISYAMN = MULD2H(ISYMB,ISYMC2) 1280 ISYELI = MULD2H(ISYMD,ISYINT) 1281C 1282 LENGTH = NCKIJ(ISELMN) 1283C 1284 KTMAT = 1 1285 KINT1 = KTMAT + NCKIJ(ISELMN) 1286 KINT2 = KINT1 + NCKI(ISYELI) 1287 KC2TEMP = KINT2 + NCKI(ISYELI) 1288 KINTIL = KC2TEMP + NMAIJA(ISYAMN) 1289 KEND1 = KINTIL + NMATIJ(ISYMIL) 1290 LWRK1 = LWORK - KEND1 1291C 1292 IF (LWRK1 .LT. 0) THEN 1293 CALL QUIT('Out of memory in T3_ONEL2W (sort)') 1294 ENDIF 1295C 1296C Calculate the terms 1297C 1) 1298C W^{df}(enml) T^{ad}_{mn} g_{ifle} 1299C 2) 1300C W^{df}(emnl) T^{ad}_{mn} g_{ielf} 1301C 1302C----------------------------- 1303C Sort C2 1304C----------------------------- 1305C 1306 DO ISYMN = 1, NSYM 1307 ISYMAM = MULD2H(ISYAMN,ISYMN) 1308 DO ISYMA = 1, NSYM 1309 ISYMM = MULD2H(ISYMAM,ISYMA) 1310 ISYMMN = MULD2H(ISYMM,ISYMN) 1311C 1312 DO M = 1, NRHF(ISYMM) 1313 DO N = 1, NRHF(ISYMN) 1314C 1315 KOFF1 = IT2SP(ISYAMN,ISYMB) 1316 * + NCKI(ISYAMN)*(B-1) 1317 * + ICKI(ISYMAM,ISYMN) 1318 * + NT1AM(ISYMAM)*(N-1) 1319 * + IT1AM(ISYMA,ISYMM) 1320 * + NVIR(ISYMA)*(M-1) 1321 * + 1 1322C 1323 KOFF2 = KC2TEMP - 1 1324 * + IMAIJA(ISYMMN,ISYMA) 1325 * + IMATIJ(ISYMM,ISYMN) 1326 * + NRHF(ISYMM)*(N-1) 1327 * + M 1328 1329C 1330 CALL DCOPY(NVIR(ISYMA),C2TP(KOFF1),1, 1331 * WORK(KOFF2),NMATIJ(ISYMMN)) 1332C 1333 ENDDO 1334 ENDDO 1335 ENDDO 1336 ENDDO 1337C 1338C--------------------------- 1339C Sort g integrals. 1340C--------------------------- 1341C 1342 DO ISYML = 1, NSYM 1343 ISYMEI = MULD2H(ISYELI,ISYML) 1344 ISYMDL = MULD2H(ISYML,ISYMD) 1345 DO ISYME = 1, NSYM 1346 ISYMI = MULD2H(ISYMEI,ISYME) 1347 ISYMEL = MULD2H(ISYME,ISYML) 1348C 1349 DO L = 1, NRHF(ISYML) 1350 NDL = IT1AM(ISYMD,ISYML) + NVIR(ISYMD)*(L-1) + D 1351 DO E = 1, NVIR(ISYME) 1352 DO I = 1, NRHF(ISYMI) 1353 NEI = IT1AM(ISYME,ISYMI) + NVIR(ISYME)*(I-1) + E 1354C 1355 KOFF1 = IT2AM(ISYMDL,ISYMEI) + INDEX(NDL,NEI) 1356 KOFF2 = KINT1 - 1 1357 * + ICKI(ISYMEL,ISYMI) 1358 * + NT1AM(ISYMEL)*(I-1) 1359 * + IT1AM(ISYME,ISYML) 1360 * + NVIR(ISYME)*(L-1) 1361 * + E 1362 KOFF3 = KINT2 - 1 1363 * + ICKI(ISYMEI,ISYML) 1364 * + NT1AM(ISYMEI)*(L-1) 1365 * + IT1AM(ISYME,ISYMI) 1366 * + NVIR(ISYME)*(I-1) 1367 * + E 1368C 1369 WORK(KOFF2) = XIAJB(KOFF1) 1370 WORK(KOFF3) = XIAJB(KOFF1) 1371C 1372 ENDDO 1373 ENDDO 1374 ENDDO 1375C 1376 ENDDO 1377 ENDDO 1378C 1379C---------------------- 1380C Construct TMAT 1381C---------------------- 1382C 1383 DO I = 1, LENGTH 1384 TMAT(I) = WMAT(INDSQ(I,2)) 1385C 1386 WORK(KTMAT-1+I) = WMAT(INDSQ(I,5)) 1387 ENDDO 1388C 1389C--------------------------------------------- 1390C Symmetry sorting if symmetry 1391C--------------------------------------------- 1392C 1393 IF (NSYM .GT. 1) THEN 1394 CALL CC_GATHER(LENGTH,WORK(KEND1),TMAT,INDSQ(1,6)) 1395 CALL DCOPY(LENGTH,WORK(KEND1),1,TMAT,1) 1396C 1397 CALL CC_GATHER(LENGTH,WORK(KEND1),WORK(KTMAT),INDSQ(1,6)) 1398 CALL DCOPY(LENGTH,WORK(KEND1),1,WORK(KTMAT),1) 1399 ENDIF 1400C 1401C------------------------------------- 1402C Contract 1403C------------------------------------- 1404C 1405 DO ISYMA = 1, NSYM 1406 ISYMI = MULD2H(ISYRES,ISYMA) 1407 ISYMEL = MULD2H(ISYELI,ISYMI) 1408 ISYMMN = MULD2H(ISYMEL,ISELMN) 1409C 1410 KOFF1 = ISAIKL(ISYMEL,ISYMMN) + 1 1411 KOFF2 = KINT1 1412 * + ICKI(ISYMEL,ISYMI) 1413 KOFF3 = KEND1 1414C 1415 NUMBEL = MAX(1,NT1AM(ISYMEL)) 1416 NUMBMN = MAX(1,NMATIJ(ISYMMN)) 1417C 1418 CALL DGEMM('T','N',NMATIJ(ISYMMN),NRHF(ISYMI),NT1AM(ISYMEL), 1419 * ONE,TMAT(KOFF1),NUMBEL,WORK(KOFF2),NUMBEL, 1420 * ZERO,WORK(KOFF3),NUMBMN) 1421C 1422 KOFF1 = KTMAT 1423 * + ISAIKL(ISYMEL,ISYMMN) 1424 KOFF2 = KINT2 1425 * + ICKI(ISYMEL,ISYMI) 1426 KOFF3 = KEND1 1427C 1428 NUMBEL = MAX(1,NT1AM(ISYMEL)) 1429 NUMBMN = MAX(1,NMATIJ(ISYMMN)) 1430C 1431 CALL DGEMM('T','N',NMATIJ(ISYMMN),NRHF(ISYMI),NT1AM(ISYMEL), 1432 * ONE,WORK(KOFF1),NUMBEL,WORK(KOFF2),NUMBEL, 1433 * ONE,WORK(KOFF3),NUMBMN) 1434C 1435 KOFF1 = KC2TEMP 1436 * + IMAIJA(ISYMMN,ISYMA) 1437 KOFF2 = KEND1 1438 KOFF3 = IT1AM(ISYMA,ISYMI) + 1 1439C 1440 NUMBMN = MAX(1,NMATIJ(ISYMMN)) 1441 NUMBA = MAX(1,NVIR(ISYMA)) 1442C 1443 CALL DGEMM('T','N',NVIR(ISYMA),NRHF(ISYMI),NMATIJ(ISYMMN), 1444 * -ONE,WORK(KOFF1),NUMBMN,WORK(KOFF2),NUMBMN, 1445 * ONE,OMEGA1(KOFF3),NUMBA) 1446C 1447 ENDDO 1448C 1449C-------------------------------------- 1450C 3) 1451C W^{ef}(dlnm) T^{ad}_{mn} g_{ielf} 1452C-------------------------------------- 1453C 1454C------------------------------------- 1455C Sort integrals 1456C g(BiDl) sorted as g^{BD}_{li} 1457C-------------------------------------- 1458C 1459 ISYMIL = MULD2H(ISYMBD,ISYINT) 1460 DO ISYMI = 1,NSYM 1461 ISYML = MULD2H(ISYMIL,ISYMI) 1462 ISYMBI = MULD2H(ISYMB,ISYMI) 1463 ISYMDL = MULD2H(ISYMD,ISYML) 1464 DO I = 1,NRHF(ISYMI) 1465 NBI = IT1AM(ISYMB,ISYMI) + NVIR(ISYMB)*(I-1) + B 1466 DO L = 1,NRHF(ISYML) 1467 NDL = IT1AM(ISYMD,ISYML) + NVIR(ISYMD)*(L-1) + D 1468 KOFF1 = IT2AM(ISYMBI,ISYMDL) + INDEX(NBI,NDL) 1469 KOFF2 = KINTIL + IMATIJ(ISYML,ISYMI) 1470 * + NRHF(ISYML)*(I-1) + L - 1 1471 WORK(KOFF2) = XIAJB(KOFF1) 1472 ENDDO 1473 ENDDO 1474 ENDDO 1475C 1476C---------------------- 1477C Construct TMAT 1478C---------------------- 1479C 1480 DO I = 1, LENGTH 1481 TMAT(I) = WMAT(INDSQ(I,4)) 1482 ENDDO 1483C 1484C------------------------------------- 1485C Omega(ai) = Omega(ai) + 1486C W^{BD}(elnm) T^{ea}_{nm} g_{iBlD} 1487C 1488C Tmat(enml) * T_{enma} * Xiajb^{BD}_{li} 1489C-------------------------------------- 1490C 1491 DO ISYMA = 1,NSYM 1492 ISYMI = MULD2H(ISYRES,ISYMA) 1493 ISYENM = MULD2H(ISYMC2,ISYMA) 1494 ISYML = MULD2H(ISELMN,ISYENM) 1495C 1496 KOFF1 = ISAIKJ(ISYENM,ISYML) + 1 1497 KOFF2 = KINTIL + IMATIJ(ISYML,ISYMI) 1498 KOFF3 = KEND1 + ISAIKJ(ISYENM,ISYMI) 1499 NUMENM = MAX(1,NCKI(ISYENM)) 1500 NUML = MAX(1,NRHF(ISYML)) 1501 CALL DGEMM('N','N',NCKI(ISYENM),NRHF(ISYMI),NRHF(ISYML), 1502 * ONE,TMAT(KOFF1),NUMENM,WORK(KOFF2),NUML, 1503 * ZERO,WORK(KOFF3),NUMENM) 1504C 1505 KOFF1 = IT2SP(ISYENM,ISYMA) + 1 1506 KOFF2 = KEND1 + ISAIKJ(ISYENM,ISYMI) 1507 KOFF3 = IT1AM(ISYMA,ISYMI) + 1 1508 NUMENM = MAX(1,NCKI(ISYENM)) 1509 NUMA = MAX(1,NVIR(ISYMA)) 1510C 1511 CALL DGEMM('T','N',NVIR(ISYMA),NRHF(ISYMI),NCKI(ISYENM), 1512 * -ONE,C2TP(KOFF1),NUMENM,WORK(KOFF2),NUMENM, 1513 * ONE,OMEGA1(KOFF3),NUMA) 1514 ENDDO 1515C 1516C------------------------------------- 1517C 4) 1518C W^{df}(elmn) T^{ad}_{mn} L_{ifle} 1519C 5) 1520C W^{df}(emln) T^{ad}_{mn} L_{ielf} 1521C-------------------------------------- 1522C 1523C--------------------------- 1524C Sort L integrals. 1525C--------------------------- 1526C 1527 DO ISYML = 1, NSYM 1528 ISYMEI = MULD2H(ISYELI,ISYML) 1529 ISYMDL = MULD2H(ISYML,ISYMD) 1530 DO ISYME = 1, NSYM 1531 ISYMI = MULD2H(ISYMEI,ISYME) 1532 ISYMEL = MULD2H(ISYME,ISYML) 1533C 1534 DO L = 1, NRHF(ISYML) 1535 NDL = IT1AM(ISYMD,ISYML) + NVIR(ISYMD)*(L-1) + D 1536 DO E = 1, NVIR(ISYME) 1537 DO I = 1, NRHF(ISYMI) 1538 NEI = IT1AM(ISYME,ISYMI) + NVIR(ISYME)*(I-1) + E 1539C 1540 KOFF1 = IT2AM(ISYMDL,ISYMEI) + INDEX(NDL,NEI) 1541 KOFF2 = KINT1 - 1 1542 * + ICKI(ISYMEL,ISYMI) 1543 * + NT1AM(ISYMEL)*(I-1) 1544 * + IT1AM(ISYME,ISYML) 1545 * + NVIR(ISYME)*(L-1) 1546 * + E 1547 KOFF3 = KINT2 - 1 1548 * + ICKI(ISYMEI,ISYML) 1549 * + NT1AM(ISYMEI)*(L-1) 1550 * + IT1AM(ISYME,ISYMI) 1551 * + NVIR(ISYME)*(I-1) 1552 * + E 1553C 1554 WORK(KOFF2) = YIAJB(KOFF1) 1555 WORK(KOFF3) = YIAJB(KOFF1) 1556C 1557 ENDDO 1558 ENDDO 1559 ENDDO 1560C 1561 ENDDO 1562 ENDDO 1563C 1564C---------------------- 1565C Construct TMAT 1566C---------------------- 1567C 1568 DO I = 1, LENGTH 1569 TMAT(I) = WMAT(INDSQ(I,1)) 1570C 1571 WORK(KTMAT-1+I) = WMAT(I) 1572 ENDDO 1573C 1574C--------------------------------------------- 1575C Symmetry sorting if symmetry 1576C--------------------------------------------- 1577C 1578 IF (NSYM .GT. 1) THEN 1579 CALL CC_GATHER(LENGTH,WORK(KEND1),TMAT,INDSQ(1,6)) 1580 CALL DCOPY(LENGTH,WORK(KEND1),1,TMAT,1) 1581C 1582 CALL CC_GATHER(LENGTH,WORK(KEND1),WORK(KTMAT),INDSQ(1,6)) 1583 CALL DCOPY(LENGTH,WORK(KEND1),1,WORK(KTMAT),1) 1584 ENDIF 1585C 1586C------------------------------------- 1587C Contract 1588C------------------------------------- 1589C 1590 DO ISYMA = 1, NSYM 1591 ISYMI = MULD2H(ISYRES,ISYMA) 1592 ISYMEL = MULD2H(ISYELI,ISYMI) 1593 ISYMMN = MULD2H(ISYMEL,ISELMN) 1594C 1595 KOFF1 = ISAIKL(ISYMEL,ISYMMN) + 1 1596 KOFF2 = KINT1 1597 * + ICKI(ISYMEL,ISYMI) 1598 KOFF3 = KEND1 1599C 1600 NUMBEL = MAX(1,NT1AM(ISYMEL)) 1601 NUMBMN = MAX(1,NMATIJ(ISYMMN)) 1602C 1603 CALL DGEMM('T','N',NMATIJ(ISYMMN),NRHF(ISYMI),NT1AM(ISYMEL), 1604 * ONE,TMAT(KOFF1),NUMBEL,WORK(KOFF2),NUMBEL, 1605 * ZERO,WORK(KOFF3),NUMBMN) 1606C 1607 KOFF1 = KTMAT 1608 * + ISAIKL(ISYMEL,ISYMMN) 1609 KOFF2 = KINT2 1610 * + ICKI(ISYMEL,ISYMI) 1611 KOFF3 = KEND1 1612C 1613 NUMBEL = MAX(1,NT1AM(ISYMEL)) 1614 NUMBMN = MAX(1,NMATIJ(ISYMMN)) 1615C 1616 CALL DGEMM('T','N',NMATIJ(ISYMMN),NRHF(ISYMI),NT1AM(ISYMEL), 1617 * ONE,WORK(KOFF1),NUMBEL,WORK(KOFF2),NUMBEL, 1618 * ONE,WORK(KOFF3),NUMBMN) 1619C 1620 KOFF1 = KC2TEMP 1621 * + IMAIJA(ISYMMN,ISYMA) 1622 KOFF2 = KEND1 1623 KOFF3 = IT1AM(ISYMA,ISYMI) + 1 1624C 1625 NUMBMN = MAX(1,NMATIJ(ISYMMN)) 1626 NUMBA = MAX(1,NVIR(ISYMA)) 1627C 1628 CALL DGEMM('T','N',NVIR(ISYMA),NRHF(ISYMI),NMATIJ(ISYMMN), 1629 * ONE,WORK(KOFF1),NUMBMN,WORK(KOFF2),NUMBMN, 1630 * ONE,OMEGA1(KOFF3),NUMBA) 1631C 1632 ENDDO 1633C-------------------------------------- 1634C 6) 1635C W^{ef}(dnlm) T^{ad}_{mn} L_{ielf} 1636C-------------------------------------- 1637C 1638C------------------------------------- 1639C Sort integrals 1640C L(BiDl) sorted as L^{BD}_{li} 1641C-------------------------------------- 1642C 1643 ISYMIL = MULD2H(ISYMBD,ISYINT) 1644 DO ISYMI = 1,NSYM 1645 ISYML = MULD2H(ISYMIL,ISYMI) 1646 ISYMBI = MULD2H(ISYMB,ISYMI) 1647 ISYMDL = MULD2H(ISYMD,ISYML) 1648 DO I = 1,NRHF(ISYMI) 1649 NBI = IT1AM(ISYMB,ISYMI) + NVIR(ISYMB)*(I-1) + B 1650 DO L = 1,NRHF(ISYML) 1651 NDL = IT1AM(ISYMD,ISYML) + NVIR(ISYMD)*(L-1) + D 1652 KOFF1 = IT2AM(ISYMBI,ISYMDL) + INDEX(NBI,NDL) 1653 KOFF2 = KINTIL + IMATIJ(ISYML,ISYMI) 1654 * + NRHF(ISYML)*(I-1) + L - 1 1655 WORK(KOFF2) = YIAJB(KOFF1) 1656 ENDDO 1657 ENDDO 1658 ENDDO 1659C 1660 1661C 1662C---------------------- 1663C Construct TMAT 1664C---------------------- 1665C 1666 DO I = 1, LENGTH 1667 TMAT(I) = WMAT(INDSQ(I,3)) 1668 ENDDO 1669C 1670C------------------------------------- 1671C Omega(ai) = Omega(ai) + 1672C W^{BD}(enlm) T^{ea}_{nm} L_{iBlD} 1673C 1674C Tmat(enml) * T_{enma} * Yiajb^{BD}_{li} 1675C-------------------------------------- 1676C 1677 DO ISYMA = 1,NSYM 1678 ISYMI = MULD2H(ISYRES,ISYMA) 1679 ISYENM = MULD2H(ISYMC2,ISYMA) 1680 ISYML = MULD2H(ISELMN,ISYENM) 1681C 1682 KOFF1 = ISAIKJ(ISYENM,ISYML) + 1 1683 KOFF2 = KINTIL + IMATIJ(ISYML,ISYMI) 1684 KOFF3 = KEND1 1685 NUMENM = MAX(1,NCKI(ISYENM)) 1686 NUML = MAX(1,NRHF(ISYML)) 1687C 1688 CALL DGEMM('N','N',NCKI(ISYENM),NRHF(ISYMI),NRHF(ISYML), 1689 * ONE,TMAT(KOFF1),NUMENM,WORK(KOFF2),NUML, 1690 * ZERO,WORK(KOFF3),NUMENM) 1691C 1692 KOFF1 = IT2SP(ISYENM,ISYMA) + 1 1693 KOFF2 = KEND1 1694 KOFF3 = IT1AM(ISYMA,ISYMI) + 1 1695 NUMENM = MAX(1,NCKI(ISYENM)) 1696 NUMA = MAX(1,NVIR(ISYMA)) 1697C 1698 CALL DGEMM('T','N',NVIR(ISYMA),NRHF(ISYMI),NCKI(ISYENM), 1699 * ONE,C2TP(KOFF1),NUMENM,WORK(KOFF2),NUMENM, 1700 * ONE,OMEGA1(KOFF3),NUMA) 1701C 1702C 1703 ENDDO 1704C 1705C---------------------------- 1706C End. 1707C---------------------------- 1708C 1709 CALL QEXIT('T3_ONEL2W') 1710C 1711 RETURN 1712 END 1713C /* Deck t3_onel3w */ 1714 SUBROUTINE T3_ONEL3W(OMEGA1,WMAT,TMAT,ISYMIM,XIAJB,XAIBJ,ISYINT, 1715 * C2TP,ISYMC2,INDSQ,LENSQ,WORK,LWORK, 1716 * ISYMB,B,ISYMD,D) 1717C 1718C Written by Filip Pawlowski and Poul Jorgensen based on 1719C SUBROUTINE T3_ONEL3 by K. Hald, Spring 2002. 1720C 1721C Calculate the term (t^{def}_{lmn} - t^{def}_{lnm}) T^{de}_{lm} L_{ianf} 1722C based on W intermediates 1723C 1724C 1) 1725C ( W^{fe}(dlmn) - W^{fe}(dlnm) ) T^{de}_{lm} L_{ianf} 1726C 1727C + ( W^{fd}(emln) - W^{fd}(enlm) ) T^{de}_{lm} L_{ianf} 1728C 1729C 2) 1730C ( W^{de}(fnml) - W^{de}(fmnl) ) T^{de}_{lm} L_{ianf} 1731C 1732C 1733C 1734C Note : XIAJB is coming in as L and not g. 1735C 1736 IMPLICIT NONE 1737C 1738#include "priunit.h" 1739#include "ccsdsym.h" 1740#include "ccorb.h" 1741C 1742 INTEGER ISYMIM, ISYINT, ISYMC2, LENSQ, LWORK, ISYMB, ISYMD 1743 INTEGER INDSQ(LENSQ,6), INDEX 1744 INTEGER ISYRE1, ISYRES, ISYMBD, ISFLMN, ISYAIN, ISYFLM, LENGTH 1745 INTEGER KTMAT, KC2TEMP, KINT, KEND1, LWRK1, ISYMM, ISYMFL 1746 INTEGER ISYMF, ISYML, ISYMLM, ISYMFM, KOFF1, KOFF2, KOFF3 1747 INTEGER ISYMI, ISYMAN, ISYMA, ISYMN, ISYMAI, ISYMBN, NBN 1748 INTEGER NAI, NUMFLM, NUMBAI 1749 INTEGER KTLM, ISYBLM, ISYMBL, ISYMFN 1750 INTEGER NUMFN, NUMAI 1751 INTEGER ISYRE2, IS1FLM 1752C 1753#if defined (SYS_CRAY) 1754 REAL OMEGA1(*), WMAT(*), TMAT(*), XIAJB(*) 1755 REAL C2TP(*), WORK(LWORK), ZERO, ONE, XAIBJ(*) 1756 REAL DDOT 1757#else 1758 DOUBLE PRECISION OMEGA1(*), WMAT(*), TMAT(*), XIAJB(*) 1759 DOUBLE PRECISION C2TP(*), WORK(LWORK), ZERO, ONE, XAIBJ(*) 1760 DOUBLE PRECISION DDOT 1761#endif 1762C 1763 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 1764C 1765 INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 1766C 1767 CALL QENTER('T3_ONEL3W') 1768C 1769 ISYMBD = MULD2H(ISYMB,ISYMD) 1770C 1771 ISYRE1 = MULD2H(ISYMIM,ISYMC2) 1772 ISYRE2 = MULD2H(ISYRE1,ISYMBD) 1773 ISYRES = MULD2H(ISYRE2,ISYINT) 1774C 1775 ISFLMN = ISYMIM 1776 ISYAIN = MULD2H(ISYMB,ISYINT) 1777 ISYFLM = MULD2H(ISYMC2,ISYMD) 1778 ISYMLM = MULD2H(ISYMBD,ISYMC2) 1779C 1780 LENGTH = NCKIJ(ISFLMN) 1781C 1782 KTMAT = 1 1783 KC2TEMP = KTMAT + NCKIJ(ISFLMN) 1784 KINT = KC2TEMP + NCKI(ISYFLM) 1785 KTLM = KINT + NCKI(ISYAIN) 1786 KEND1 = KTLM + NMATIJ(ISYMLM) 1787 LWRK1 = LWORK - KEND1 1788C 1789C-------------------------------------------- 1790C Calculate the terms: 1791C 1) 1792C ( W^{fe}(dlmn) - W^{fe}(dlnm) ) T^{de}_{lm} L_{ianf} 1793C 1794C + ( W^{fd}(emln) - W^{fd}(enlm) ) T^{de}_{lm} L_{ianf} 1795C--------------------------------------------- 1796C 1797 IF (LWRK1 .LT. 0) THEN 1798 CALL QUIT('Out of memory in T3_ONEL3W (sort)') 1799 ENDIF 1800C 1801C----------------------------- 1802C Sort C2 1803C----------------------------- 1804C 1805 DO ISYMM = 1, NSYM 1806 ISYMFL = MULD2H(ISYFLM,ISYMM) 1807 DO ISYMF = 1, NSYM 1808 ISYML = MULD2H(ISYMFL,ISYMF) 1809 ISYMLM = MULD2H(ISYMM,ISYML) 1810 ISYMFM = MULD2H(ISYMF,ISYMM) 1811C 1812 DO M = 1, NRHF(ISYMM) 1813 DO L = 1, NRHF(ISYML) 1814C 1815 KOFF1 = IT2SP(ISYFLM,ISYMD) 1816 * + NCKI(ISYFLM)*(D-1) 1817 * + ICKI(ISYMFM,ISYML) 1818 * + NT1AM(ISYMFM)*(L-1) 1819 * + IT1AM(ISYMF,ISYMM) 1820 * + NVIR(ISYMF)*(M-1) 1821 * + 1 1822C 1823 KOFF2 = KC2TEMP - 1 1824 * + ICKI(ISYMFL,ISYMM) 1825 * + NT1AM(ISYMFL)*(M-1) 1826 * + IT1AM(ISYMF,ISYML) 1827 * + NVIR(ISYMF)*(L-1) 1828 * + 1 1829 1830C 1831 CALL DCOPY(NVIR(ISYMF),C2TP(KOFF1),1,WORK(KOFF2),1) 1832C 1833 ENDDO 1834 ENDDO 1835 ENDDO 1836 ENDDO 1837C 1838C--------------------------- 1839C Sort integrals. 1840C--------------------------- 1841C 1842 DO ISYMI = 1, NSYM 1843 ISYMAN = MULD2H(ISYAIN,ISYMI) 1844 DO ISYMA = 1, NSYM 1845 ISYMN = MULD2H(ISYMAN,ISYMA) 1846 ISYMAI = MULD2H(ISYMA,ISYMI) 1847 ISYMBN = MULD2H(ISYMB,ISYMN) 1848C 1849 DO N = 1, NRHF(ISYMN) 1850 NBN = IT1AM(ISYMB,ISYMN) + NVIR(ISYMB)*(N-1) + B 1851 DO A = 1, NVIR(ISYMA) 1852 DO I = 1, NRHF(ISYMI) 1853 NAI = IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I-1) + A 1854C 1855 KOFF1 = IT2AM(ISYMBN,ISYMAI) + INDEX(NBN,NAI) 1856 KOFF2 = KINT - 1 1857 * + ICKI(ISYMAI,ISYMN) 1858 * + NT1AM(ISYMAI)*(N-1) 1859 * + IT1AM(ISYMA,ISYMI) 1860 * + NVIR(ISYMA)*(I-1) 1861 * + A 1862C 1863 WORK(KOFF2) = XIAJB(KOFF1) 1864C 1865 ENDDO 1866 ENDDO 1867 ENDDO 1868C 1869 ENDDO 1870 ENDDO 1871C 1872C---------------------- 1873C Construct TMAT 1874C---------------------- 1875C 1876 DO I = 1, LENGTH 1877 TMAT(I) = WMAT(I) 1878 * - WMAT(INDSQ(I,3)) 1879C 1880 WORK(KTMAT-1+I) = WMAT(INDSQ(I,1)) 1881 * - WMAT(INDSQ(I,4)) 1882 ENDDO 1883C 1884C 1885C------------------------------------- 1886C Contract 1887C------------------------------------- 1888C 1889 ISYMAI = ISYRES 1890C 1891 ISYMN = MULD2H(ISYRES,ISYAIN) 1892C 1893 IS1FLM = MULD2H(ISYMIM,ISYMN) 1894 KOFF1 = ISAIKJ(IS1FLM,ISYMN) + 1 1895 KOFF2 = IT2SP(ISYFLM,ISYMD) 1896 * + NCKI(ISYFLM)*(D-1) 1897 * + 1 1898 KOFF3 = KEND1 1899C 1900 CALL DZERO(WORK(KOFF3),NRHF(ISYMN)) 1901C 1902 NUMFLM = MAX(1,NCKI(ISYFLM)) 1903C 1904C 1905 CALL DGEMV('T',NCKI(ISYFLM),NRHF(ISYMN),ONE, 1906 * TMAT(KOFF1),NUMFLM,C2TP(KOFF2),1, 1907 * ZERO,WORK(KOFF3),1) 1908C 1909C 1910 KOFF1 = KTMAT 1911 * + ISAIKJ(IS1FLM,ISYMN) 1912 KOFF2 = KC2TEMP 1913 KOFF3 = KEND1 1914C 1915 NUMFLM = MAX(1,NCKI(ISYFLM)) 1916C 1917C 1918 CALL DGEMV('T',NCKI(ISYFLM),NRHF(ISYMN),ONE, 1919 * WORK(KOFF1),NUMFLM,WORK(KOFF2),1, 1920 * ONE,WORK(KOFF3),1) 1921C 1922C 1923 KOFF1 = KINT 1924 * + ICKI(ISYMAI,ISYMN) 1925 KOFF2 = KEND1 1926 KOFF3 = 1 1927C 1928 NUMBAI = MAX(1,NT1AM(ISYMAI)) 1929C 1930C 1931 CALL DGEMV('N',NT1AM(ISYMAI),NRHF(ISYMN),-ONE, 1932 * WORK(KOFF1),NUMBAI,WORK(KOFF2),1, 1933 * ONE,OMEGA1(KOFF3),1) 1934C 1935C------------------------------------------------------------ 1936C Calculate the terms: 1937C 2) 1938C ( W^{de}(fnml) - W^{de}(fmnl) ) T^{de}_{lm} L_{ianf} 1939C------------------------------------------------------------ 1940C 1941C------------------------------------- 1942C Sort multipliers 1943C T^{BD}_{lm} sorted as T_{lm} 1944C-------------------------------------- 1945C 1946 ISYMLM = MULD2H(ISYMBD,ISYMC2) 1947 DO ISYMM = 1,NSYM 1948 ISYML = MULD2H(ISYMLM,ISYMM) 1949 ISYBLM = MULD2H(ISYMB,ISYMLM) 1950 ISYMBL = MULD2H(ISYMB,ISYML) 1951 DO L = 1,NRHF(ISYML) 1952 DO M = 1,NRHF(ISYMM) 1953C 1954 KOFF1 = IT2SP(ISYBLM,ISYMD) 1955 * + NCKI(ISYBLM)*(D-1) 1956 * + ICKI(ISYMBL,ISYMM) 1957 * + NT1AM(ISYMBL)*(M-1) 1958 * + IT1AM(ISYMB,ISYML) 1959 * + NVIR(ISYMB)*(L-1) 1960 * + B 1961C 1962 KOFF2 = KTLM + IMATIJ(ISYML,ISYMM) 1963 * + NRHF(ISYML)*(M-1) + L - 1 1964C 1965 1966 WORK(KOFF2) = C2TP(KOFF1) 1967 ENDDO 1968 ENDDO 1969 ENDDO 1970C write(lupri,*)'Norm T_{lm} = ',ddot(NMATIJ(ISYMLM), 1971C * WORK(KTLM),1,WORK(KTLM),1) 1972C 1973C---------------------- 1974C Construct TMAT 1975C---------------------- 1976C 1977 DO I = 1, LENGTH 1978 TMAT(I) = WMAT(INDSQ(I,3)) - WMAT(INDSQ(I,4)) 1979 ENDDO 1980C write(lupri,*)'Norm TMAT = ',ddot(LENGTH, 1981C * TMAT,1,TMAT,1) 1982C 1983C--------------------------------------------- 1984C Symmetry sorting if symmetry 1985C--------------------------------------------- 1986C 1987 IF (NSYM .GT. 1) THEN 1988 CALL CC_GATHER(LENGTH,WORK(KEND1),TMAT,INDSQ(1,6)) 1989 CALL DCOPY(LENGTH,WORK(KEND1),1,TMAT,1) 1990 ENDIF 1991C 1992C------------------------------------------------------ 1993C Omega(ai) = Omega(ai) + 1994C ( W^{BD}(fnml) - W^{BD}(fmnl) ) T^{BD}_{lm} L_{ianf} 1995C 1996C Tmat(fnlm) * T_{lm} * Xaibj^{aifn} 1997C------------------------------------------------------ 1998C 1999 ISYMFN = MULD2H(ISFLMN,ISYMLM) 2000 ISYMAI = MULD2H(ISYMFN,ISYINT) 2001C 2002 KOFF1 = ISAIKL(ISYMFN,ISYMLM) + 1 2003 KOFF2 = KTLM 2004 KOFF3 = KEND1 2005 NUMFN = MAX(1,NT1AM(ISYMFN)) 2006C 2007C 2008 CALL DZERO(WORK(KOFF3),NT1AM(ISYMFN)) 2009C 2010 CALL DGEMV('N',NT1AM(ISYMFN),NMATIJ(ISYMLM),ONE, 2011 * TMAT(KOFF1),NUMFN,WORK(KOFF2),1, 2012 * ZERO,WORK(KOFF3),1) 2013C 2014 KOFF1 = IT2SQ(ISYMAI,ISYMFN) + 1 2015 KOFF2 = KEND1 2016 KOFF3 = 1 2017 NUMAI = MAX(1,NT1AM(ISYMAI)) 2018C 2019 CALL DGEMV('N',NT1AM(ISYMAI),NT1AM(ISYMFN),-ONE, 2020 * XAIBJ(KOFF1),NUMAI,WORK(KOFF2),1, 2021 * ONE,OMEGA1(KOFF3),1) 2022C 2023C---------------------------- 2024C End. 2025C---------------------------- 2026C 2027 CALL QEXIT('T3_ONEL3W') 2028 RETURN 2029 END 2030