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 get_lambdax */ 20 SUBROUTINE GET_LAMBDAX(LAMPX,LAMHX,LISTX,IDLSTX,ISYMX,LAMP0,LAMH0, 21 * WORK,LWORK) 22* 23*********************************************************************** 24* * 25* Construct Lambda^X matrices. * 26* * 27* X may actually be X,Y,Z,U,..., etc. operator. * 28* * 29* Filip Pawlowski, 05-Dec-2002, Aarhus * 30* * 31*********************************************************************** 32C 33 IMPLICIT NONE 34C 35#include "priunit.h" 36#include "ccr1rsp.h" 37#include "ccsdsym.h" 38#include "ccsdinp.h" 39#include "dummy.h" 40C 41 INTEGER IDLSTX,ISYMX 42 INTEGER KT1X 43 INTEGER LWORK,KEND1,LWRK1 44 INTEGER IOPT 45C 46 CHARACTER LISTX*3, MODEL*10 47C 48#if defined (SYS_CRAY) 49 REAL LAMPX(*),LAMHX(*),LAMP0(*),LAMH0(*),WORK(LWORK) 50#else 51 DOUBLE PRECISION LAMPX(*),LAMHX(*),LAMP0(*),LAMH0(*),WORK(LWORK) 52#endif 53C 54 CALL QENTER('GET_LAMBDAX') 55C 56C---------------------------------------- 57C Initial test of symmetry consitency 58C---------------------------------------- 59C 60C IF (ISYMX .NE. ISYLRT(IDLSTX)) THEN 61C WRITE(LUPRI,*) 'ISYLRT(IDLSTX) = ',ISYLRT(IDLSTX) 62C WRITE(LUPRI,*) 'ISYMX = ',ISYMX 63C CALL QUIT('Symmetry mismatch for operator in GET_LAMBDAX') 64C END IF 65C 66C------------------- 67C Initialization 68C------------------- 69C 70 IOPT = 1 71C 72C---------------- 73C Allocation 74C---------------- 75C 76 KT1X = 1 77 KEND1 = KT1X + NT1AM(ISYMX) 78 LWRK1 = LWORK - KEND1 79C 80 IF (LWRK1 .LT. 0) THEN 81 WRITE(LUPRI,*) 'Memory available : ',LWORK 82 WRITE(LUPRI,*) 'Memory needed : ',KEND1 83 CALL QUIT('Insufficient space in GET_LAMBDAX (1)') 84 END IF 85C 86 CALL DZERO(WORK(KT1X),NT1AM(ISYMX)) 87C 88C----------------- 89C Get Lambda_X 90C----------------- 91C 92 CALL CC_RDRSP(LISTX,IDLSTX,ISYMX,IOPT,MODEL,WORK(KT1X),DUMMY) 93C 94 CALL CCLR_LAMTRA(LAMP0,LAMPX,LAMH0,LAMHX,WORK(KT1X),ISYMX) 95C 96C------------- 97C End 98C------------- 99C 100 CALL QEXIT('GET_LAMBDAX') 101C 102 RETURN 103 END 104C /* Deck get_fock0 */ 105 SUBROUTINE GET_FOCK0(FOCK0,LAMP0,LAMH0,WORK,LWORK) 106* 107*********************************************************************** 108* * 109* Read in Fock matrix in AO basis (from file) and transform to * 110* Lambda_0 basis. * 111* * 112* Could be generalized to get MO-transformed Fock matrix or * 113* f.ex. Fa,i-bar matrix, but we do not need that in CC3. * 114* * 115* Filip Pawlowski, 05-Dec-2002, Aarhus * 116* * 117*********************************************************************** 118C 119 IMPLICIT NONE 120C 121#include "priunit.h" 122#include "ccsdsym.h" 123#include "dummy.h" 124C 125 INTEGER ISYM0,LWORK,LUFOCK 126C 127#if defined (SYS_CRAY) 128 REAL FOCK0(*),LAMP0(*),LAMH0(*),WORK(LWORK) 129#else 130 DOUBLE PRECISION FOCK0(*),LAMP0(*),LAMH0(*),WORK(LWORK) 131#endif 132C 133 CALL QENTER('GET_FOCK0') 134C 135C------------------------- 136C Some initializations 137C------------------------- 138C 139 ISYM0 = 1 140C 141C----------------------------------------------- 142C Read zeroth-order AO Fock matrix from file 143C----------------------------------------------- 144C 145 LUFOCK = -1 146 CALL GPOPEN(LUFOCK,'CC_FCKH','OLD',' ','UNFORMATTED',IDUMMY, 147 * .FALSE.) 148 REWIND(LUFOCK) 149 READ(LUFOCK) (FOCK0(I),I=1,N2BST(ISYM0)) 150 CALL GPCLOSE(LUFOCK,'KEEP') 151C 152C------------------------------------------- 153C Trasform Fock matrix to Lambda_0 basis 154C------------------------------------------- 155C 156 CALL CC_FCKMO(FOCK0,LAMP0,LAMH0,WORK,LWORK,ISYM0,ISYM0,ISYM0) 157C 158C------------- 159C End 160C------------- 161C 162 CALL QEXIT('GET_FOCK0') 163C 164 RETURN 165 END 166C /* Deck get_fockx */ 167 SUBROUTINE GET_FOCKX(FOCKX,LABELX,IDLSTX,ISYMX,LAMPY,ISYMPY, 168 * LAMHZ,ISYMHZ,WORK,LWORK) 169* 170*********************************************************************** 171* * 172* Read in matrix with property integrals (FOCKX) in AO basis * 173* and transform using Lambda_pY and Lambda_hZ. 174* * 175* Poul Jorgensen, Filip Pawlowski, Spring-2003, Aarhus 176* * 177*********************************************************************** 178C 179 IMPLICIT NONE 180C 181#include "priunit.h" 182#include "ccsdsym.h" 183#include "ccl1rsp.h" 184#include "ccorb.h" 185C 186 INTEGER ISYMX,IRREP,ISYM 187 INTEGER LWORK,KEND1,LWRK1 188 INTEGER IERR 189 INTEGER IDLSTX 190 INTEGER ISYMPY,ISYMHZ 191C 192 INTEGER ISYMYZ,ISYMXYZ,LENFCKX,KFOCKX 193C 194 CHARACTER LABELX*8 195C 196#if defined (SYS_CRAY) 197 REAL FOCKX(*),LAMPY(*),LAMHZ(*),WORK(LWORK) 198#else 199 DOUBLE PRECISION FOCKX(*),LAMPY(*),LAMHZ(*),WORK(LWORK) 200#endif 201C 202 CALL QENTER('FCKX') 203C 204 ISYMYZ = MULD2H(ISYMPY,ISYMHZ) 205 ISYMXYZ = MULD2H(ISYMX,ISYMYZ) 206C 207 LENFCKX = MAX(N2BST(ISYMX),N2BST(ISYMXYZ)) 208C 209 KFOCKX = 1 210 KEND1 = KFOCKX + LENFCKX 211 LWRK1 = LWORK - KEND1 212C 213 IF (LWRK1 .LT. 0) THEN 214 WRITE(LUPRI,*)'Memory available: LWORK = ', LWORK 215 WRITE(LUPRI,*)'Memory needed: KEND1 = ', KEND1 216 CALL QUIT('Insufficient memory in GET_FOCKX') 217 END IF 218C 219C-------------------------------------- 220C Read property integrals (FOCKX) in AO basis from file 221C-------------------------------------- 222C 223 CALL CCPRPAO(LABELX,.TRUE.,WORK(KFOCKX),IRREP,ISYM,IERR, 224 * WORK(KEND1),LWRK1) 225C 226 IF ((IERR.GT.0) .OR. (IERR.EQ.0 .AND. IRREP.NE.ISYMX)) THEN 227 CALL QUIT('GET_FOCKX: error reading oper. '//LABELX) 228 ELSE IF (IERR.LT.0) THEN 229 CALL DZERO(WORK(KFOCKX),LENFCKX) 230 END IF 231C 232C--------------------------------------------------- 233C Transform property integrals (FOCKX) using Lambda_pY and Lambda_hZ 234C--------------------------------------------------- 235C 236 CALL CC_FCKMO(WORK(KFOCKX),LAMPY,LAMHZ,WORK(KEND1),LWRK1, 237 * ISYMX,ISYMPY,ISYMHZ) 238C 239C--------------------------------- 240C Copy result to FOCKX array 241C--------------------------------- 242C 243 CALL DCOPY(N2BST(ISYMXYZ),WORK(KFOCKX),1,FOCKX,1) 244C 245C------------- 246C End 247C------------- 248C 249 CALL QEXIT('FCKX') 250C 251 RETURN 252 END 253C /* Deck get_lambda0 */ 254 SUBROUTINE GET_LAMBDA0(LAMP0,LAMH0,WORK,LWORK) 255* 256*********************************************************************** 257* * 258* Construct Lambda^0 matrices. * 259* * 260* Filip Pawlowski, 05-Dec-2002, Aarhus * 261* * 262*********************************************************************** 263C 264 IMPLICIT NONE 265C 266#include "priunit.h" 267#include "ccr1rsp.h" 268#include "ccsdsym.h" 269#include "ccsdinp.h" 270#include "dummy.h" 271C 272 INTEGER IDLST0,ISYM0 273 INTEGER KT1AMP0 274 INTEGER LWORK,KEND1,LWRK1 275 INTEGER IOPT 276C 277 CHARACTER LIST0*2, MODEL*10 278C 279#if defined (SYS_CRAY) 280 REAL LAMP0(*),LAMH0(*),WORK(LWORK) 281#else 282 DOUBLE PRECISION LAMP0(*),LAMH0(*),WORK(LWORK) 283#endif 284C 285 CALL QENTER('GET_LAMBDA0') 286C 287C------------------------- 288C Some initializations 289C------------------------- 290C 291 IOPT = 1 292 LIST0 = 'R0' 293 IDLST0 = 0 294 ISYM0 = 1 295C 296C---------------- 297C Allocation 298C---------------- 299C 300 KT1AMP0 = 1 301 KEND1 = KT1AMP0 + NT1AM(ISYM0) 302 LWRK1 = LWORK - KEND1 303C 304 IF (LWRK1 .LT. 0) THEN 305 WRITE(LUPRI,*) 'Memory available : ',LWORK 306 WRITE(LUPRI,*) 'Memory needed : ',KEND1 307 CALL QUIT('Insufficient space in GET_LAMBDA0') 308 END IF 309C 310 CALL DZERO(WORK(KT1AMP0),NT1AM(ISYM0)) 311C 312 CALL CC_RDRSP(LIST0,IDLST0,ISYM0,IOPT,MODEL,WORK(KT1AMP0),DUMMY) 313C 314 CALL LAMMAT(LAMP0,LAMH0,WORK(KT1AMP0),WORK(KEND1),LWRK1) 315C 316C------------- 317C End 318C------------- 319C 320 CALL QEXIT('GET_LAMBDA0') 321C 322 RETURN 323 END 324C /* Deck get_t1_t2 */ 325 SUBROUTINE GET_T1_T2(IOPT,DIASCL,T1,T2,LIST,IDLST,ISYM,WORK,LWORK) 326* 327*********************************************************************** 328* * 329* Get T1 and/or T2 amplitudes/multipliers for a given list. * 330* * 331* IOPT = 1 : Get only T1 * 332* IOPT = 2 : Get only T2 * 333* IOPT = 3 : Get both T1 and T2 * 334* * 335* Filip Pawlowski, 06-Dec-2002, Aarhus * 336* * 337*********************************************************************** 338C 339 IMPLICIT NONE 340C 341#include "priunit.h" 342#include "ccr1rsp.h" 343#include "ccsdsym.h" 344#include "ccsdinp.h" 345#include "dummy.h" 346C 347 INTEGER IDLST,ISYM,LWORK,IOPT 348C 349 CHARACTER LIST*(*), MODEL*10 350C 351 LOGICAL DIASCL 352C 353#if defined (SYS_CRAY) 354 REAL T1(*),T2(*),WORK(LWORK) 355 REAL TWO 356#else 357 DOUBLE PRECISION T1(*),T2(*),WORK(LWORK) 358 DOUBLE PRECISION TWO 359#endif 360C 361 PARAMETER (TWO = 2.0D0) 362C 363 CALL QENTER('GET_T1_T2') 364C 365C------------------------- 366C Initial test of IOPT 367C------------------------- 368C 369 IF ((IOPT.NE.1).AND.(IOPT.NE.2).AND.(IOPT.NE.3)) THEN 370 WRITE(LUPRI,*)'IOPT = ', IOPT 371 CALL QUIT('Illegal option in GET_T1_T2') 372 END IF 373C 374C------------------------------------------------- 375C If iopt=1 get only T1 amplitudes/multipliers 376C------------------------------------------------- 377C 378 IF (IOPT.EQ.1) THEN 379C 380 CALL CC_RDRSP(LIST,IDLST,ISYM,IOPT,MODEL,T1,DUMMY) 381C 382 END IF 383C 384C------------------------------------------------- 385C If iopt=2 get only T2 amplitudes/multipliers 386C------------------------------------------------- 387C 388 IF (IOPT.EQ.2) THEN 389C 390 CALL CC_RDRSP(LIST,IDLST,ISYM,IOPT,MODEL,DUMMY,T2) 391C 392 IF (DIASCL) THEN 393 CALL CCLR_DIASCL(T2,TWO,ISYM) 394 END IF 395C 396 CALL CC_T2SQ(T2,WORK,ISYM) 397C 398 CALL CC3_T2TP(T2,WORK,ISYM) 399C 400 END IF 401C 402C--------------------------------------------------- 403C If iopt=3 get T1 and T2 amplitudes/multipliers 404C--------------------------------------------------- 405C 406 IF (IOPT.EQ.3) THEN 407C 408 CALL CC_RDRSP(LIST,IDLST,ISYM,IOPT,MODEL,T1,T2) 409C 410 IF (DIASCL) THEN 411 CALL CCLR_DIASCL(T2,TWO,ISYM) 412 END IF 413C 414 CALL CC_T2SQ(T2,WORK,ISYM) 415C 416 CALL CC3_T2TP(T2,WORK,ISYM) 417C 418 END IF 419C 420C------------- 421C End 422C------------- 423C 424 CALL QEXIT('GET_T1_T2') 425C 426 RETURN 427 END 428C /* Deck sort_fockck */ 429 SUBROUTINE SORT_FOCKCK(FOCKCK,FOCK0,ISYFCK) 430* 431*********************************************************************** 432* * 433* Sort the Fock matrix to get F(ck) block * 434* * 435* Filip Pawlowski, 06-Dec-2002, Aarhus * 436* * 437*********************************************************************** 438C 439 IMPLICIT NONE 440C 441#include "priunit.h" 442#include "ccsdsym.h" 443#include "ccorb.h" 444C 445 INTEGER ISYFCK,ISYMC,ISYMK,KOFF1,KOFF2 446C 447#if defined (SYS_CRAY) 448 REAL FOCK0(*),FOCKCK(*) 449#else 450 DOUBLE PRECISION FOCK0(*),FOCKCK(*) 451#endif 452C 453 CALL QENTER('SORT_FOCKCK') 454C 455 DO ISYMC = 1,NSYM 456C 457 ISYMK = MULD2H(ISYMC,ISYFCK) 458C 459 DO K = 1,NRHF(ISYMK) 460C 461 DO C = 1,NVIR(ISYMC) 462C 463 KOFF1 = IFCVIR(ISYMK,ISYMC) + 464 * NORB(ISYMK)*(C - 1) + K 465 KOFF2 = IT1AM(ISYMC,ISYMK) 466 * + NVIR(ISYMC)*(K - 1) + C 467C 468 FOCKCK(KOFF2) = FOCK0(KOFF1) 469C 470 ENDDO ! C 471 ENDDO ! K 472 ENDDO ! ISYMC 473C 474 CALL QEXIT('SORT_FOCKCK') 475C 476 RETURN 477 END 478C /* Deck get_orben */ 479 SUBROUTINE GET_ORBEN(FOCKD,WORK,LWORK) 480* 481*********************************************************************** 482* * 483* Get orbital energies. * 484* * 485* Filip Pawlowski, 06-Dec-2002, Aarhus * 486* * 487*********************************************************************** 488C 489 IMPLICIT NONE 490C 491#include "priunit.h" 492#include "ccsdsym.h" 493#include "dummy.h" 494#include "inftap.h" 495#include "ccorb.h" 496#include "ccsdinp.h" 497C 498 INTEGER LWORK 499C 500#if defined (SYS_CRAY) 501 REAL FOCKD(*),WORK(LWORK) 502#else 503 DOUBLE PRECISION FOCKD(*),WORK(LWORK) 504#endif 505C 506 CALL QENTER('GET_ORBEN') 507C 508C------------------------------------- 509C Read canonical orbital energies: 510C------------------------------------- 511C 512 LUSIFC = -1 513C 514 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY, 515 * .FALSE.) 516 REWIND LUSIFC 517C 518 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 519 READ (LUSIFC) 520 READ (LUSIFC) (FOCKD(I), I=1,NORBTS) 521C 522 CALL GPCLOSE(LUSIFC,'KEEP') 523C 524C--------------------------------------------- 525C Delete frozen orbitals in Fock diagonal. 526C--------------------------------------------- 527C 528 IF (FROIMP .OR. FROEXP) 529 * CALL CCSD_DELFRO(FOCKD,WORK,LWORK) 530C 531C------------- 532C End 533C------------- 534C 535 CALL QEXIT('GET_ORBEN') 536C 537 RETURN 538 END 539C /* Deck intocc_t3bar0 */ 540 SUBROUTINE INTOCC_T3BAR0(LUTOC,FNTOC,LAMH0,ISYINT,T3BOG1,T3BOL1, 541 * T3BOG2,T3BOL2,WORK,LWORK) 542* 543*********************************************************************** 544* * 545* Construct occupied integrals which are required to calculate * 546* t3bar_0 multipliers * 547* * 548* Construct 2*C-E of the integrals. * 549* Have integral for both (ij,k,a) and (a,k,j,i) * 550* * 551* * 552* Filip Pawlowski, 09-Dec-2002, Aarhus * 553* * 554*********************************************************************** 555C 556 IMPLICIT NONE 557C 558#include "priunit.h" 559#include "ccsdsym.h" 560C 561 INTEGER ISYINT,LWORK,KINTOC,KEND1,LWRK1,IOFF 562 INTEGER LUTOC 563C 564 CHARACTER*(*) FNTOC 565C 566#if defined (SYS_CRAY) 567 REAL LAMH0(*),T3BOG1(*),T3BOL1(*),T3BOG2(*),T3BOL2(*),WORK(LWORK) 568#else 569 DOUBLE PRECISION LAMH0(*),T3BOG1(*),T3BOL1(*),T3BOG2(*),T3BOL2(*) 570 DOUBLE PRECISION WORK(LWORK) 571#endif 572C 573 CALL QENTER('INTOCC_T3BAR0') 574C 575C--------------- 576C Allocation 577C--------------- 578C 579 KINTOC = 1 580 KEND1 = KINTOC + NTOTOC(ISYINT) 581 LWRK1 = LWORK - KEND1 582C 583 IF (LWRK1 .LT. 0) THEN 584 WRITE(LUPRI,*) 'Memory available : ',LWORK 585 WRITE(LUPRI,*) 'Memory needed : ',KEND1 586 CALL QUIT('Insufficient space in INTOCC_T3BAR0') 587 END IF 588C 589 CALL DZERO(WORK(KINTOC),NTOTOC(ISYINT)) 590C 591C----------------------------------------------------------- 592C Construct 2*C-E of the integrals. 593C----------------------------------------------------------- 594C 595 IOFF = 1 596C ! Read in integrals (ia|j delta) sitting as KINTOC(ah ip | jp delta) 597 IF (NTOTOC(ISYINT) .GT. 0) THEN 598 CALL GETWA2(LUTOC,FNTOC,WORK(KINTOC),IOFF,NTOTOC(ISYINT)) 599 ENDIF 600C 601C ! Transform the integrals with hole index to (ia|jk) and sort as 602C ! T3BOG1(kija) 603 CALL CC3_TROCC(WORK(KINTOC),T3BOG1,LAMH0,WORK(KEND1),LWRK1,ISYINT) 604 605C ! Calculate 2*(ia|jk) - (ka|ji) sorted as T3BOL1(kija) 606 CALL CCSDT_TCMEOCC(T3BOG1,T3BOL1,ISYINT) 607C 608C ! (ia|jk) sorted as T3BOG1(kija) and now sorted as T3BOG2(ajik) 609 CALL CCFOP_SORT(T3BOG1,T3BOG2,ISYINT,1) 610C 611 CALL CCFOP_SORT(T3BOL1,T3BOL2,ISYINT,1) 612C 613C------------- 614C End 615C------------- 616C 617 CALL QEXIT('INTOCC_T3BAR0') 618C 619 RETURN 620 END 621C /* Deck intocc_t30 */ 622 SUBROUTINE INTOCC_T30(LUCKJD,FNCKJD,LAMP0,ISYINT,T3OG1,T3OG2, 623 * WORK,LWORK) 624* 625*********************************************************************** 626* * 627* Construct occupied integrals which are required to calculate * 628* t3_0 amplitudes (in terms of SMAT and UMAT). * 629* * 630* * 631* Filip Pawlowski, 09-Dec-2002, Aarhus * 632* * 633*********************************************************************** 634C 635 IMPLICIT NONE 636C 637#include "priunit.h" 638#include "ccsdsym.h" 639C 640 INTEGER ISYINT,LWORK,KINTOC,KEND1,LWRK1,IOFF 641 INTEGER LUCKJD 642C 643 CHARACTER*(*) FNCKJD 644C 645#if defined (SYS_CRAY) 646 REAL LAMP0(*),T3OG1(*),T3OG2(*),WORK(LWORK) 647#else 648 DOUBLE PRECISION LAMP0(*),T3OG1(*),T3OG2(*),WORK(LWORK) 649#endif 650C 651 CALL QENTER('INTOCC_T30') 652C 653C--------------- 654C Allocation 655C--------------- 656C 657 KINTOC = 1 658 KEND1 = KINTOC + NTOTOC(ISYINT) 659 LWRK1 = LWORK - KEND1 660C 661 IF (LWRK1 .LT. 0) THEN 662 WRITE(LUPRI,*) 'Memory available : ',LWORK 663 WRITE(LUPRI,*) 'Memory needed : ',KEND1 664 CALL QUIT('Insufficient space in INTOCC_T30') 665 END IF 666C 667 CALL DZERO(WORK(KINTOC),NTOTOC(ISYINT)) 668C 669C------------------------ 670C Occupied integrals for S3MAT and U3MAT. 671C----------------------- 672C 673 IOFF = 1 674 IF (NTOTOC(ISYINT) .GT. 0) THEN 675 CALL GETWA2(LUCKJD,FNCKJD,WORK(KINTOC),IOFF,NTOTOC(ISYINT)) 676 ENDIF 677C 678C---------------------------------------------------------------------- 679C Transform (ai| delta j) integrals to (ai|kj) and sort as I(k,i,j,a) 680C---------------------------------------------------------------------- 681C here the xlamdp transformation need to be used 682C (delta is particle index) 683C 684 CALL CC3_TROCC(WORK(KINTOC),T3OG1,LAMP0,WORK(KEND1),LWRK1,ISYINT) 685 686C 687C---------------------------------------------------------------------- 688C (ai|kj) sorted as T3OG1(kija) becomes T3OG2(ajik) needed for U3MAT 689C---------------------------------------------------------------------- 690C 691 CALL CCFOP_SORT(T3OG1,T3OG2,ISYINT,1) 692C 693C------------- 694C End 695C------------- 696C 697 CALL QEXIT('INTOCC_T30') 698C 699 RETURN 700 END 701C /* Deck intocc_t3barx */ 702 SUBROUTINE INTOCC_T3BARX(LSKIPL1R, 703 * LUTOC,FNTOC,ISYHAM,LAMH0,ISYM0,ISYINT0, 704 * LAMHY,ISYMY,ISYINTY,W3BXOG1,W3BXOL1, 705 * W3BXOGX1,W3BXOLX1,WORK,LWORK) 706* 707*********************************************************************** 708* * 709* Construct occupied integrals which are required to calculate * 710* t3bar_X multipliers: * 711* * 712* g(ia|j k-) and L(ia|j k-) * 713* * 714* Both LambdaH_X and LambdaH_0 transformed integrals are needed * 715* * 716* ISYHAM - symmetry of Hamiltonian * 717* ISYM0 - symmetry of operator in LambdaH_0 transformation * 718* ISYINT0 - symmetry of LambdaH_0-transformed integrals * 719* ISYMY - symmetry of operator in LambdaH_Y transformation * 720* ISYINTY - symmetry of LambdaH_Y-transformed integrals * 721* * 722* * 723* Filip Pawlowski, 09-Dec-2002, Aarhus * 724* * 725*********************************************************************** 726C 727 IMPLICIT NONE 728C 729#include "priunit.h" 730#include "ccsdsym.h" 731C 732 LOGICAL LSKIPL1R 733C 734 INTEGER ISYHAM,ISYM0,ISYINT0,ISYMY,ISYINTY,LWORK,LUTOC 735 INTEGER KINTOC,KEND1,LWRK1,IOFF 736C 737 CHARACTER*(*) FNTOC 738C 739#if defined (SYS_CRAY) 740 REAL LAMH0(*),LAMHY(*),W3BXOG1(*),W3BXOL1(*),W3BXOGX1(*) 741 REAL W3BXOLX1(*),WORK(LWORK) 742#else 743 DOUBLE PRECISION LAMH0(*),LAMHY(*),W3BXOG1(*),W3BXOL1(*) 744 DOUBLE PRECISION W3BXOGX1(*),W3BXOLX1(*) 745 DOUBLE PRECISION WORK(LWORK) 746#endif 747C 748 CALL QENTER('INTOCC_T3BARX') 749C 750C--------------- 751C Allocation 752C--------------- 753C 754 KINTOC = 1 755 KEND1 = KINTOC + NTOTOC(ISYHAM) 756 LWRK1 = LWORK - KEND1 757C 758 IF (LWRK1 .LT. 0) THEN 759 WRITE(LUPRI,*) 'Memory available : ',LWORK 760 WRITE(LUPRI,*) 'Memory needed : ',KEND1 761 CALL QUIT('Insufficient space in INTOCC_T3BARX2') 762 END IF 763C 764 CALL DZERO(WORK(KINTOC),NTOTOC(ISYHAM)) 765C 766C------------------------------- 767C Read in occupied integrals 768C------------------------------- 769C 770 IOFF = 1 771 IF (NTOTOC(ISYHAM) .GT. 0) THEN 772 CALL GETWA2(LUTOC,FNTOC,WORK(KINTOC),IOFF,NTOTOC(ISYHAM)) 773 ENDIF 774C 775C----------------------------- 776C LambdaH_0 transformation 777C----------------------------- 778C 779 CALL INTOCC_T3BARX2(WORK(KINTOC),LAMH0,ISYM0,ISYINT0,W3BXOG1, 780 * W3BXOL1,WORK(KEND1),LWRK1) 781C 782C----------------------------- 783C LambdaH_Y transformation 784C----------------------------- 785C 786 IF (.NOT. LSKIPL1R) THEN 787 CALL INTOCC_T3BARX2(WORK(KINTOC),LAMHY,ISYMY,ISYINTY,W3BXOGX1, 788 * W3BXOLX1,WORK(KEND1),LWRK1) 789 END IF 790C 791C------------- 792C End 793C------------- 794C 795 CALL QEXIT('INTOCC_T3BARX') 796C 797 RETURN 798 END 799C /* Deck intocc_t3barx2 */ 800 SUBROUTINE INTOCC_T3BARX2(OCCINT,LAMH,ISYMX,ISYINT,TROCCG,TROCCL, 801 * WORK,LWORK) 802*********************************************************************** 803* * 804* Construct occupied integrals which are required to calculate * 805* t3bar_X multipliers. * 806* * 807* ISYMX - symmetry of an operator used in the transformation * 808* ISYINT - symmetry of transformed integrals * 809* * 810* * 811* Filip Pawlowski, 09-Dec-2002, Aarhus * 812* * 813*********************************************************************** 814C 815 IMPLICIT NONE 816C 817#include "priunit.h" 818C 819 INTEGER ISYMX,ISYINT,LWORK 820C 821#if defined (SYS_CRAY) 822 REAL OCCINT(*),LAMH(*),TROCCG(*),TROCCL(*),WORK(LWORK) 823#else 824 DOUBLE PRECISION OCCINT(*),LAMH(*),TROCCG(*),TROCCL(*),WORK(LWORK) 825#endif 826C 827 CALL QENTER('INTOCC_T3BARX2') 828C 829C---------------------------------------------------------------------- 830C Transform (ia|j delta) integrals to (ia|j k-) and sort G(j,i,k-,a) 831C---------------------------------------------------------------------- 832C 833 CALL CCLR_TROCC(OCCINT,TROCCG,LAMH,ISYMX,WORK,LWORK) 834C 835C---------------------------------------------------------------------- 836C integrals g(ia|j k-) sorted as G(j,i,k-,a) resorted as G(k-,i,j,a) 837C on input KTROCCG on output KTROCCG 838C L(ia|j k-) KTROCCL 839C---------------------------------------------------------------------- 840C 841 CALL CC3_SRTOCC(TROCCG,TROCCL,ISYINT) 842C 843C 844C------------- 845C End 846C------------- 847C 848 CALL QEXIT('INTOCC_T3BARX2') 849C 850 RETURN 851 END 852C /* Deck intocc_t3x */ 853 SUBROUTINE INTOCC_T3X(LUCKJDR,FNCKJDR,LAMP0,ISYINT,T3OGX,WORK, 854 * LWORK) 855* 856*********************************************************************** 857* * 858* Besides two occupied integrals constructed in INTOCC_T30 * 859* routine the additional integrals T3OGX are needed to calculate * 860* t3_X amplitudes. * 861* Those additional ocupied integrals (T3OGX) are constructed here.* 862* * 863* * 864* Filip Pawlowski, 06-Jan-2003, Aarhus * 865* * 866*********************************************************************** 867C 868 IMPLICIT NONE 869C 870#include "priunit.h" 871#include "ccsdsym.h" 872C 873 INTEGER ISYINT,LWORK,KINTOC,KEND1,LWRK1,IOFF 874 INTEGER LUCKJDR 875C 876 CHARACTER*(*) FNCKJDR 877C 878#if defined (SYS_CRAY) 879 REAL LAMP0(*),T3OGX(*),WORK(LWORK) 880#else 881 DOUBLE PRECISION LAMP0(*),T3OGX(*),WORK(LWORK) 882#endif 883C 884 CALL QENTER('INTOCC_T3X') 885C 886C--------------- 887C Allocation 888C--------------- 889C 890 KINTOC = 1 891 KEND1 = KINTOC + NTOTOC(ISYINT) 892 LWRK1 = LWORK - KEND1 893C 894 IF (LWRK1 .LT. 0) THEN 895 WRITE(LUPRI,*) 'Memory available : ',LWORK 896 WRITE(LUPRI,*) 'Memory needed : ',KEND1 897 CALL QUIT('Insufficient space in INTOCC_T3X') 898 END IF 899C 900 CALL DZERO(WORK(KINTOC),NTOTOC(ISYINT)) 901C 902C------------------------ 903C Occupied integrals for S3MAT and U3MAT. 904C----------------------- 905C 906 IOFF = 1 907 IF (NTOTOC(ISYINT) .GT. 0) THEN 908 CALL GETWA2(LUCKJDR,FNCKJDR,WORK(KINTOC),IOFF,NTOTOC(ISYINT)) 909 ENDIF 910C 911C---------------------------------------------------------------------- 912C Transform (ai|j delta) integrals to (ai|j k) and sort as (ij,k,a) 913C---------------------------------------------------------------------- 914C here the xlamdp transformation need to be used 915C (delta is particle index) 916C 917 CALL CC3_TROCC(WORK(KINTOC),T3OGX,LAMP0,WORK(KEND1),LWRK1,ISYINT) 918C 919C------------- 920C End 921C------------- 922C 923 CALL QEXIT('INTOCC_T3X') 924C 925 RETURN 926 END 927C /* Deck intvir_t30_d */ 928 SUBROUTINE INTVIR_T30_D(LUDKBC,FNDKBC,LUDELD,FNDELD,ISYINT,T3VDG1, 929 * T3VDG2,T3VDG3,LAMH,ISYMD,D,WORK,LWORK) 930* 931*********************************************************************** 932* * 933* Construct virtual integrals which are required to calculate * 934* t3_0 amplitudes (in terms of SMAT and UMAT). * 935* * 936* ISYINT - symmetry of integrals * 937* ISYMD - symmetry of external (fixed) D index) * 938* * 939* **************************************************** * 940* !!! THIS ROUTINE IS TO BE USED INSIDE THE D-LOOP !!! * 941* **************************************************** * 942* * 943* * 944* Filip Pawlowski, 09-Dec-2002, Aarhus * 945* * 946*********************************************************************** 947C 948 IMPLICIT NONE 949C 950#include "priunit.h" 951#include "ccorb.h" 952#include "ccsdsym.h" 953C 954 INTEGER ISYINT,ISYMD,LWORK,KINTVI,KEND1,LWRK1,IOFF,ISYCKB 955 INTEGER LUDKBC,LUDELD 956C 957 CHARACTER*(*) FNDKBC,FNDELD 958C 959#if defined (SYS_CRAY) 960 REAL T3VDG1(*),T3VDG2(*),T3VDG3(*),LAMH(*),WORK(LWORK) 961#else 962 DOUBLE PRECISION T3VDG1(*),T3VDG2(*),T3VDG3(*),LAMH(*),WORK(LWORK) 963#endif 964C 965 CALL QENTER('INTVIR_T30_D') 966C 967 ISYCKB = MULD2H(ISYINT,ISYMD) 968C 969C--------------- 970C Allocation 971C--------------- 972C 973 KINTVI = 1 974 KEND1 = KINTVI + NCKA(ISYCKB) 975 LWRK1 = LWORK - KEND1 976C 977 IF (LWRK1 .LT. 0) THEN 978 WRITE(LUPRI,*) 'Memory available : ',LWORK 979 WRITE(LUPRI,*) 'Memory needed : ',KEND1 980 CALL QUIT('Insufficient space in INTVIR_T30_D') 981 END IF 982C 983 CALL DZERO(WORK(KINTVI),NCKA(ISYCKB)) 984C 985C----------------------------------------------------- 986C Read virtual integrals used in first S3MAT 987C----------------------------------------------------- 988C 989 IOFF = ICKBD(ISYCKB,ISYMD) + NCKATR(ISYCKB)*(D - 1) + 1 990 IF (NCKATR(ISYCKB) .GT. 0) THEN 991 CALL GETWA2(LUDKBC,FNDKBC,T3VDG1,IOFF,NCKATR(ISYCKB)) 992 ENDIF 993C 994C----------------------------------------------------------- 995C Sort the integrals for S3MAT 996C----------------------------------------------------------- 997C 998 CALL CCSDT_SRTVIR(T3VDG1,T3VDG2,WORK(KEND1),LWRK1,ISYMD,ISYINT) 999C 1000C------------------------------------------------------ 1001C Read virtual integrals used in first U3MAT. 1002C------------------------------------------------------ 1003C 1004 IOFF = ICKAD(ISYCKB,ISYMD) + NCKA(ISYCKB)*(D - 1) + 1 1005 IF (NCKA(ISYCKB) .GT. 0) THEN 1006 CALL GETWA2(LUDELD,FNDELD,WORK(KINTVI),IOFF,NCKA(ISYCKB)) 1007 ENDIF 1008C 1009 CALL CCSDT_TRVIR(WORK(KINTVI),T3VDG3,LAMH,ISYMD,D,ISYINT, 1010 * WORK(KEND1),LWRK1) 1011C 1012C------------- 1013C End 1014C------------- 1015C 1016 CALL QEXIT('INTVIR_T30_D') 1017C 1018 RETURN 1019 END 1020C /* Deck intvir_t3bar0_d */ 1021 SUBROUTINE INTVIR_T3BAR0_D(LU3FOPX,FN3FOPX,LU3FOP2X,FN3FOP2X, 1022 * LUDKBC3,FNDKBC3,LU3VI,FN3VI,ISYINT, 1023 * T3BVDL1,T3BVDG1,T3BVDG2,T3BVDL2, 1024 * T3BVDG3,T3BVDL3,LAMP,ISYMD,D, 1025 * WORK,LWORK) 1026*********************************************************************** 1027* * 1028* Construct virtual integrals which are required to calculate * 1029* t3bar_0 amplitudes (in terms of SMAT and UMAT). * 1030* * 1031* ISYINT - symmetry of integrals * 1032* ISYMD - symmetry of external (fixed) D index) * 1033* * 1034* * 1035* ********************************************************** * 1036* !!! THIS ROUTINE IS TO BE USED INSIDE THE D- or B-LOOP !!! * 1037* ********************************************************** * 1038* * 1039* * 1040* Filip Pawlowski, 09-Dec-2002, Aarhus * 1041* * 1042*********************************************************************** 1043C 1044 IMPLICIT NONE 1045C 1046#include "priunit.h" 1047#include "ccorb.h" 1048#include "ccsdsym.h" 1049C 1050 INTEGER ISYINT,ISYMD,LWORK 1051 INTEGER LU3FOPX,LU3FOP2X,LUDKBC3,LU3VI 1052 INTEGER ISYCKB,IOFF 1053 INTEGER KINTVI,KTRVI6,KEND1,LWRK1 1054C 1055 CHARACTER*(*) FN3FOPX,FN3FOP2X,FNDKBC3,FN3VI 1056C 1057#if defined (SYS_CRAY) 1058 REAL T3BVDL1(*),T3BVDG1(*),T3BVDG2(*),T3BVDL2(*),T3BVDG3(*) 1059 REAL T3BVDL3(*),WORK(LWORK) 1060 REAL LAMP(*) 1061#else 1062 DOUBLE PRECISION T3BVDL1(*),T3BVDG1(*),T3BVDG2(*),T3BVDL2(*) 1063 DOUBLE PRECISION T3BVDG3(*),T3BVDL3(*),WORK(LWORK) 1064 DOUBLE PRECISION LAMP(*) 1065#endif 1066C 1067 CALL QENTER('INTVIR_T3BAR0_D') 1068C 1069 ISYCKB = MULD2H(ISYINT,ISYMD) 1070C 1071C ------------------------------------------- 1072C Read 2*C-E of integral used for t3-bar 1073C ------------------------------------------- 1074C 1075 IOFF = ICKBD(ISYCKB,ISYMD) + NCKATR(ISYCKB)*(D - 1) + 1 1076 IF (NCKATR(ISYCKB) .GT. 0) THEN 1077 CALL GETWA2(LU3FOP2X,FN3FOP2X,T3BVDL1,IOFF,NCKATR(ISYCKB)) 1078 ENDIF 1079C 1080C ------------------------------------------------- 1081C Integrals used for t3-bar for cc3 1082C ------------------------------------------------- 1083C 1084 IOFF = ICKBD(ISYCKB,ISYMD) + NCKATR(ISYCKB)*(D - 1) + 1 1085 IF (NCKATR(ISYCKB) .GT. 0) THEN 1086 CALL GETWA2(LUDKBC3,FNDKBC3,T3BVDG1,IOFF,NCKATR(ISYCKB)) 1087 ENDIF 1088C 1089 CALL CCSDT_SRVIR3(T3BVDG1,WORK,ISYMD,D,ISYINT) !Probably not neccesary 1090C 1091 CALL CCSDT_SRTVIR(T3BVDG1,T3BVDG2,WORK,LWORK,ISYMD,ISYINT) 1092C 1093C ------------------------------------------------ 1094C Sort the integrals for t3-bar 1095C ------------------------------------------------ 1096C 1097 CALL CCSDT_SRTVIR(T3BVDL1,T3BVDL2,WORK,LWORK,ISYMD,ISYINT) 1098C 1099C--------------- 1100C Allocation 1101C--------------- 1102C 1103 KINTVI = 1 1104 KEND1 = KINTVI + NCKA(ISYCKB) 1105 LWRK1 = LWORK - KEND1 1106C 1107 IF (LWRK1 .LT. 0) THEN 1108 CALL QUIT('Insufficient space for allocation in '// 1109 * 'INTVIR_T3BAR0_D (1)') 1110 END IF 1111C 1112 CALL DZERO(WORK(KINTVI),NCKATR(ISYCKB)) 1113C 1114C ---------------------------------------------------- 1115C Read virtual integrals used in u3am for t3-bar. 1116C ---------------------------------------------------- 1117C 1118 IOFF = ICKAD(ISYCKB,ISYMD) + NCKA(ISYCKB)*(D - 1) + 1 1119 IF (NCKA(ISYCKB) .GT. 0) THEN 1120 CALL GETWA2(LU3VI,FN3VI,WORK(KINTVI),IOFF, 1121 * NCKA(ISYCKB)) 1122 ENDIF 1123C 1124 CALL CCSDT_TRVIR(WORK(KINTVI),T3BVDG3,LAMP,ISYMD,D,ISYINT, 1125 * WORK(KEND1),LWRK1) 1126C 1127C-------------------------------------- 1128C Allocation (regain the workspace) 1129C-------------------------------------- 1130C 1131 KTRVI6 = 1 1132 KEND1 = KTRVI6 + NCKATR(ISYCKB) 1133 LWRK1 = LWORK - KEND1 1134C 1135 IF (LWRK1 .LT. 0) THEN 1136 CALL QUIT('Insufficient space for allocation in '// 1137 * 'INTVIR_T3BAR0_D (2)') 1138 END IF 1139C 1140 CALL DZERO(WORK(KTRVI6),NCKATR(ISYCKB)) 1141C 1142C---------------------- 1143C Read in integrals 1144C---------------------- 1145C 1146 IOFF = ICKBD(ISYCKB,ISYMD) + NCKATR(ISYCKB)*(D - 1) + 1 1147 IF (NCKATR(ISYCKB) .GT. 0) THEN 1148 CALL GETWA2(LU3FOPX,FN3FOPX,WORK(KTRVI6),IOFF,NCKATR(ISYCKB)) 1149 ENDIF 1150C 1151 IF (LWRK1 .LT. NCKATR(ISYCKB)) THEN 1152 CALL QUIT('Insufficient space for allocation in '// 1153 & 'INTVIR_T3BAR0_D (2)') 1154 END IF 1155 1156 CALL DCOPY(NCKATR(ISYCKB),WORK(KTRVI6),1,T3BVDL3,1) 1157C 1158C------------- 1159C End 1160C------------- 1161C 1162 CALL QEXIT('INTVIR_T3BAR0_D') 1163C 1164 RETURN 1165 END 1166C /* Deck intvir_t3barx_d */ 1167 SUBROUTINE INTVIR_T3BARX_D(LSKIPL1R, 1168 * ISYINT,LU3VI,FN3VI,LU3VI2,FN3VI2, 1169 * LU3FOP,FN3FOP,LU3FOP2,FN3FOP2, 1170 * W3BXVDGX1,W3BXVDG1,W3BXVDGX2,W3BXVDG2, 1171 * W3BXVDLX1,W3BXVDL1,W3BXVDLX2,W3BXVDL2, 1172 * LAMPY,ISYMY,LAMP0,ISYM0,ISYMD,D, 1173 * WORK,LWORK) 1174* 1175*********************************************************************** 1176* * 1177* Construct virtual integrals which are required to calculate * 1178* t3bar_X multipliers: * 1179* * 1180* * 1181* Both LambdaH_X and LambdaH_0 transformed integrals are needed * 1182* * 1183* ISYM0 - symmetry of operator in LambdaH_0 transformation * 1184* ISYINT - symmetry of integrals before transformation 1185* ISYMY - symmetry of operator in LambdaH_Y transformation * 1186* * 1187* * 1188* USE ONLY INSIDE D-LOOP 1189* * 1190* Filip Pawlowski, 09-Dec-2002, Aarhus * 1191* * 1192*********************************************************************** 1193C 1194C 1195 IMPLICIT NONE 1196C 1197#include "priunit.h" 1198#include "ccsdsym.h" 1199#include "ccorb.h" 1200C 1201 LOGICAL LSKIPL1R 1202C 1203 INTEGER ISYINT,ISYM0,ISYMY,ISYMD,LWORK 1204 INTEGER LU3VI2,LU3FOP,LU3FOP2,LU3VI 1205 INTEGER KINTVI,KEND1,LWRK1,IOFF 1206 INTEGER ISYCKB,ISYCKBY 1207C 1208 CHARACTER*(*) FN3VI2,FN3FOP,FN3FOP2,FN3VI 1209C 1210#if defined (SYS_CRAY) 1211 REAL LAMP0(*),LAMPY(*),W3BXVDGX1(*),W3BXVDG1(*),W3BXVDGX2(*) 1212 REAL W3BXVDG2(*) 1213 REAL W3BXVDLX1(*),W3BXVDL1(*),W3BXVDLX2(*),W3BXVDL2(*) 1214 REAL WORK(LWORK) 1215#else 1216 DOUBLE PRECISION LAMP0(*),LAMPY(*),W3BXVDGX1(*),W3BXVDG1(*) 1217 DOUBLE PRECISION W3BXVDLX1(*),W3BXVDL1(*),W3BXVDLX2(*),W3BXVDL2(*) 1218 DOUBLE PRECISION W3BXVDGX2(*),W3BXVDG2(*) 1219 DOUBLE PRECISION WORK(LWORK) 1220#endif 1221C 1222 CALL QENTER('INTVIR_T3BARX_D') 1223C 1224 ISYCKB = MULD2H(ISYINT,ISYMD) 1225 IF (.NOT. LSKIPL1R) THEN 1226 ISYCKBY = MULD2H(ISYMY,ISYMD) 1227 END IF 1228C 1229 KINTVI = 1 1230 IF (.NOT. LSKIPL1R) THEN 1231 KEND1 = KINTVI + MAX(NCKA(ISYCKB),NCKA(ISYCKBY)) 1232 ELSE 1233 KEND1 = KINTVI + NCKA(ISYCKB) 1234 ENDIF 1235 LWRK1 = LWORK - KEND1 1236C 1237 IF (LWRK1 .LT. 0) THEN 1238 WRITE(LUPRI,*) 'Memory available : ',LWORK 1239 WRITE(LUPRI,*) 'Memory needed : ',KEND1 1240 CALL QUIT('Insufficient space in INTVIR_T3BARX_D') 1241 END IF 1242C 1243C------------------------------------------------------------------------ 1244C 1245 IOFF = ICKAD(ISYCKB,ISYMD) + NCKA(ISYCKB)*(D - 1) + 1 1246 IF (NCKA(ISYCKB) .GT. 0) THEN 1247 CALL GETWA2(LU3VI,FN3VI,WORK(KINTVI),IOFF, 1248 * NCKA(ISYCKB)) 1249 ENDIF 1250 1251C 1252C ----------------- 1253C LAMPY TRANSFORMED 1254C ----------------- 1255C 1256C Transform g(c1^h k1^p | delta b2^h) integrals 1257C to g(c1^h k1^p | d2^pY b2^h) with lampY 1258C 1259 IF (.NOT. LSKIPL1R) THEN 1260 CALL INTVIR_T3BARX_D2(WORK(KINTVI),LAMPY,ISYMY,ISYINT,W3BXVDGX1, 1261 * ISYMD,D,WORK(KEND1),LWRK1) 1262 END IF 1263C 1264C ----------------- 1265C LAMP0 TRANSFORMED 1266C ----------------- 1267C 1268C Transform g(c1^h k1^p | delta b2^h) integrals 1269C to g(c1^h k1^p | d2^p b2^h) with lamp0 1270C 1271 CALL INTVIR_T3BARX_D2(WORK(KINTVI),LAMP0,ISYM0,ISYINT,W3BXVDG1, 1272 * ISYMD,D,WORK(KEND1),LWRK1) 1273C------------------------------------------------------------------------ 1274C 1275 IOFF = ICKAD(ISYCKB,ISYMD) + NCKA(ISYCKB)*(D - 1) + 1 1276 IF (NCKA(ISYCKB) .GT. 0) THEN 1277C 1278C Read in g(b2^h k1^p | delta c1^h) integrals and 1279C transform and sort --- see above 1280C 1281 CALL GETWA2(LU3VI2,FN3VI2,WORK(KINTVI),IOFF,NCKA(ISYCKB)) 1282 ENDIF 1283C 1284C ----------------- 1285C LAMPY TRANSFORMED 1286C ----------------- 1287C 1288C 1289 IF (.NOT. LSKIPL1R) THEN 1290 CALL INTVIR_T3BARX_D2(WORK(KINTVI),LAMPY,ISYMY,ISYINT,W3BXVDGX2, 1291 * ISYMD,D,WORK(KEND1),LWRK1) 1292 END IF 1293C 1294C ----------------- 1295C LAMP0 TRANSFORMED 1296C ----------------- 1297C 1298 CALL INTVIR_T3BARX_D2(WORK(KINTVI),LAMP0,ISYM0,ISYINT,W3BXVDG2, 1299 * ISYMD,D,WORK(KEND1),LWRK1) 1300C------------------------------------------------------------------------ 1301C 1302C 1303C Read in L(c1^h k1^p | delta b2^h) integrals 1304C 1305 CALL GETWA2(LU3FOP,FN3FOP,WORK(KINTVI),IOFF, 1306 & NCKA(ISYCKB)) 1307C 1308C ----------------- 1309C LAMPY TRANSFORMED 1310C ----------------- 1311C 1312C Transform L(c1^h k1^p | delta b2^h) integrals 1313C to L(c1^h k1^p | d2^pY b2^h) 1314C 1315 IF (.NOT. LSKIPL1R) THEN 1316 CALL INTVIR_T3BARX_D2(WORK(KINTVI),LAMPY,ISYMY,ISYINT,W3BXVDLX1, 1317 * ISYMD,D,WORK(KEND1),LWRK1) 1318 END IF 1319C 1320C ----------------- 1321C LAMP0 TRANSFORMED 1322C ----------------- 1323C 1324C Transform L(c1^h k1^p | delta b2^h) integrals 1325C to L(c1^h k1^p | d2^p b2^h) 1326C 1327 CALL INTVIR_T3BARX_D2(WORK(KINTVI),LAMP0,ISYM0,ISYINT,W3BXVDL1, 1328 * ISYMD,D,WORK(KEND1),LWRK1) 1329C------------------------------------------------------------------------ 1330C 1331C Read in L(b2^h k1^p | delta c1^h) integrals and 1332C transform and sort --- see above 1333C 1334 CALL GETWA2(LU3FOP2,FN3FOP2,WORK(KINTVI),IOFF,NCKA(ISYCKB)) 1335C 1336C 1337C ----------------- 1338C LAMPY TRANSFORMED 1339C ----------------- 1340C 1341 IF (.NOT. LSKIPL1R) THEN 1342 CALL INTVIR_T3BARX_D2(WORK(KINTVI),LAMPY,ISYMY,ISYINT,W3BXVDLX2, 1343 * ISYMD,D,WORK(KEND1),LWRK1) 1344 END IF 1345C 1346C 1347C ----------------- 1348C LAMP0 TRANSFORMED 1349C ----------------- 1350C 1351 CALL INTVIR_T3BARX_D2(WORK(KINTVI),LAMP0,ISYM0,ISYINT,W3BXVDL2, 1352 * ISYMD,D,WORK(KEND1),LWRK1) 1353C------------------------------------------------------------------------ 1354C 1355C------------- 1356C End 1357C------------- 1358C 1359 CALL QEXIT('INTVIR_T3BARX_D') 1360C 1361 RETURN 1362 END 1363C /* Deck intvir_t3barx_d2 */ 1364 SUBROUTINE INTVIR_T3BARX_D2(VIRINT,LAMP,ISYMX,ISYINT,TRVIR,ISYMD, 1365 * D,WORK,LWORK) 1366*********************************************************************** 1367* * 1368* Construct virtual integrals which are required to calculate * 1369* t3bar_X multipliers. * 1370* * 1371* ISYINT - symmetry of integrals before transformation * 1372* ISYMX - symmetry of operator * 1373* ISYINTX - symmetry of integrals after transformation * 1374* * 1375* Filip Pawlowski, 09-Dec-2002, Aarhus * 1376* * 1377*********************************************************************** 1378C 1379 IMPLICIT NONE 1380C 1381#include "priunit.h" 1382#include "ccsdsym.h" 1383#include "ccorb.h" 1384C 1385 INTEGER ISYMX,ISYINT,ISYMD,LWORK 1386 INTEGER ISYINTX 1387C 1388#if defined (SYS_CRAY) 1389 REAL VIRINT(*),LAMP(*),TRVIR(*),WORK(LWORK) 1390#else 1391 DOUBLE PRECISION VIRINT(*),LAMP(*),TRVIR(*),WORK(LWORK) 1392#endif 1393C 1394 CALL QENTER('INTVIR_T3BARX_D2') 1395C 1396 ISYINTX = MULD2H(ISYINT,ISYMX) 1397C 1398C---------------------------------------------------------------------- 1399C Transform virtual g-integrals with LambdaP 1400C---------------------------------------------------------------------- 1401C 1402 CALL CCLR_TRVIR(VIRINT,TRVIR,LAMP,ISYMX,ISYMD,D,ISYINT,WORK,LWORK) 1403C 1404C------------------- 1405C Sort integrals 1406C------------------- 1407C 1408 CALL CCSDT_SRVIR3(TRVIR,WORK,ISYMD,D,ISYINTX) 1409C 1410C 1411C------------- 1412C End 1413C------------- 1414C 1415 CALL QEXIT('INTVIR_T3BARX_D2') 1416C 1417 RETURN 1418 END 1419C /* Deck get_t30_bd */ 1420 SUBROUTINE GET_T30_BD(ISYMT3,ISYINT,T2TP,ISYMT2,T3MAT,FOCKD,DIAG, 1421 * INDSQ,LENSQ,S3MAT,T3VDG1,T3VDG2,T3OG1, 1422 * IINDEX,S3MAT3,T3VBG1,T3VBG2,INDEX2,U3MAT, 1423 * T3VDG3,T3OG2,U3MAT3,T3VBG3,ISYMB,B,ISYMD,D, 1424 * ISCKIJ,WORK,LWORK) 1425* 1426*********************************************************************** 1427* * 1428* Construct t3_0 ampllitudes for fixed B an D usin S and U * 1429* intermediates 1430* * 1431* USE ONLY INSIDE BD-LOOPS 1432* * 1433* Filip Pawlowski, 10-Dec-2002, Aarhus * 1434* * 1435*********************************************************************** 1436C 1437C 1438 IMPLICIT NONE 1439C 1440#include "priunit.h" 1441#include "ccsdsym.h" 1442#include "ccorb.h" 1443C 1444 INTEGER ISYMT3,ISYINT,ISYMT2,LENSQ,INDSQ(LENSQ,6),IINDEX,INDEX2 1445 INTEGER ISYMB,ISYMD,ISCKIJ,LWORK 1446 INTEGER ISYMBD 1447C 1448#if defined (SYS_CRAY) 1449 REAL T2TP(*),T3MAT(*),FOCKD(*),DIAG(*) 1450 REAL S3MAT(*),T3VDG1(*),T3VDG2(*),T3OG1(*) 1451 REAL S3MAT3(*),T3VBG1(*),T3VBG2(*) 1452 REAL U3MAT(*),T3VDG3(*),T3OG2(*) 1453 REAL U3MAT3(*),T3VBG3(*) 1454 REAL WORK(LWORK) 1455#else 1456 DOUBLE PRECISION T2TP(*),T3MAT(*),FOCKD(*),DIAG(*) 1457 DOUBLE PRECISION S3MAT(*),T3VDG1(*),T3VDG2(*),T3OG1(*) 1458 DOUBLE PRECISION S3MAT3(*),T3VBG1(*),T3VBG2(*) 1459 DOUBLE PRECISION U3MAT(*),T3VDG3(*),T3OG2(*) 1460 DOUBLE PRECISION U3MAT3(*),T3VBG3(*) 1461 DOUBLE PRECISION WORK(LWORK) 1462#endif 1463C 1464 CALL QENTER('GET_T30_BD') 1465C 1466C------------------------------------------ 1467C Initial check of symmetry consistency 1468C------------------------------------------ 1469C 1470 ISYMBD = MULD2H(ISYMB,ISYMD) 1471C 1472 IF (ISCKIJ .NE. MULD2H(ISYMBD,ISYMT3)) THEN 1473 WRITE(LUPRI,*)'ISCKIJ = ', ISCKIJ 1474 WRITE(LUPRI,*)'MULD2H(ISYMBD,ISYMT3) = ', MULD2H(ISYMBD,ISYMT3) 1475 CALL QUIT('Symmetry inconsistency in GET_T30_BD') 1476 END IF 1477C 1478C------------------------------------------------------------------------ 1479C Calculate the S(ci,bk,dj) matrix for T3 for B,D. (S3MAT) 1480C------------------------------------------------------------------- 1481C 1482 CALL CC3_SMAT(0.0d0,T2TP,ISYMT2,T3MAT,T3VDG1,T3VDG2,T3OG1,ISYINT, 1483 * FOCKD,DIAG,S3MAT,WORK,LWORK,IINDEX,INDSQ,LENSQ, 1484 * ISYMB,B,ISYMD,D) 1485C 1486 CALL T3_FORBIDDEN(S3MAT,ISYMT3,ISYMB,B,ISYMD,D) 1487C 1488C------------------------------------------------------------------- 1489C Calculate the S(ci,bk,dj) matrix for T3 for D,B. (S3MAT3) 1490C------------------------------------------------------------------- 1491C 1492 CALL CC3_SMAT(0.0d0,T2TP,ISYMT2,T3MAT,T3VBG1,T3VBG2,T3OG1, 1493 * ISYINT,FOCKD,DIAG,S3MAT3,WORK,LWORK,INDEX2, 1494 * INDSQ,LENSQ,ISYMD,D,ISYMB,B) 1495C 1496 CALL T3_FORBIDDEN(S3MAT3,ISYMT3,ISYMD,D,ISYMB,B) 1497C 1498C--------------------------------------------------------------------------- 1499C Calculate U(ci,jk) for fixed b,d. (U3MAT) 1500C-------------------------------------------------- 1501C 1502 CALL CC3_UMAT(0.0D0,T2TP,ISYMT2,T3VDG3,T3OG2,ISYINT,FOCKD,DIAG, 1503 * U3MAT,T3MAT,WORK,LWORK,INDSQ,LENSQ,ISYMB,B,ISYMD,D) 1504C 1505 CALL T3_FORBIDDEN(U3MAT,ISYMT3,ISYMB,B,ISYMD,D) 1506C 1507C-------------------------------------------------- 1508C Calculate U(ci,jk) for fixed d,b. (U3MAT3) 1509C-------------------------------------------------- 1510C 1511 CALL CC3_UMAT(0.0D0,T2TP,ISYMT2,T3VBG3,T3OG2,ISYINT,FOCKD,DIAG, 1512 * U3MAT3,T3MAT,WORK,LWORK,INDSQ,LENSQ,ISYMD,D,ISYMB,B) 1513C 1514 CALL T3_FORBIDDEN(U3MAT3,ISYMT3,ISYMD,D,ISYMB,B) 1515C 1516C-------------------------------------------------- 1517C Sum up S and U intermediates to get T3 BD amplitudes 1518C-------------------------------------------------- 1519C 1520 CALL CC3_T3BD(ISCKIJ,S3MAT,S3MAT3,U3MAT,U3MAT3,T3MAT,INDSQ,LENSQ) 1521C 1522C------------- 1523C End 1524C------------- 1525C 1526 CALL QEXIT('GET_T30_BD') 1527C 1528 RETURN 1529 END 1530C /* Deck get_t3bar0_bd */ 1531 SUBROUTINE GET_T3BAR0_BD(ISYMT3,T1AM,ISYMT1,T2TP,ISYMT2,TMAT, 1532 * FCKBA,FOCKD,DIAG,XIAJB,ISINT1,ISINT2, 1533 * INDSQ,LENSQ,SMAT2,T3BVDG1,T3BVDG2, 1534 * T3BVDL1,T3BVDL2,T3BOG1,T3BOL1,IINDEX, 1535 * SMAT4,T3BVBG1,T3BVBG2,T3BVBL1,T3BVBL2, 1536 * INDEX2,UMAT2,T3BVDG3,T3BVDL3,T3BOG2, 1537 * T3BOL2,UMAT4,T3BVBG3,T3BVBL3,ISYMB,B, 1538 * ISYMD,D,ISCKIJ,WORK,LWORK) 1539* 1540*********************************************************************** 1541* * 1542* Construct t3bar_0 multipliers for fixed B an D usin S and U * 1543* intermediates 1544* * 1545* USE ONLY INSIDE BD-LOOPS 1546* * 1547* Filip Pawlowski, 10-Dec-2002, Aarhus * 1548* * 1549*********************************************************************** 1550C 1551C 1552 IMPLICIT NONE 1553C 1554#include "priunit.h" 1555#include "ccsdsym.h" 1556#include "ccorb.h" 1557C 1558 INTEGER ISYMT3,ISYMT1,ISYMT2,ISINT1,ISINT2 1559 INTEGER LENSQ,INDSQ(LENSQ,6),IINDEX,INDEX2 1560 INTEGER ISYMB,ISYMD,ISCKIJ,LWORK 1561 INTEGER ISYMBD 1562C 1563#if defined (SYS_CRAY) 1564 REAL T1AM(*),T2TP(*),TMAT(*) 1565 REAL FCKBA(*),FOCKD(*),DIAG(*),XIAJB(*) 1566 REAL SMAT2(*),T3BVDG1(*),T3BVDG2(*),T3BVDL1(*),T3BVDL2(*) 1567 REAL T3BOG1(*),T3BOL1(*) 1568 REAL SMAT4(*),T3BVBG1(*),T3BVBG2(*),T3BVBL1(*),T3BVBL2(*) 1569 REAL UMAT2(*),T3BVDG3(*),T3BVDL3(*),T3BOG2(*),T3BOL2(*) 1570 REAL UMAT4(*),T3BVBG3(*),T3BVBL3(*) 1571 REAL WORK(LWORK) 1572 REAL HALF 1573#else 1574 DOUBLE PRECISION T1AM(*),T2TP(*),TMAT(*) 1575 DOUBLE PRECISION FCKBA(*),FOCKD(*),DIAG(*),XIAJB(*) 1576 DOUBLE PRECISION SMAT2(*),T3BVDG1(*),T3BVDG2(*),T3BVDL1(*) 1577 DOUBLE PRECISION T3BVDL2(*),T3BOG1(*),T3BOL1(*) 1578 DOUBLE PRECISION SMAT4(*),T3BVBG1(*),T3BVBG2(*),T3BVBL1(*) 1579 DOUBLE PRECISION T3BVBL2(*) 1580 DOUBLE PRECISION UMAT2(*),T3BVDG3(*),T3BVDL3(*),T3BOG2(*) 1581 DOUBLE PRECISION T3BOL2(*) 1582 DOUBLE PRECISION UMAT4(*),T3BVBG3(*),T3BVBL3(*) 1583 DOUBLE PRECISION WORK(LWORK) 1584 DOUBLE PRECISION HALF 1585#endif 1586C 1587 PARAMETER(HALF = 0.5D0) 1588C 1589 CALL QENTER('GET_T3BAR0_BD') 1590C 1591C------------------------------------------ 1592C Initial check of symmetry consistency 1593C------------------------------------------ 1594C 1595 ISYMBD = MULD2H(ISYMB,ISYMD) 1596C 1597 IF (ISCKIJ .NE. MULD2H(ISYMBD,ISYMT3)) THEN 1598 WRITE(LUPRI,*)'ISCKIJ = ', ISCKIJ 1599 WRITE(LUPRI,*)'MULD2H(ISYMBD,ISYMT3) = ', MULD2H(ISYMBD,ISYMT3) 1600 CALL QUIT('Symmetry inconsistency in GET_T3BAR0_BD') 1601 END IF 1602C ---------------------------------------------------- 1603C Calculate the S(ci,bk,dj) matrix for B,D for T3-BAR: 1604C ---------------------------------------------------- 1605 CALL DZERO(SMAT2,NCKIJ(ISCKIJ)) 1606C 1607 CALL CCFOP_SMAT(0.0D0,T1AM,ISYMT1,T2TP,ISYMT2,TMAT,FCKBA,XIAJB, 1608 * ISINT1,T3BVDG1,T3BVDG2,T3BVDL1,T3BVDL2,T3BOG1, 1609 * T3BOL1,ISINT2,FOCKD,DIAG,SMAT2,WORK,LWORK,IINDEX, 1610 * INDSQ,LENSQ,ISYMB,B,ISYMD,D) 1611C 1612 CALL DSCAL(NCKIJ(ISCKIJ),HALF,SMAT2,1) 1613C 1614 CALL T3_FORBIDDEN(SMAT2,ISYMT3,ISYMB,B,ISYMD,D) 1615C 1616C ---------------------------------------------------- 1617C Calculate the S(ci,bk,dj) matrix for D,B for T3-BAR: 1618C ---------------------------------------------------- 1619 CALL DZERO(SMAT4,NCKIJ(ISCKIJ)) 1620 1621 CALL CCFOP_SMAT(0.0D0,T1AM,ISYMT1,T2TP,ISYMT2,TMAT,FCKBA,XIAJB, 1622 * ISINT1,T3BVBG1,T3BVBG2,T3BVBL1,T3BVBL2,T3BOG1, 1623 * T3BOL1,ISINT2,FOCKD,DIAG,SMAT4,WORK,LWORK,INDEX2, 1624 * INDSQ,LENSQ,ISYMD,D,ISYMB,B) 1625C 1626 CALL DSCAL(NCKIJ(ISCKIJ),HALF,SMAT4,1) 1627C 1628 CALL T3_FORBIDDEN(SMAT4,ISYMT3,ISYMD,D,ISYMB,B) 1629C 1630C ------------------------------------------------ 1631C Calculate U(ci,jk) for fixed b,d for t3-bar. 1632C ------------------------------------------------ 1633 CALL DZERO(UMAT2,NCKIJ(ISCKIJ)) 1634 1635 CALL CCFOP_UMAT(0.0D0,T1AM,ISYMT1,T2TP,ISYMT2,XIAJB,ISINT1, 1636 * FCKBA,T3BVDG3,T3BVDL3,T3BOG2,T3BOL2,ISINT2,FOCKD, 1637 * DIAG,UMAT2,TMAT,WORK,LWORK,INDSQ,LENSQ,ISYMB,B, 1638 * ISYMD,D) 1639C 1640 CALL DSCAL(NCKIJ(ISCKIJ),HALF,UMAT2,1) 1641C 1642 CALL T3_FORBIDDEN(UMAT2,ISYMT3,ISYMB,B,ISYMD,D) 1643C 1644C ------------------------------------------------ 1645C Calculate U(ci,jk) for fixed d,b for t3-bar. 1646C ------------------------------------------------ 1647 CALL DZERO(UMAT4,NCKIJ(ISCKIJ)) 1648C 1649 CALL CCFOP_UMAT(0.0D0,T1AM,ISYMT1,T2TP,ISYMT2,XIAJB,ISINT1, 1650 * FCKBA,T3BVBG3,T3BVBL3,T3BOG2,T3BOL2,ISINT2,FOCKD, 1651 * DIAG,UMAT4,TMAT,WORK,LWORK,INDSQ,LENSQ,ISYMD,D, 1652 * ISYMB,B) 1653C 1654 CALL DSCAL(NCKIJ(ISCKIJ),HALF,UMAT4,1) 1655 1656 CALL T3_FORBIDDEN(UMAT4,ISYMT3,ISYMD,D,ISYMB,B) 1657C 1658C ------------------------------------------- 1659C Sum up the S-bar and U-bar to get a real T3-bar 1660C ------------------------------------------- 1661 CALL CC3_T3BD(ISCKIJ,SMAT2,SMAT4,UMAT2,UMAT4,TMAT,INDSQ,LENSQ) 1662C 1663C------------- 1664C End 1665C------------- 1666C 1667 CALL QEXIT('GET_T3BAR0_BD') 1668C 1669 RETURN 1670 END 1671C /* Deck get_t3barx_bd */ 1672 SUBROUTINE GET_T3BARX_BD(NOVIRT, 1673 * TMAT,ISTMAT,FOCKY,ISYFKY,WMAT,ISWMAT, 1674 * L2TP,ISYML2,FCKYCK,ISYFCKY,W3BXVDLX2, 1675 * W3BXVDLX1,W3BXVDGX2,W3BXVDGX1,W3BXOLX1, 1676 * W3BXOGX1,ISINT2Y,INDAJLB,INDAJLC,INDSQ, 1677 * LENSQ,L2TPY,ISYML2Y,FCKBA,ISYFCKBA, 1678 * W3BXVDL2,W3BXVDL1,W3BXVDG2,W3BXVDG1, 1679 * W3BXOL1,W3BXOG1,ISINT2,INDAJLBY,INDAJLCY, 1680 * L1Y,ISYML1Y,XIAJB,ISINT1,FREQ,DIAGW, 1681 * FOCKD,B,ISYMB,D,ISYMD,ISYMT3,WORK,LWORK) 1682* 1683*********************************************************************** 1684* * 1685* Construct t3bar_X multipliers for fixed B an D using W * 1686* intermediates 1687* * 1688* USE ONLY INSIDE BD-LOOPS 1689* * 1690* Filip Pawlowski, 10-Dec-2002, Aarhus * 1691* * 1692*********************************************************************** 1693C 1694 IMPLICIT NONE 1695C 1696#include "priunit.h" 1697#include "ccsdsym.h" 1698#include "ccorb.h" 1699C 1700 LOGICAL NOVIRT 1701 INTEGER ISTMAT,ISYFKY,ISWMAT,ISYML2,ISYFCKY,ISINT2Y,ISYML2Y 1702 INTEGER ISYFCKBA,ISINT2,ISYML1Y,ISINT1,ISYMB,ISYMD,ISYMT3 1703 INTEGER INDAJLB,INDAJLC,LENSQ,INDSQ(LENSQ,6),INDAJLBY,INDAJLCY 1704 INTEGER LWORK 1705 INTEGER ISYMBD 1706C 1707#if defined (SYS_CRAY) 1708 REAL TMAT(*),FOCKY(*),WMAT(*),L2TP(*),FCKYCK(*) 1709 REAL W3BXVDLX2(*),W3BXVDLX1(*),W3BXVDGX2(*),W3BXVDGX1(*) 1710 REAL W3BXOLX1(*),W3BXOGX1(*),L2TPY(*),FCKBA(*) 1711 REAL W3BXVDL2(*),W3BXVDL1(*),W3BXVDG2(*),W3BXVDG1(*) 1712 REAL W3BXOL1(*),W3BXOG1(*),L1Y(*),XIAJB(*) 1713 REAL DIAGW(*),FOCKD(*),WORK(LWORK) 1714 REAL FREQ 1715#else 1716 DOUBLE PRECISION TMAT(*),FOCKY(*),WMAT(*),L2TP(*),FCKYCK(*) 1717 DOUBLE PRECISION W3BXVDLX2(*),W3BXVDLX1(*) 1718 DOUBLE PRECISION W3BXVDGX2(*),W3BXVDGX1(*) 1719 DOUBLE PRECISION W3BXOLX1(*),W3BXOGX1(*),L2TPY(*),FCKBA(*) 1720 DOUBLE PRECISION W3BXVDL2(*),W3BXVDL1(*),W3BXVDG2(*),W3BXVDG1(*) 1721 DOUBLE PRECISION W3BXOL1(*),W3BXOG1(*),L1Y(*),XIAJB(*) 1722 DOUBLE PRECISION DIAGW(*),FOCKD(*),WORK(LWORK) 1723 DOUBLE PRECISION FREQ 1724#endif 1725C 1726 CALL QENTER('GET_T3BARX_BD') 1727C 1728C------------------------------------------ 1729C Initial check of symmetry consistency 1730C------------------------------------------ 1731C 1732 ISYMBD = MULD2H(ISYMB,ISYMD) 1733C 1734 IF (ISWMAT .NE. MULD2H(ISYMBD,ISYMT3)) THEN 1735 WRITE(LUPRI,*)'ISWMAT = ', ISWMAT 1736 WRITE(LUPRI,*)'MULD2H(ISYMBD,ISYMT3) = ', MULD2H(ISYMBD,ISYMT3) 1737 CALL QUIT('Symmetry inconsistency in GET_T3BARX_BD') 1738 END IF 1739C 1740C calculate <L3|[Y^,tau3]|HF> 1741C 1742 IF (.NOT.NOVIRT) THEN 1743 !virtual part 1744 CALL WBARBD_V(TMAT,ISTMAT,FOCKY,ISYFKY,WMAT,ISWMAT,WORK,LWORK) 1745 END IF 1746C 1747 !occupied part 1748 CALL WBARBD_O(TMAT,ISTMAT,FOCKY,ISYFKY,WMAT,ISWMAT,WORK,LWORK) 1749C 1750C calculate <L2|[Y,tau3]|HF> 1751C 1752 CALL WBARBD_T2(B,ISYMB,D,ISYMD,L2TP,ISYML2,FOCKY,ISYFKY,WMAT, 1753 * ISWMAT) 1754C 1755C calculate <L2|[H^Y,tau3]|HF> 1756C 1757 CALL WBARBD_TMAT(L2TP,ISYML2,WMAT,TMAT,ISWMAT,FCKYCK,ISYFCKY, 1758 * W3BXVDLX2,W3BXVDLX1,W3BXVDGX2,W3BXVDGX1,W3BXOLX1, 1759 * W3BXOGX1,ISINT2Y,WORK,LWORK,INDAJLB,INDAJLC, 1760 * INDSQ,LENSQ,ISYMB,B,ISYMD,D) 1761C 1762C calculate <L2Y|[H^,tau3]|HF> 1763C 1764 CALL WBARBD_TMAT(L2TPY,ISYML2Y,WMAT,TMAT,ISWMAT,FCKBA,ISYFCKBA, 1765 * W3BXVDL2,W3BXVDL1,W3BXVDG2,W3BXVDG1,W3BXOL1, 1766 * W3BXOG1,ISINT2,WORK,LWORK,INDAJLBY,INDAJLCY, 1767 * INDSQ,LENSQ,ISYMB,B,ISYMD,D) 1768C 1769C calculate <L1Y|[H^,tau3]|HF> 1770C 1771 CALL WBARBD_L1(L1Y,ISYML1Y,TMAT,XIAJB,ISINT1,WMAT,WORK,LWORK, 1772 * INDSQ,LENSQ,ISYMB,B,ISYMD,D) 1773C 1774C------------------------------------------------ 1775C Divide by the energy difference and 1776C remove the forbidden elements 1777C------------------------------------------------ 1778C 1779 CALL WBD_DIA(B,ISYMB,D,ISYMD,FREQ,ISWMAT,WMAT,DIAGW,FOCKD) 1780C 1781 CALL T3_FORBIDDEN(WMAT,ISYMT3,ISYMB,B,ISYMD,D) 1782C 1783C 1784C------------- 1785C End 1786C------------- 1787C 1788 CALL QEXIT('GET_T3BARX_BD') 1789C 1790 RETURN 1791 END 1792C /* Deck get_t3x_bd */ 1793 SUBROUTINE GET_T3X_BD(ISYMT3,WMAT,ISWMAT,TMAT,ISTMAT,FOCKX,ISYFKX, 1794 * T2TP,ISYMT2,T2TPX,ISYMT2X, 1795 * W3XVDG1,W3XVDG2,W3XOG1,ISYINT, 1796 * W3XVDGX1,W3XVDGX2,W3XOGX1,ISYINTX, 1797 * FREQX,DIAG,FOCKD, 1798 * INDSQ,LENSQ,B,ISYMB,D,ISYMD,WORK,LWORK) 1799* 1800*********************************************************************** 1801* * 1802* Construct t3_X amplitudes for fixed B an D using W * 1803* intermediates 1804* * 1805* USE ONLY INSIDE BD-LOOPS 1806* * 1807* Filip Pawlowski, 13-Dec-2002, Aarhus * 1808* * 1809*********************************************************************** 1810C 1811 IMPLICIT NONE 1812C 1813#include "priunit.h" 1814#include "ccsdsym.h" 1815#include "ccorb.h" 1816C 1817 INTEGER ISYMT3,ISWMAT,ISTMAT,ISYFKX,ISYMT2,ISYMT2X,ISYINT 1818 INTEGER ISYINTX,LENSQ,INDSQ(LENSQ,6),ISYMB,ISYMD,LWORK 1819 INTEGER ISYMBD 1820C 1821#if defined (SYS_CRAY) 1822 REAL WMAT(*),TMAT(*),FOCKX(*),T2TP(*),T2TPX(*),DIAG(*),FOCKD(*) 1823 REAL W3XVDG1(*),W3XVDG2(*),W3XOG1(*) 1824 REAL W3XVDGX1(*),W3XVDGX2(*),W3XOGX1(*) 1825 REAL WORK(LWORK) 1826 REAL FREQX 1827#else 1828 DOUBLE PRECISION WMAT(*),TMAT(*),FOCKX(*),T2TP(*),T2TPX(*) 1829 DOUBLE PRECISION DIAG(*),FOCKD(*) 1830 DOUBLE PRECISION W3XVDG1(*),W3XVDG2(*),W3XOG1(*) 1831 DOUBLE PRECISION W3XVDGX1(*),W3XVDGX2(*),W3XOGX1(*) 1832 DOUBLE PRECISION WORK(LWORK) 1833 DOUBLE PRECISION FREQX 1834#endif 1835C 1836 CALL QENTER('GET_T3X_BD') 1837C 1838C------------------------------------------ 1839C Initital test of symmetry consistency 1840C------------------------------------------ 1841C 1842 ISYMBD = MULD2H(ISYMB,ISYMD) 1843C 1844 IF (MULD2H(ISYMT3,ISYMBD) .NE. ISWMAT) THEN 1845 WRITE(LUPRI,*)'ISYMT3 = ', ISYMT3 1846 WRITE(LUPRI,*)'ISWMAT = ', ISWMAT 1847 CALL QUIT('Symmetry inconsistency in GET_T3X_BD') 1848 END IF 1849C 1850C------------------------------------------------------ 1851C Calculate the term <mu3|[X,T3]|HF> virtual contribution 1852C added in W^BD(KWMAT) 1853C------------------------------------------------------ 1854C 1855 CALL WBD_V(TMAT,ISTMAT,FOCKX,ISYFKX,WMAT,ISWMAT,WORK,LWORK) 1856C 1857C------------------------------------------------------ 1858C Calculate the term <mu3|[X,T3]|HF> occupied contribution 1859C added in W^BD(KWMAT) 1860C------------------------------------------------------ 1861C 1862 CALL WBD_O(TMAT,ISTMAT,FOCKX,ISYFKX,WMAT,ISWMAT,WORK,LWORK) 1863C 1864C------------------------------------------------------ 1865C Calculate the term <mu3|[[X,T2],T2]|HF> 1866C added in W^BD(KWMAT) 1867C------------------------------------------------------ 1868C 1869 CALL WBD_T2(.FALSE.,B,ISYMB,D,ISYMD,T2TP,ISYMT2,T2TP,ISYMT2, 1870 * FOCKX,ISYFKX, 1871 * INDSQ,LENSQ,WMAT,ISWMAT,WORK,LWORK) 1872C 1873C------------------------------------------------------ 1874C To get the entire T3^X add the two terms 1875C------------------------------------------------------ 1876C 1877 CALL WBD_GROUND(T2TPX,ISYMT2X,TMAT,W3XVDG1,W3XVDG2,W3XOG1,ISYINT, 1878 * WMAT,WORK,LWORK,INDSQ,LENSQ,ISYMB,B,ISYMD,D) 1879C 1880 CALL WBD_GROUND(T2TP,ISYMT2,TMAT,W3XVDGX1,W3XVDGX2,W3XOGX1, 1881 * ISYINTX,WMAT,WORK,LWORK,INDSQ,LENSQ,ISYMB,B, 1882 * ISYMD,D) 1883C 1884C------------------------------------------------ 1885C Divide by the energy difference and 1886C remove the forbidden elements 1887C------------------------------------------------ 1888C 1889 CALL WBD_DIA(B,ISYMB,D,ISYMD,FREQX,ISWMAT,WMAT,DIAG,FOCKD) 1890C 1891 CALL T3_FORBIDDEN(WMAT,ISYMT3,ISYMB,B,ISYMD,D) 1892C 1893C------------- 1894C End 1895C------------- 1896C 1897 CALL QEXIT('GET_T3X_BD') 1898C 1899 RETURN 1900 END 1901C /* Deck get_e3_bd */ 1902 SUBROUTINE GET_E3_BD(ISYMT3,WMAT,ISWMAT,TMAT, 1903 * T2TP,ISYMT2,T2TPX,ISYMT2X, 1904 * W3XVDG1,W3XVDG2,W3XOG1,ISYINT, 1905 * W3XVDGX1,W3XVDGX2,W3XOGX1,ISYINTX, 1906 * FREQX,DIAG,FOCKD, 1907 * INDSQ,LENSQ,B,ISYMB,D,ISYMD,WORK,LWORK) 1908* 1909*********************************************************************** 1910* * 1911* Construct e3 eigenvectors for fixed B an D using W * 1912* intermediates 1913* * 1914* USE ONLY INSIDE BD-LOOPS 1915* * 1916* Filip Pawlowski, 12-Nov-2003, Aarhus * 1917* * 1918*********************************************************************** 1919C 1920 IMPLICIT NONE 1921C 1922#include "priunit.h" 1923#include "ccsdsym.h" 1924#include "ccorb.h" 1925C 1926 INTEGER ISYMT3,ISWMAT,ISYMT2,ISYMT2X,ISYINT 1927 INTEGER ISYINTX,LENSQ,INDSQ(LENSQ,6),ISYMB,ISYMD,LWORK 1928 INTEGER ISYMBD 1929C 1930#if defined (SYS_CRAY) 1931 REAL WMAT(*),TMAT(*),T2TP(*),T2TPX(*),DIAG(*),FOCKD(*) 1932 REAL W3XVDG1(*),W3XVDG2(*),W3XOG1(*) 1933 REAL W3XVDGX1(*),W3XVDGX2(*),W3XOGX1(*) 1934 REAL WORK(LWORK) 1935 REAL FREQX 1936#else 1937 DOUBLE PRECISION WMAT(*),TMAT(*),T2TP(*),T2TPX(*) 1938 DOUBLE PRECISION DIAG(*),FOCKD(*) 1939 DOUBLE PRECISION W3XVDG1(*),W3XVDG2(*),W3XOG1(*) 1940 DOUBLE PRECISION W3XVDGX1(*),W3XVDGX2(*),W3XOGX1(*) 1941 DOUBLE PRECISION WORK(LWORK) 1942 DOUBLE PRECISION FREQX 1943#endif 1944C 1945 CALL QENTER('E3BD') 1946C 1947C------------------------------------------ 1948C Initital test of symmetry consistency 1949C------------------------------------------ 1950C 1951 ISYMBD = MULD2H(ISYMB,ISYMD) 1952C 1953 IF (MULD2H(ISYMT3,ISYMBD) .NE. ISWMAT) THEN 1954 WRITE(LUPRI,*)'ISYMT3 = ', ISYMT3 1955 WRITE(LUPRI,*)'ISWMAT = ', ISWMAT 1956 CALL QUIT('Symmetry inconsistency in GET_E3_BD') 1957 END IF 1958C 1959C------------------------------------------------------ 1960C <mu3|[H,E2]|HF> + <mu3|[[H,E1],T20]|HF> 1961C------------------------------------------------------ 1962C 1963 CALL WBD_GROUND(T2TPX,ISYMT2X,TMAT,W3XVDG1,W3XVDG2,W3XOG1,ISYINT, 1964 * WMAT,WORK,LWORK,INDSQ,LENSQ,ISYMB,B,ISYMD,D) 1965C 1966 CALL WBD_GROUND(T2TP,ISYMT2,TMAT,W3XVDGX1,W3XVDGX2,W3XOGX1, 1967 * ISYINTX,WMAT,WORK,LWORK,INDSQ,LENSQ,ISYMB,B, 1968 * ISYMD,D) 1969C 1970C------------------------------------------------ 1971C Divide by the energy difference and 1972C remove the forbidden elements 1973C------------------------------------------------ 1974C 1975 CALL WBD_DIA(B,ISYMB,D,ISYMD,FREQX,ISWMAT,WMAT,DIAG,FOCKD) 1976C 1977 CALL T3_FORBIDDEN(WMAT,ISYMT3,ISYMB,B,ISYMD,D) 1978C 1979C------------- 1980C End 1981C------------- 1982C 1983 CALL QEXIT('E3BD') 1984C 1985 RETURN 1986 END 1987C /* Deck get_m3bar_bd */ 1988 SUBROUTINE GET_M3BAR_BD(TMAT,WMAT,ISWMAT,L2TP,ISYML2,FCKYCK, 1989 * ISYFCKY,W3BXVDLX2,W3BXVDLX1,W3BXVDGX2, 1990 * W3BXVDGX1,W3BXOLX1,W3BXOGX1,ISINT2Y, 1991 * INDAJLB,INDAJLC,INDSQ,LENSQ,M2TP,ISYMM2, 1992 * FCKBA,ISYFCKBA,W3BXVDL2,W3BXVDL1,W3BXVDG2, 1993 * W3BXVDG1,W3BXOL1,W3BXOG1,ISINT2,INDAJLBY, 1994 * INDAJLCY,M1,ISYMM1,XIAJB,ISINT1,FREQ, 1995 * DIAGW,FOCKD,B,ISYMB,D,ISYMD,ISYMT3,WORK, 1996 * LWORK) 1997* 1998*********************************************************************** 1999* * 2000* Construct the triples part of the Mbar auxilary vector (used in * 2001* the calculation of the transition moments) fixed B an D using W * 2002* intermediates * 2003* * 2004* USE ONLY INSIDE BD-LOOPS * 2005* * 2006* Filip Pawlowski, 07-Jan-2002, Aarhus * 2007* * 2008*********************************************************************** 2009C 2010 IMPLICIT NONE 2011C 2012#include "priunit.h" 2013#include "ccsdsym.h" 2014#include "ccorb.h" 2015C 2016 INTEGER ISWMAT,ISYML2,ISYFCKY,ISINT2Y,ISYMM2 2017 INTEGER ISYFCKBA,ISINT2,ISYMM1,ISINT1,ISYMB,ISYMD,ISYMT3 2018 INTEGER INDAJLB,INDAJLC,LENSQ,INDSQ(LENSQ,6),INDAJLBY,INDAJLCY 2019 INTEGER LWORK 2020 INTEGER ISYMBD 2021C 2022#if defined (SYS_CRAY) 2023 REAL TMAT(*),WMAT(*),L2TP(*),FCKYCK(*) 2024 REAL W3BXVDLX2(*),W3BXVDLX1(*),W3BXVDGX2(*),W3BXVDGX1(*) 2025 REAL W3BXOLX1(*),W3BXOGX1(*),M2TP(*),FCKBA(*) 2026 REAL W3BXVDL2(*),W3BXVDL1(*),W3BXVDG2(*),W3BXVDG1(*) 2027 REAL W3BXOL1(*),W3BXOG1(*),M1(*),XIAJB(*) 2028 REAL DIAGW(*),FOCKD(*),WORK(LWORK) 2029 REAL FREQ 2030#else 2031 DOUBLE PRECISION TMAT(*),WMAT(*),L2TP(*),FCKYCK(*) 2032 DOUBLE PRECISION W3BXVDLX2(*),W3BXVDLX1(*) 2033 DOUBLE PRECISION W3BXVDGX2(*),W3BXVDGX1(*) 2034 DOUBLE PRECISION W3BXOLX1(*),W3BXOGX1(*),M2TP(*),FCKBA(*) 2035 DOUBLE PRECISION W3BXVDL2(*),W3BXVDL1(*),W3BXVDG2(*),W3BXVDG1(*) 2036 DOUBLE PRECISION W3BXOL1(*),W3BXOG1(*),M1(*),XIAJB(*) 2037 DOUBLE PRECISION DIAGW(*),FOCKD(*),WORK(LWORK) 2038 DOUBLE PRECISION FREQ 2039#endif 2040C 2041 CALL QENTER('GET_M3BAR_BD') 2042C 2043C------------------------------------------ 2044C Initial check of symmetry consistency 2045C------------------------------------------ 2046C 2047 ISYMBD = MULD2H(ISYMB,ISYMD) 2048C 2049 IF (ISWMAT .NE. MULD2H(ISYMBD,ISYMT3)) THEN 2050 WRITE(LUPRI,*)'ISWMAT = ', ISWMAT 2051 WRITE(LUPRI,*)'MULD2H(ISYMBD,ISYMT3) = ', MULD2H(ISYMBD,ISYMT3) 2052 CALL QUIT('Symmetry inconsistency in GET_M3BAR_BD') 2053 END IF 2054C 2055C calculate <L2|[[H,R1],tau3]|HF> 2056C 2057 CALL WBARBD_TMAT(L2TP,ISYML2,WMAT,TMAT,ISWMAT,FCKYCK,ISYFCKY, 2058 * W3BXVDLX2,W3BXVDLX1,W3BXVDGX2,W3BXVDGX1,W3BXOLX1, 2059 * W3BXOGX1,ISINT2Y,WORK,LWORK,INDAJLB,INDAJLC, 2060 * INDSQ,LENSQ,ISYMB,B,ISYMD,D) 2061C 2062C calculate <M2|[H^,tau3]|HF> 2063C 2064 CALL WBARBD_TMAT(M2TP,ISYMM2,WMAT,TMAT,ISWMAT,FCKBA,ISYFCKBA, 2065 * W3BXVDL2,W3BXVDL1,W3BXVDG2,W3BXVDG1,W3BXOL1, 2066 * W3BXOG1,ISINT2,WORK,LWORK,INDAJLBY,INDAJLCY, 2067 * INDSQ,LENSQ,ISYMB,B,ISYMD,D) 2068C 2069C calculate <M1|[H^,tau3]|HF> 2070C 2071 CALL WBARBD_L1(M1,ISYMM1,TMAT,XIAJB,ISINT1,WMAT,WORK,LWORK, 2072 * INDSQ,LENSQ,ISYMB,B,ISYMD,D) 2073 2074C 2075C------------------------------------------------ 2076C Divide by the energy difference and 2077C remove the forbidden elements 2078C------------------------------------------------ 2079C 2080 CALL WBD_DIA(B,ISYMB,D,ISYMD,FREQ,ISWMAT,WMAT,DIAGW,FOCKD) 2081C 2082 CALL T3_FORBIDDEN(WMAT,ISYMT3,ISYMB,B,ISYMD,D) 2083C 2084C------------- 2085C End 2086C------------- 2087C 2088 CALL QEXIT('GET_M3BAR_BD') 2089C 2090 RETURN 2091 END 2092C /* Deck get_l3bar_bd */ 2093 SUBROUTINE GET_L3BAR_BD(TMAT,WMAT,ISWMAT,INDSQ,LENSQ,L2TP,ISYML2, 2094 * FCKBA,ISYFCKBA,W3BXVDL2,W3BXVDL1,W3BXVDG2, 2095 * W3BXVDG1,W3BXOL1,W3BXOG1,ISINT2,INDAJLBY, 2096 * INDAJLCY,L1,ISYML1,XIAJB,ISINT1,FREQ, 2097 * DIAGW,FOCKD,B,ISYMB,D,ISYMD,ISYMT3,WORK, 2098 * LWORK) 2099* 2100*********************************************************************** 2101* * 2102* Construct the triples part of the left eigenvector of the * 2103* Jacobian for fixed B an D using W intermediates * 2104* * 2105* USE ONLY INSIDE BD-LOOPS * 2106* * 2107* Filip Pawlowski, 07-Jan-2002, Aarhus * 2108* * 2109*********************************************************************** 2110C 2111 IMPLICIT NONE 2112C 2113#include "priunit.h" 2114#include "ccsdsym.h" 2115#include "ccorb.h" 2116C 2117 INTEGER ISWMAT,ISYML2 2118 INTEGER ISYFCKBA,ISINT2,ISYML1,ISINT1,ISYMB,ISYMD,ISYMT3 2119 INTEGER LENSQ,INDSQ(LENSQ,6),INDAJLBY,INDAJLCY 2120 INTEGER LWORK 2121 INTEGER ISYMBD 2122C 2123#if defined (SYS_CRAY) 2124 REAL TMAT(*),WMAT(*) 2125 REAL L2TP(*),FCKBA(*) 2126 REAL W3BXVDL2(*),W3BXVDL1(*),W3BXVDG2(*),W3BXVDG1(*) 2127 REAL W3BXOL1(*),W3BXOG1(*),L1(*),XIAJB(*) 2128 REAL DIAGW(*),FOCKD(*),WORK(LWORK) 2129 REAL FREQ 2130#else 2131 DOUBLE PRECISION TMAT(*),WMAT(*) 2132 DOUBLE PRECISION L2TP(*),FCKBA(*) 2133 DOUBLE PRECISION W3BXVDL2(*),W3BXVDL1(*),W3BXVDG2(*),W3BXVDG1(*) 2134 DOUBLE PRECISION W3BXOL1(*),W3BXOG1(*),L1(*),XIAJB(*) 2135 DOUBLE PRECISION DIAGW(*),FOCKD(*),WORK(LWORK) 2136 DOUBLE PRECISION FREQ 2137#endif 2138C 2139 CALL QENTER('GET_L3BAR_BD') 2140C 2141C------------------------------------------ 2142C Initial check of symmetry consistency 2143C------------------------------------------ 2144C 2145 ISYMBD = MULD2H(ISYMB,ISYMD) 2146C 2147 IF (ISWMAT .NE. MULD2H(ISYMBD,ISYMT3)) THEN 2148 WRITE(LUPRI,*)'ISWMAT = ', ISWMAT 2149 WRITE(LUPRI,*)'MULD2H(ISYMBD,ISYMT3) = ', MULD2H(ISYMBD,ISYMT3) 2150 CALL QUIT('Symmetry inconsistency in GET_L3BAR_BD') 2151 END IF 2152C 2153C calculate <L2|[H^,tau3]|HF> 2154C 2155 CALL WBARBD_TMAT(L2TP,ISYML2,WMAT,TMAT,ISWMAT,FCKBA,ISYFCKBA, 2156 * W3BXVDL2,W3BXVDL1,W3BXVDG2,W3BXVDG1,W3BXOL1, 2157 * W3BXOG1,ISINT2,WORK,LWORK,INDAJLBY,INDAJLCY, 2158 * INDSQ,LENSQ,ISYMB,B,ISYMD,D) 2159C 2160C calculate <L1|[H^,tau3]|HF> 2161C 2162 CALL WBARBD_L1(L1,ISYML1,TMAT,XIAJB,ISINT1,WMAT,WORK,LWORK, 2163 * INDSQ,LENSQ,ISYMB,B,ISYMD,D) 2164C 2165C------------------------------------------------ 2166C Divide by the energy difference and 2167C remove the forbidden elements 2168C------------------------------------------------ 2169C 2170 CALL WBD_DIA(B,ISYMB,D,ISYMD,FREQ,ISWMAT,WMAT,DIAGW,FOCKD) 2171C 2172 CALL T3_FORBIDDEN(WMAT,ISYMT3,ISYMB,B,ISYMD,D) 2173C 2174C------------- 2175C End 2176C------------- 2177C 2178 CALL QEXIT('GET_L3BAR_BD') 2179C 2180 RETURN 2181 END 2182C /* Deck write_t3_dl */ 2183 SUBROUTINE WRITE_T3_DL(LUFILE,FNFILE,T3,ISYMT3,ISYMD,ISYMB,B) 2184*************************************************** 2185* 2186* Write T3 amplitudes as T3^D(ai,bj,l) to disc. 2187* 2188* F. Pawlowski, 25-02-2003, Aarhus. 2189*************************************************** 2190C 2191 IMPLICIT NONE 2192C 2193#include "ccsdsym.h" 2194#include "ccorb.h" 2195#include "cc3t3d.h" 2196C 2197 INTEGER LUFILE 2198 INTEGER ISYMT3,ISYMD,ISYMB 2199 INTEGER ISYML,ISYMDL,ISAIBJ,ISYMJ,ISYMBJ,ISYMAI,ISYAIL 2200 INTEGER KOFF1,NBJ,IADR 2201C 2202 CHARACTER*(*) FNFILE 2203C 2204#if defined (SYS_CRAY) 2205 REAL T3(*) 2206#else 2207 DOUBLE PRECISION T3(*) 2208#endif 2209C 2210 CALL QENTER('WRITE_T3_DL') 2211C 2212 DO ISYML = 1, NSYM 2213 ISYMDL = MULD2H(ISYMD,ISYML) 2214 ISAIBJ = MULD2H(ISYMT3,ISYMDL) 2215 DO L = 1, NRHF(ISYML) 2216 DO ISYMJ = 1, NSYM 2217 ISYMBJ = MULD2H(ISYMJ,ISYMB) 2218 ISYMAI = MULD2H(ISAIBJ,ISYMBJ) 2219 ISYAIL = MULD2H(ISYMAI,ISYML) 2220 DO J = 1, NRHF(ISYMJ) 2221 2222 KOFF1 = ISAIKJ(ISYAIL,ISYMJ)+ NCKI(ISYAIL)*(J-1) 2223 * + ISAIK(ISYMAI, ISYML)+NT1AM(ISYMAI)*(L-1) 2224 * + 1 2225 2226 NBJ = IT1AM(ISYMB,ISYMJ)+NVIR(ISYMB)*(J-1)+B 2227 IADR = ISWTL(ISAIBJ,ISYML)+NT2SQ(ISAIBJ)*(L-1) 2228 * + IT2SQ(ISYMAI,ISYMBJ)+NT1AM(ISYMAI)*(NBJ-1) 2229 * + 1 2230 2231 CALL PUTWA2(LUFILE,FNFILE,T3(KOFF1),IADR,NT1AM(ISYMAI)) 2232C 2233 END DO 2234 END DO 2235 END DO 2236 END DO 2237C 2238C-------------- 2239C End. 2240C-------------- 2241C 2242 CALL QEXIT('WRITE_T3_DL') 2243C 2244 RETURN 2245 END 2246C 2247C /* Deck sort_t2_ij */ 2248 SUBROUTINE SORT_T2_IJ(XIJ,ISYMA,A,ISYMB,B,T2TP,ISYMT2) 2249C 2250C XIJ^AB(ij) = t2tp(aijb) 2251C 2252 IMPLICIT NONE 2253C 2254#include "priunit.h" 2255#include "ccsdsym.h" 2256#include "ccorb.h" 2257C 2258 INTEGER ISYMA, ISYMB, ISYMT2 2259 INTEGER ISYMAB, ISYMIJ, ISYMI, ISYMJ, ISYMAIJ, ISYMAI 2260 INTEGER KOFF1, KOFF2 2261#if defined (SYS_CRAY) 2262 REAL XIJ(*), T2TP(*) 2263#else 2264 DOUBLE PRECISION XIJ(*), T2TP(*) 2265#endif 2266C 2267 CALL QENTER('SORT_T2_IJ') 2268C 2269 ISYMAB = MULD2H(ISYMA,ISYMB) 2270 ISYMIJ = MULD2H(ISYMT2,ISYMAB) 2271 DO ISYMI = 1,NSYM 2272 ISYMJ = MULD2H(ISYMIJ,ISYMI) 2273 ISYMAIJ = MULD2H(ISYMA,ISYMIJ) 2274 ISYMAI = MULD2H(ISYMA,ISYMI) 2275 DO I = 1,NRHF(ISYMI) 2276 DO J = 1,NRHF(ISYMJ) 2277 KOFF1 = IMATIJ(ISYMI,ISYMJ) 2278 * + NRHF(ISYMI)*(J-1) 2279 * + I 2280C 2281 KOFF2 = IT2SP(ISYMAIJ,ISYMB) 2282 * + NCKI(ISYMAIJ)*(B-1) 2283 * + ICKI(ISYMAI,ISYMJ) 2284 * + NT1AM(ISYMAI)*(J-1) 2285 * + IT1AM(ISYMA,ISYMI) 2286 * + NVIR(ISYMA)*(I-1) 2287 * + A 2288C 2289 XIJ(KOFF1) = T2TP(KOFF2) 2290 END DO 2291 END DO 2292 END DO 2293C 2294 CALL QEXIT('SORT_T2_IJ') 2295C 2296 RETURN 2297 END 2298C /* Deck sort_t2_aji */ 2299 SUBROUTINE SORT_T2_AJI(XAJI,ISYMB,B,T2TP,ISYMT2) 2300C 2301C t2tp(aijb) as I^I(AJI) 2302C 2303 IMPLICIT NONE 2304C 2305#include "priunit.h" 2306#include "ccsdsym.h" 2307#include "ccorb.h" 2308C 2309 INTEGER ISYMI, ISYMT2 2310 INTEGER ISYMAJI, ISYMJ, ISYMAJ, ISYMB, ISYMA, ISYMAIJ, ISYMAI 2311 INTEGER KOFF1, KOFF2 2312#if defined (SYS_CRAY) 2313 REAL XAJI(*), T2TP(*) 2314#else 2315 DOUBLE PRECISION XAJI(*), T2TP(*) 2316#endif 2317C 2318 CALL QENTER('SORT_T2_AJI') 2319C 2320 ISYMAIJ = MULD2H(ISYMT2,ISYMB) 2321 ISYMAJI = MULD2H(ISYMT2,ISYMB) 2322 DO ISYMJ = 1,NSYM 2323 ISYMAI = MULD2H(ISYMAJI,ISYMJ) 2324 DO ISYMA = 1,NSYM 2325 ISYMAJ = MULD2H(ISYMJ,ISYMA) 2326 ISYMI = MULD2H(ISYMAI,ISYMA) 2327 DO J = 1,NRHF(ISYMJ) 2328 DO I = 1,NRHF(ISYMI) 2329 DO A = 1,NVIR(ISYMA) 2330 KOFF1 = ISAIK(ISYMAJ,ISYMI) 2331 * + NT1AM(ISYMAJ)*(I-1) 2332 * + IT1AM(ISYMA,ISYMJ) 2333 * + NVIR(ISYMA)*(J-1) + A 2334C 2335 KOFF2 = IT2SP(ISYMAIJ,ISYMB) 2336 * + NCKI(ISYMAIJ)*(B-1) 2337 * + ICKI(ISYMAI,ISYMJ) 2338 * + NT1AM(ISYMAI)*(J-1) 2339 * + IT1AM(ISYMA,ISYMI) 2340 * + NVIR(ISYMA)*(I-1) 2341 * + A 2342C 2343 XAJI(KOFF1) = T2TP(KOFF2) 2344 END DO 2345 END DO 2346 END DO 2347 END DO 2348 END DO 2349C 2350 CALL QEXIT('SORT_T2_AJI') 2351C 2352 RETURN 2353 END 2354