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 cc_den */ 20 SUBROUTINE CC_DEN(POTNUC,AODEN,ZKDIA,WORK,LWORK,IOPT) 21C 22C Written by Asger Halkier February/March 1996. 23C 24C Version: 1.0 25C 26C Purpose: To calculate various different matrices, 27C depending on the value of IOPT, all constructed 28C from the Coupled Cluster density matrices! 29C 30C Current models: CCD, CCSD 31C 32C For IOPT = 1 the right hand side of the equation for the 33C zero'th order orbital rotation lagrange multiplier response 34C (Zeta-Kappa-0) is returned in the array AODEN! 35C 36C For IOPT = 2 both the one and two electron Coupled Cluster lamda 37C density is set up and used to calculate the ccsd-energy ECCSD! 38C (Used to debug the construction of the density-matrices) 39C 40C For IOPT = 3 both the one and two electron effective Coupled Cluster 41C density is set up and used to calculate the ccsd-energy ECCSD! 42C (Used to debug the construction of the density-matrices) 43C 44C Modified Summer 1998 to incorporate canonical surface, so the 45C diagonal blocks of Kappa-bar-0 is returned in ZKDIA for 46C frozen core calculations. 47C 48#include "implicit.h" 49#include "priunit.h" 50#include "dummy.h" 51#include "maxash.h" 52#include "mxcent.h" 53#include "maxorb.h" 54#include "aovec.h" 55#include "iratdef.h" 56 PARAMETER (ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, TWO = 2.0D0) 57 PARAMETER (FOUR = 4.0D0) 58 DIMENSION INDEXA(MXCORB_CC) 59 DIMENSION AODEN(*), ZKDIA(*), WORK(LWORK) 60#include "ccorb.h" 61#include "ccisao.h" 62#include "r12int.h" 63#include "inftap.h" 64#include "blocks.h" 65#include "ccfield.h" 66#include "ccsdinp.h" 67#include "ccinftap.h" 68#include "ccsdsym.h" 69#include "ccsdio.h" 70#include "distcl.h" 71#include "cbieri.h" 72#include "eritap.h" 73#include "ccfro.h" 74C 75 CHARACTER MODEL*10 76 CHARACTER NAME1*8 77 CHARACTER NAME2*8 78C 79 CALL QENTER('CC_DEN') 80C 81 IF ((FROIMP) .AND. (IOPT .EQ. 1)) THEN 82C 83 NAME1 = 'CCFRO1IN' 84 NAME2 = 'CCFRO2IN' 85C 86 LUN1 = -1 87 LUN2 = -1 88C 89 CALL WOPEN2(LUN1,NAME1,64,0) 90 CALL WOPEN2(LUN2,NAME2,64,0) 91C 92 ENDIF 93C 94 IF (IOPT .LT. 2) THEN 95 CALL HEADER('Constructing right-hand-side for CCSD-kappa-0',-1) 96 ENDIF 97C 98C----------------------------------------- 99C Initialization of timing parameters. 100C----------------------------------------- 101C 102 TIMTOT = ZERO 103 TIMTOT = SECOND() 104 TIMDEN = ZERO 105 TIMRES = ZERO 106 TIRDAO = ZERO 107 TIMHE2 = ZERO 108 TIMONE = ZERO 109 TIMONE = SECOND() 110C 111C---------------------------------------------------- 112C Both zeta- and t-vectors are totally symmetric. 113C---------------------------------------------------- 114C 115 ISYMTR = 1 116 ISYMOP = 1 117C 118C----------------------------------- 119C Initial work space allocation. 120C----------------------------------- 121C 122 KZ2AM = 1 123 KT2AM = KZ2AM + NT2SQ(1) 124 KT2AMT = KT2AM + NT2AMX 125 KLAMDP = KT2AMT + NT2AMX 126 KLAMDH = KLAMDP + NLAMDT 127 KT1AM = KLAMDH + NLAMDT 128 KZ1AM = KT1AM + NT1AMX 129 KEND1 = KZ1AM + NT1AMX 130 LWRK1 = LWORK - KEND1 131C 132 IF (LWRK1 .LT. 0) THEN 133 WRITE(LUPRI,*) 'CC_DEN: Available:', LWORK, 'Needed:', KEND1 134 CALL QUIT('Insufficient memory for initial allocation '// 135 & 'in CC_DEN') 136 ENDIF 137C 138C---------------------------------------- 139C Read zero'th order zeta amplitudes. 140C---------------------------------------- 141C 142 IOPTRW = 3 143 CALL CC_RDRSP('L0',0,1,IOPTRW,MODEL,WORK(KZ1AM),WORK(KZ2AM)) 144C 145C-------------------------------- 146C Square up zeta2 amplitudes. 147C-------------------------------- 148C 149 CALL DCOPY(NT2AMX,WORK(KZ2AM),1,WORK(KT2AM),1) 150 CALL CC_T2SQ(WORK(KT2AM),WORK(KZ2AM),1) 151C 152C------------------------------------------- 153C Read zero'th order cluster amplitudes. 154C------------------------------------------- 155C 156 IOPTRW = 3 157 CALL CC_RDRSP('R0',0,1,IOPTRW,MODEL,WORK(KT1AM),WORK(KT2AM)) 158C 159C------------------------------------------------ 160C Zero out single vectors in CCD-calculation. 161C------------------------------------------------ 162C 163 IF (CCD) THEN 164 CALL DZERO(WORK(KT1AM),NT1AMX) 165 CALL DZERO(WORK(KZ1AM),NT1AMX) 166 ENDIF 167C 168C---------------------------------- 169C Calculate the lamda matrices. 170C---------------------------------- 171C 172 CALL LAMMAT(WORK(KLAMDP),WORK(KLAMDH),WORK(KT1AM),WORK(KEND1), 173 * LWRK1) 174C 175C--------------------------------------- 176C Set up 2C-E of cluster amplitudes. 177C--------------------------------------- 178C 179 ISYOPE = 1 180C 181 CALL DCOPY(NT2AMX,WORK(KT2AM),1,WORK(KT2AMT),1) 182 IOPTTCME = 1 183 CALL CCSD_TCMEPK(WORK(KT2AMT),1.0D0,ISYOPE,IOPTTCME) 184C 185C------------------------------- 186C Work space allocation one. 187C Note that D(ai) = ZETA(ai) 188C and both D(ia) and h(ia) 189C are stored transposed! 190C------------------------------- 191C 192 KONEAI = KZ1AM 193 KONEAB = KONEAI + NT1AMX 194 KONEIJ = KONEAB + NMATAB(1) 195 KONEIA = KONEIJ + NMATIJ(1) 196 KXMAT = KONEIA + NT1AMX 197 KYMAT = KXMAT + NMATIJ(1) 198 KMINT = KYMAT + NMATAB(1) 199 KONINT = KMINT + N3ORHF(1) 200 KMIRES = KONINT + N2BST(ISYMOP) 201 IF (IOPT .GT. 1) THEN 202 IF (IOPT .EQ. 2) THEN 203 KEND1 = KMIRES + N3ORHF(1) 204 LWRK1 = LWORK - KEND1 205 ELSE IF (IOPT .EQ. 3) THEN 206 LENBAR = 2*NT1AMX + NMATIJ(1) + NMATAB(1) 207 * + 2*NCOFRO(1) + 2*NT1FRO(1) 208 KKABAR = KMIRES + N3ORHF(1) 209 KDHFAO = KKABAR + LENBAR 210 KKABAO = KDHFAO + N2BST(1) 211 KCMO = KKABAO + N2BST(1) 212 KEND1 = KCMO + NLAMDS 213 LWRK1 = LWORK - KEND1 214 ENDIF 215 ELSE IF (IOPT .LT. 2) THEN 216 KINTAI = KMIRES + N3ORHF(1) 217 KINTAB = KINTAI + NT1AMX 218 KINTIJ = KINTAB + NMATAB(1) 219 KINTIA = KINTIJ + NMATIJ(1) 220 KINABT = KINTIA + NT1AMX 221 KINIJT = KINABT + NMATAB(1) 222 KD1ABT = KINIJT + NMATIJ(1) 223 KD1IJT = KD1ABT + NMATAB(1) 224 KEND1 = KD1IJT + NMATIJ(1) 225 LWRK1 = LWORK - KEND1 226 IF (FROIMP) THEN 227 KFROII = KEND1 228 KFROIJ = KFROII + NFROFR(1) 229 KFROJI = KFROIJ + NCOFRO(1) 230 KFROAI = KFROJI + NCOFRO(1) 231 KFROIA = KFROAI + NT1FRO(1) 232 KFD1II = KFROIA + NT1FRO(1) 233 KEND1 = KFD1II + NFROFR(1) 234 LWRK1 = LWORK - KEND1 235 ENDIF 236 ENDIF 237C 238 IF (FROIMP) THEN 239 KCMOF = KEND1 240 KEND1 = KCMOF + NLAMDS 241 LWRK1 = LWORK - KEND1 242 ENDIF 243C 244 IF (LWRK1 .LT. 0) THEN 245 WRITE(LUPRI,*) 'CC_DEN: Available:', LWORK, 'Needed:', KEND1 246 CALL QUIT('Insufficient memory for allocation 1 CC_DEN') 247 ENDIF 248C 249 IF (FROIMP) THEN 250C 251C------------------------------------------- 252C Get the FULL MO coefficient matrix. 253C------------------------------------------- 254C 255 CALL CMO_ALL(WORK(KCMOF),WORK(KEND1),LWRK1) 256C 257 ENDIF 258C 259C------------------------------------------------------ 260C Initialize remaining one electron density arrays. 261C------------------------------------------------------ 262C 263 CALL DZERO(WORK(KONEAB),NMATAB(1)) 264 CALL DZERO(WORK(KONEIJ),NMATIJ(1)) 265 CALL DZERO(WORK(KONEIA),NT1AMX) 266C 267C-------------------------------------------------------- 268C Calculate X-intermediate of zeta- and t-amplitudes. 269C-------------------------------------------------------- 270C 271 CALL CC_XI(WORK(KXMAT),WORK(KZ2AM),ISYMTR,WORK(KT2AM),ISYMOP, 272 * WORK(KEND1),LWRK1) 273C 274C-------------------------------------------------------- 275C Calculate Y-intermediate of zeta- and t-amplitudes. 276C-------------------------------------------------------- 277C 278 CALL CC_YI(WORK(KYMAT),WORK(KZ2AM),ISYMTR,WORK(KT2AM),ISYMOP, 279 * WORK(KEND1),LWRK1) 280C 281C-------------------------------------------------------------- 282C Construct three remaining blocks of one electron density. 283C-------------------------------------------------------------- 284C 285 CALL DCOPY(NMATAB(1),WORK(KYMAT),1,WORK(KONEAB),1) 286 CALL CC_EITR(WORK(KONEAB),WORK(KONEIJ),WORK(KEND1),LWRK1,1) 287 CALL DIJGEN(WORK(KONEIJ),WORK(KXMAT)) 288 CALL DIAGEN(WORK(KONEIA),WORK(KT2AMT),WORK(KONEAI)) 289C 290C--------------------------------- 291C Read one-electron integrals. 292C--------------------------------- 293C 294 CALL CCRHS_ONEAO(WORK(KONINT),WORK(KEND1),LWRK1) 295C 296 IF (IOPT .LT. 2) THEN 297C 298C--------------------------------------------------------- 299C Ove 24-20-1996 300C Read one-electron integrals into Fock-matrix for 301C finite field. 302C--------------------------------------------------------- 303C 304 DO 13 IF = 1, NFIELD 305c DTIME = SECOND() 306 FF = EFIELD(IF) 307 CALL CC_ONEP(WORK(KONINT),WORK(KEND1),LWRK1,FF,1,LFIELD(IF)) 308c DTIME = DTIME - SECOND() 309c TIMFCK = TIMFCK + DTIME 310 13 CONTINUE 311C 312C-------------------------------------------------- 313C Transform one electron integrals to MO-basis. 314C-------------------------------------------------- 315C 316 ISYM = 1 317 CALL CCDINTMO(WORK(KINTIJ),WORK(KINTIA),WORK(KINTAB), 318 * WORK(KINTAI),WORK(KONINT),WORK(KLAMDP), 319 * WORK(KLAMDH),WORK(KEND1),LWRK1,ISYM) 320C 321 IF (IPRINT .GT. 50) THEN 322C 323 CALL AROUND('One electron integrals in MO-basis') 324C 325 DO 10 ISYM = 1,NSYM 326C 327 WRITE(LUPRI,333) 'Symmetry block number:', ISYM 328 333 FORMAT(/3X,A22,2X,I1) 329C 330 KOFFIJ = KINTIJ + IMATIJ(ISYM,ISYM) 331 KOFFAI = KINTAI + IT1AM(ISYM,ISYM) 332 KOFFIA = KINTIA + IT1AM(ISYM,ISYM) 333 KOFFAB = KINTAB + IMATAB(ISYM,ISYM) 334C 335 CALL AROUND('occ.occ part') 336 CALL OUTPUT(WORK(KOFFIJ),1,NRHF(ISYM),1,NRHF(ISYM), 337 * NRHF(ISYM),NRHF(ISYM),1,LUPRI) 338C 339 CALL AROUND('vir.occ part') 340 CALL OUTPUT(WORK(KOFFAI),1,NVIR(ISYM),1,NRHF(ISYM), 341 * NVIR(ISYM),NRHF(ISYM),1,LUPRI) 342C 343 CALL AROUND('occ.vir part') 344 CALL OUTPUT(WORK(KOFFIA),1,NVIR(ISYM),1,NRHF(ISYM), 345 * NVIR(ISYM),NRHF(ISYM),1,LUPRI) 346C 347 CALL AROUND('vir.vir part') 348 CALL OUTPUT(WORK(KOFFAB),1,NVIR(ISYM),1,NVIR(ISYM), 349 * NVIR(ISYM),NVIR(ISYM),1,LUPRI) 350C 351 10 CONTINUE 352C 353 HIJNO = DDOT(NMATIJ(1),WORK(KINTIJ),1,WORK(KINTIJ),1) 354 HAINO = DDOT(NT1AMX,WORK(KINTAI),1,WORK(KINTAI),1) 355 HIANO = DDOT(NT1AMX,WORK(KINTIA),1,WORK(KINTIA),1) 356 HABNO = DDOT(NMATAB(1),WORK(KINTAB),1,WORK(KINTAB),1) 357C 358 WRITE(LUPRI,*) ' ' 359 WRITE(LUPRI,*) 'Norm of one electron integrals in MO-basis:' 360 WRITE(LUPRI,*) HIJNO, HAINO, HIANO, HABNO 361C 362 ENDIF 363C 364 IF (FROIMP) THEN 365C 366 ISYM = 1 367 CALL CCIFROMO(WORK(KFROIJ),WORK(KFROJI),WORK(KFROAI), 368 * WORK(KFROIA),WORK(KFROII),WORK(KONINT), 369 * WORK(KLAMDP),WORK(KLAMDH),WORK(KCMOF), 370 * WORK(KEND1),LWRK1,ISYM) 371C 372 CALL CCFD1II(WORK(KFD1II)) 373C 374 ENDIF 375C 376C-------------------------------------------------- 377C Set up transposed integrals and densities. 378C-------------------------------------------------- 379C 380 CALL DCOPY(NMATAB(1),WORK(KINTAB),1,WORK(KINABT),1) 381 CALL DCOPY(NMATIJ(1),WORK(KINTIJ),1,WORK(KINIJT),1) 382 CALL DCOPY(NMATAB(1),WORK(KONEAB),1,WORK(KD1ABT),1) 383 CALL DCOPY(NMATIJ(1),WORK(KONEIJ),1,WORK(KD1IJT),1) 384C 385 CALL CC_EITR(WORK(KINABT),WORK(KINIJT),WORK(KEND1),LWRK1,1) 386 CALL CC_EITR(WORK(KD1ABT),WORK(KD1IJT),WORK(KEND1),LWRK1,1) 387C 388C------------------------------------------------------------ 389C Calculate one electron contribution to Zeta-kappa-0. 390C------------------------------------------------------------ 391C 392 CALL DZERO(AODEN,NT1AM(1)) 393C 394 ISYM = 1 395 CALL CCDENZK0(AODEN,WORK(KINTIJ),WORK(KINTAI),WORK(KINTIA), 396 * WORK(KINTAB),WORK(KONEIJ),WORK(KONEAI), 397 * WORK(KONEIA),WORK(KONEAB),WORK(KINIJT), 398 * WORK(KINABT),WORK(KD1IJT),WORK(KD1ABT), 399 * WORK(KT1AM),WORK(KEND1),LWRK1,ISYM) 400C 401 IF (FROIMP) THEN 402C 403C-------------------------------------------------------- 404C Calculate one-electron contribution to right- 405C hand-side of correlated-frozen multipliers. 406C-------------------------------------------------------- 407C 408 ISYM = 1 409 ICON = 1 410 KOFRES = 2*NT1AMX + NMATIJ(1) + NMATAB(1) + 2*NT1FRO(1) + 1 411 CALL CC_ETIJF(ZKDIA(KOFRES),WORK(KONEIJ),WORK(KONEAB), 412 * WORK(KONEAI),WORK(KONEIA),WORK(KD1IJT), 413 * WORK(KFD1II),SK1,SK2,SK3,SK4,WORK(KINTIJ), 414 * WORK(KINTAI),WORK(KINTIA),WORK(KINIJT), 415 * WORK(KINABT),WORK(KFROIJ),WORK(KFROJI), 416 * WORK(KFROAI),WORK(KFROIA),WORK(KFROII), 417 * WORK(KT1AM),WORK(KEND1),LWRK1,ISYM,ICON) 418C 419C-------------------------------------------------------- 420C Calculate one-electron contribution to right- 421C hand-side of virtual-frozen multipliers. 422C-------------------------------------------------------- 423C 424 ISYM = 1 425 ICON = 1 426 KOFRES = 2*NT1AMX + NMATIJ(1) + NMATAB(1) + 1 427 CALL CC_ETAIF(ZKDIA(KOFRES),WORK(KONEAB),WORK(KONEAI), 428 * WORK(KONEIA),WORK(KD1IJT),WORK(KD1ABT), 429 * WORK(KFD1II),SK1,SK2,SK3,SK4,WORK(KINTIJ), 430 * WORK(KINTAB),WORK(KINTAI),WORK(KINTIA), 431 * WORK(KINABT),WORK(KFROIJ),WORK(KFROJI), 432 * WORK(KFROAI),WORK(KFROIA),WORK(KFROII), 433 * WORK(KT1AM),WORK(KEND1),LWRK1,ISYM,ICON) 434C 435C-------------------------------------------------------- 436C Calculate one-electron contribution to right- 437C hand-side of frozen-frozen multipliers. 438C-------------------------------------------------------- 439C 440c ISYM = 1 441c ICON = 1 442c KOFRES = 2*NT1AMX + NMATIJ(1) + NMATAB(1) + 2*NT1FRO(1) 443c * + 2*NCOFRO(1) + 1 444c CALL CC_ETIFJF(ZKDIA(KOFRES),WORK(KFD1II),SK1,SK2,SK3,SK4, 445c * WORK(KFROIJ),WORK(KFROJI),WORK(KFROAI), 446c * WORK(KFROIA),WORK(KFROII),ISYM,ICON) 447C 448C---------------------------------------------------- 449C Calculate one-electron contribution to 450C right-hand-side for diagonal multipliers. 451C---------------------------------------------------- 452C 453c ISYM = 1 454c CALL CCDIAZK0(ZKDIA,WORK(KINTIJ),WORK(KINTAI),WORK(KINTIA), 455c * WORK(KINTAB),WORK(KONEIJ),WORK(KONEAI), 456c * WORK(KONEIA),WORK(KONEAB),WORK(KINIJT), 457c * WORK(KINABT),WORK(KD1IJT),WORK(KD1ABT), 458c * WORK(KT1AM),WORK(KEND1),LWRK1,ISYM) 459C 460 ENDIF 461 ENDIF 462C 463 IF (IOPT .GT. 1) THEN 464C 465C------------------------------------------------- 466C Read orbital relaxation vector from disc. 467C------------------------------------------------- 468C 469 IF (IOPT .EQ. 3) THEN 470C 471 CALL DZERO(WORK(KKABAR),LENBAR) 472C 473 LUBAR0 = -978 474 CALL GPOPEN(LUBAR0,'CCKABAR0','UNKNOWN',' ','UNFORMATTED', 475 * IDUMMY,.FALSE.) 476 REWIND(LUBAR0) 477 READ(LUBAR0) (WORK(KKABAR+I-1), I = 1,LENBAR) 478 CALL GPCLOSE(LUBAR0,'KEEP') 479C 480C---------------------------------------------------------------- 481C Read MO-coefficients from interface file and reorder. 482C---------------------------------------------------------------- 483C 484 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY, 485 & .FALSE.) 486 REWIND LUSIFC 487 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 488 READ (LUSIFC) 489 READ (LUSIFC) 490 READ (LUSIFC) (WORK(KCMO+I-1), I=1,NLAMDS) 491 CALL GPCLOSE(LUSIFC,'KEEP') 492C 493 CALL CMO_REORDER(WORK(KCMO),WORK(KEND1),LWRK1) 494C 495C-------------------------------------------------------------------- 496C Calculate ao-transformed zeta-kappa-bar-0 and HF density. 497C-------------------------------------------------------------------- 498C 499 KOFDIJ = KKABAR 500 KOFDAB = KOFDIJ + NMATIJ(1) 501 KOFDAI = KOFDAB + NMATAB(1) 502 KOFDIA = KOFDAI + NT1AMX 503C 504 ISDEN = 1 505 CALL DZERO(WORK(KKABAO),N2BST(1)) 506 CALL CC_DENAO(WORK(KKABAO),ISDEN,WORK(KOFDAI),WORK(KOFDAB), 507 * WORK(KOFDIJ),WORK(KOFDIA),ISDEN, 508 * WORK(KCMO),1,WORK(KCMO),1,WORK(KEND1),LWRK1) 509C 510 CALL CCS_D1AO(WORK(KDHFAO),WORK(KEND1),LWRK1) 511 IF (FROIMP .OR. FROEXP) THEN 512 MODEL = 'DUMMY' 513 CALL CC_FCD1AO(WORK(KDHFAO),WORK(KEND1),LWRK1,MODEL) 514 ENDIF 515C 516C--------------------------------------------------------------- 517C Add orbital relaxation for effective density matrix. 518C--------------------------------------------------------------- 519C 520 CALL DCOPY(N2BST(1),WORK(KKABAO),1,AODEN,1) 521C 522 ENDIF 523C 524C----------------------------------------------------------- 525C Backtransform the one electron density to AO-basis. 526C----------------------------------------------------------- 527C 528 IF (IOPT .EQ. 2) CALL DZERO(AODEN,N2BST(1)) 529C 530 ISDEN = 1 531 CALL CC_DENAO(AODEN,ISDEN,WORK(KONEAI),WORK(KONEAB), 532 * WORK(KONEIJ),WORK(KONEIA),ISDEN,WORK(KLAMDP),1, 533 * WORK(KLAMDH),1,WORK(KEND1),LWRK1) 534C 535C------------------------------------------------------ 536C Add frozen core contributions to AO densities. 537C------------------------------------------------------ 538C 539 IF (FROIMP) THEN 540 IF (IOPT .EQ. 3) THEN 541C 542 KOFFAI = KKABAR + NMATIJ(1) + NMATAB(1) + 2*NT1AMX 543 KOFFIA = KOFFAI + NT1FRO(1) 544 KOFFIJ = KOFFIA + NT1FRO(1) 545 KOFFJI = KOFFIJ + NCOFRO(1) 546C 547 ISDEN = 1 548 ICON = 1 549 CALL CC_D1FCB(AODEN,WORK(KOFFIJ),WORK(KOFFJI), 550 * WORK(KOFFAI),WORK(KOFFIA),WORK(KEND1), 551 * LWRK1,ISDEN,ICON) 552C 553 ISDEN = 1 554 ICON = 2 555 CALL CC_D1FCB(WORK(KKABAO),WORK(KOFFIJ),WORK(KOFFJI), 556 * WORK(KOFFAI),WORK(KOFFIA),WORK(KEND1), 557 * LWRK1,ISDEN,ICON) 558C 559 ELSE 560C 561 MODEL = 'DUMMY' 562 CALL CC_FCD1AO(AODEN,WORK(KEND1),LWRK1,MODEL) 563C 564 ENDIF 565 ENDIF 566C 567C----------------------------------------------------------- 568C Calculate one electron contribution to ccsd-energy. 569C----------------------------------------------------------- 570C 571 WRITE (LUPRI,*) 'AODEN in CC_DEN for IOPT=',IOPT 572 CALL CC_PRONELAO(AODEN,1) 573C 574 ECCSD1 = DDOT(N2BST(ISYMOP),AODEN,1,WORK(KONINT),1) 575 ECCSD2 = ZERO 576C 577 ENDIF 578C 579C-------------------------------------------------------- 580C Calculate M-intermediate of zeta- and t-amplitudes. 581C-------------------------------------------------------- 582C 583 CALL CC_MI(WORK(KMINT),WORK(KZ2AM),ISYMTR,WORK(KT2AM),ISYMOP, 584 * WORK(KEND1),LWRK1) 585C 586C-------------------------------------------------------- 587C Calculate resorted M-intermediate M(imjk)->M(mkij). 588C-------------------------------------------------------- 589C 590 CALL CC_MIRS(WORK(KMIRES),WORK(KMINT)) 591C 592 TIMONE = SECOND() - TIMONE 593C 594C----------------------------------- 595C Start the loop over integrals. 596C----------------------------------- 597C 598 KENDS2 = KEND1 599 LWRKS2 = LWRK1 600C 601 IF (DIRECT) THEN 602 IF (HERDIR) THEN 603 CALL HERDI1(WORK(KEND1),LWRK1,IPRERI) 604 ELSE 605 KCCFB1 = KEND1 606 KINDXB = KCCFB1 + MXPRIM*MXCONT 607 KEND1 = KINDXB + (8*MXSHEL*MXCONT + 1)/IRAT 608 LWRK1 = LWORK - KEND1 609 CALL ERIDI1(KODCL1,KODCL2,KODBC1,KODBC2,KRDBC1,KRDBC2, 610 * KODPP1,KODPP2,KRDPP1,KRDPP2, 611 * KFREE,LFREE,KEND1,WORK(KCCFB1),WORK(KINDXB), 612 * WORK(KEND1),LWRK1,IPRERI) 613 KEND1 = KFREE 614 LWRK1 = LFREE 615 ENDIF 616 NTOSYM = 1 617 ELSE 618 NTOSYM = NSYM 619 ENDIF 620C 621 KENDSV = KEND1 622 LWRKSV = LWRK1 623C 624 ICDEL1 = 0 625 DO 100 ISYMD1 = 1,NTOSYM 626C 627 IF (DIRECT) THEN 628 IF (HERDIR) THEN 629 NTOT = MAXSHL 630 ELSE 631 NTOT = MXCALL 632 ENDIF 633 ELSE 634 NTOT = NBAS(ISYMD1) 635 ENDIF 636C 637 DO 110 ILLL = 1,NTOT 638C 639C--------------------------------------------- 640C If direct calculate the integrals. 641C--------------------------------------------- 642C 643 IF (DIRECT) THEN 644C 645 KEND1 = KENDSV 646 LWRK1 = LWRKSV 647C 648 DTIME = SECOND() 649 IF (HERDIR) THEN 650 CALL HERDI2(WORK(KEND1),LWRK1,INDEXA,ILLL,NUMDIS, 651 & IPRINT) 652 ELSE 653 CALL ERIDI2(ILLL,INDEXA,NUMDIS,0,0, 654 * WORK(KODCL1),WORK(KODCL2), 655 * WORK(KODBC1),WORK(KODBC2), 656 * WORK(KRDBC1),WORK(KRDBC2), 657 * WORK(KODPP1),WORK(KODPP2), 658 * WORK(KRDPP1),WORK(KRDPP2), 659 * WORK(KCCFB1),WORK(KINDXB), 660 * WORK(KEND1), LWRK1,IPRERI) 661 ENDIF 662 DTIME = SECOND() - DTIME 663 TIMHE2 = TIMHE2 + DTIME 664C 665 KRECNR = KEND1 666 KEND1 = KRECNR + (NBUFX(0) - 1)/IRAT + 1 667 LWRK1 = LWORK - KEND1 668 IF (LWRK1 .LT. 0) THEN 669 CALL QUIT('Insufficient core in CC_DEN') 670 END IF 671C 672 ELSE 673 NUMDIS = 1 674 ENDIF 675C 676C----------------------------------------------------- 677C Loop over number of distributions in disk. 678C----------------------------------------------------- 679C 680 DO 120 IDEL2 = 1,NUMDIS 681C 682 IF (DIRECT) THEN 683 IDEL = INDEXA(IDEL2) 684 IF (NOAUXB) THEN 685 IDUM = 1 686 CALL IJKAUX(IDEL,IDUM,IDUM,IDUM) 687 END IF 688 ISYMD = ISAO(IDEL) 689 ELSE 690 IDEL = IBAS(ISYMD1) + ILLL 691 ISYMD = ISYMD1 692 ENDIF 693C 694C---------------------------------------- 695C Work space allocation two. 696C---------------------------------------- 697C 698 ISYDEN = ISYMD 699C 700 KD2IJG = KEND1 701 KD2AIG = KD2IJG + ND2IJG(ISYDEN) 702 KD2IAG = KD2AIG + ND2AIG(ISYDEN) 703 KD2ABG = KD2IAG + ND2AIG(ISYDEN) 704 KEND2 = KD2ABG + ND2ABG(ISYDEN) 705 LWRK2 = LWORK - KEND2 706C 707 IF (LWRK2 .LT. 0) THEN 708 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:',KEND2 709 CALL QUIT('Insufficient space for allocation '// 710 & '2 in CC_DEN') 711 ENDIF 712C 713C------------------------------------------------------- 714C Initialize 4 two electron density arrays. 715C------------------------------------------------------- 716C 717 CALL DZERO(WORK(KD2IJG),ND2IJG(ISYDEN)) 718 CALL DZERO(WORK(KD2AIG),ND2AIG(ISYDEN)) 719 CALL DZERO(WORK(KD2IAG),ND2AIG(ISYDEN)) 720 CALL DZERO(WORK(KD2ABG),ND2ABG(ISYDEN)) 721C 722C------------------------------------------------------------------- 723C Calculate the two electron density d(pq,gamma;delta). 724C------------------------------------------------------------------- 725C 726 AUTIME = SECOND() 727C 728 CALL CC_DEN2(WORK(KD2IJG),WORK(KD2AIG),WORK(KD2IAG), 729 * WORK(KD2ABG),WORK(KZ2AM),WORK(KT2AM), 730 * WORK(KT2AMT),WORK(KMINT),WORK(KXMAT), 731 * WORK(KYMAT),WORK(KONEAB),WORK(KONEAI), 732 * WORK(KONEIA),WORK(KMIRES),WORK(KLAMDH),1, 733 * WORK(KLAMDP),1,WORK(KEND2),LWRK2,IDEL,ISYMD) 734C 735 AUTIME = SECOND() - AUTIME 736 TIMDEN = TIMDEN + AUTIME 737C 738C------------------------------------------ 739C Work space allocation three. 740C------------------------------------------ 741C 742 ISYDIS = MULD2H(ISYMD,ISYMOP) 743C 744 KXINT = KEND2 745 KEND3 = KXINT + NDISAO(ISYDIS) 746 LWRK3 = LWORK - KEND3 747C 748 IF (LWRK3 .LT. 0) THEN 749 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:',KEND3 750 CALL QUIT('Insufficient space for allocation '// 751 & '3 in CC_DEN') 752 ENDIF 753C 754C-------------------------------------------- 755C Read AO integral distribution. 756C-------------------------------------------- 757C 758 AUTIME = SECOND() 759 CALL CCRDAO(WORK(KXINT),IDEL,IDEL2,WORK(KEND3),LWRK3, 760 * WORK(KRECNR),DIRECT) 761 AUTIME = SECOND() - AUTIME 762 TIRDAO = TIRDAO + AUTIME 763C 764C---------------------------------------------------------------------- 765C Calculate integral intermediates needed for frozen core. 766C---------------------------------------------------------------------- 767C 768 IF ((FROIMP) .AND. (IOPT .EQ. 1)) THEN 769C 770 KDSRHF = KEND3 771 KOFOIN = KDSRHF + NDSRHF(ISYMD) 772 KDSFRO = KOFOIN + NOFROO(ISYDIS) 773 KSCRAI = KDSFRO + NDSFRO(ISYDIS) 774 KSCAIF = KSCRAI + NOFROO(ISYDIS) 775 KEND3 = KSCAIF + NCOFRF(ISYDIS) 776 LWRK3 = LWORK - KEND3 777C 778 IF (LWRK3 .LT. 0) THEN 779 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:',KEND3 780 CALL QUIT('Insufficient space for allocation '// 781 & 'in CC_DEN') 782 ENDIF 783C 784 CALL DZERO(WORK(KSCRAI),NOFROO(ISYDIS)) 785 CALL DZERO(WORK(KSCAIF),NCOFRF(ISYDIS)) 786C 787C------------------------------------------------------------------------- 788C Transform one index in the integral batch to correlated. 789C------------------------------------------------------------------------- 790C 791 CALL CCTRBT(WORK(KXINT),WORK(KDSRHF),WORK(KLAMDP), 792 * ISYMOP,WORK(KEND3),LWRK3,ISYDIS) 793C 794C--------------------------------------------------------------------- 795C Transform one index in the integral batch to frozen. 796C--------------------------------------------------------------------- 797C 798 CALL CC_GTOFRO(WORK(KXINT),WORK(KDSFRO),WORK(KCMOF), 799 * WORK(KEND3),LWRK3,ISYDIS) 800C 801C-------------------------------------------------------------- 802C Calculate integral batch (cor fro | cor del). 803C-------------------------------------------------------------- 804C 805 CALL CC_OFROIN(WORK(KDSRHF),WORK(KOFOIN),WORK(KCMOF), 806 * WORK(KEND3),LWRK3,ISYDIS) 807C 808C------------------------------------------------------------------ 809C Calculate terms to ai-part from KOFOIN integrals. 810C------------------------------------------------------------------ 811C 812 CALL CC_FRCOC1(WORK(KSCRAI),WORK(KOFOIN),ISYDIS) 813C 814C---------------------------------------------------------------- 815C Calculate exchange parts with KDSFRO integrals. 816C---------------------------------------------------------------- 817C 818 CALL CC_FRCOMI(WORK(KSCRAI),WORK(KSCAIF), 819 * WORK(KDSFRO),WORK(KCMOF), 820 * WORK(KEND3),LWRK3,ISYMD) 821C 822C---------------------------------------------------- 823C Calculate coulomb part of aI block. 824C---------------------------------------------------- 825C 826 CALL CC_FRCOF1(WORK(KSCAIF),WORK(KDSFRO),WORK(KCMOF), 827 * WORK(KEND3),LWRK3,ISYMD) 828C 829C----------------------------------------------------- 830C Calculate exchange part of aI block. 831C----------------------------------------------------- 832C 833 CALL CC_FRCOF2(WORK(KSCAIF),WORK(KDSRHF),WORK(KCMOF), 834 * WORK(KEND3),LWRK3,ISYMD) 835C 836C---------------------------------------------------------- 837C Dump intermediates to disc for later use. 838C---------------------------------------------------------- 839C 840 CALL CC_FRINDU(WORK(KSCRAI),WORK(KSCAIF),IDEL,ISYMD, 841 & LUN1,LUN2) 842C 843 ENDIF 844C 845C------------------------------------------------------ 846C Start loop over second AO-index (gamma). 847C------------------------------------------------------ 848C 849 AUTIME = SECOND() 850C 851 DO 130 ISYMG = 1,NSYM 852C 853 ISYMPQ = MULD2H(ISYMG,ISYDEN) 854C 855 DO 140 G = 1,NBAS(ISYMG) 856C 857 KD2GIJ = KD2IJG + ID2IJG(ISYMPQ,ISYMG) 858 * + NMATIJ(ISYMPQ)*(G - 1) 859 KD2GAI = KD2AIG + ID2AIG(ISYMPQ,ISYMG) 860 * + NT1AM(ISYMPQ)*(G - 1) 861 KD2GAB = KD2ABG + ID2ABG(ISYMPQ,ISYMG) 862 * + NMATAB(ISYMPQ)*(G - 1) 863 KD2GIA = KD2IAG + ID2AIG(ISYMPQ,ISYMG) 864 * + NT1AM(ISYMPQ)*(G - 1) 865C 866C----------------------------------------------- 867C Work space allocation four. 868C----------------------------------------------- 869C 870 KINTAO = KEND3 871 KD2AOB = KINTAO + N2BST(ISYMPQ) 872 KEND4 = KD2AOB + N2BST(ISYMPQ) 873 LWRK4 = LWORK - KEND4 874C 875 IF (LWRK4 .LT. 0) THEN 876 WRITE(LUPRI,*) 'Available:', LWORK 877 WRITE(LUPRI,*) 'Needed:', KEND4 878 CALL QUIT('Insufficient work space in CC_DEN') 879 ENDIF 880C 881 CALL DZERO(WORK(KD2AOB),N2BST(ISYMPQ)) 882C 883C------------------------------------------------------------- 884C Calculate frozen core contributions to d. 885C------------------------------------------------------------- 886C 887 IF (FROIMP) THEN 888C 889 KFD2IJ = KEND4 890 KFD2JI = KFD2IJ + NCOFRO(ISYMPQ) 891 KFD2AI = KFD2JI + NCOFRO(ISYMPQ) 892 KFD2IA = KFD2AI + NT1FRO(ISYMPQ) 893 KFD2II = KFD2IA + NT1FRO(ISYMPQ) 894 KEND4 = KFD2II + NFROFR(ISYMPQ) 895 LWRK4 = LWORK - KEND4 896C 897 IF (LWRK4 .LT. 0) THEN 898 WRITE (LUPRI,*) 'Available:', LWORK 899 WRITE (LUPRI,*) 'Needed:', KEND4 900 CALL QUIT( 901 * 'Insufficient work space in CC_DEN') 902 ENDIF 903C 904 CALL DZERO(WORK(KFD2IJ),NCOFRO(ISYMPQ)) 905 CALL DZERO(WORK(KFD2JI),NCOFRO(ISYMPQ)) 906 CALL DZERO(WORK(KFD2AI),NT1FRO(ISYMPQ)) 907 CALL DZERO(WORK(KFD2IA),NT1FRO(ISYMPQ)) 908 CALL DZERO(WORK(KFD2II),NFROFR(ISYMPQ)) 909C 910 CALL CC_FD2BL(WORK(KFD2II),WORK(KFD2IJ), 911 * WORK(KFD2JI),WORK(KFD2AI), 912 * WORK(KFD2IA),WORK(KONEIJ), 913 * WORK(KONEAB),WORK(KONEAI), 914 * WORK(KONEIA),WORK(KCMOF), 915 * WORK(KLAMDH),WORK(KLAMDP), 916 * WORK(KEND4),LWRK4,IDEL, 917 * ISYMD,G,ISYMG) 918C 919 CALL CC_FD2AO(WORK(KD2AOB),WORK(KFD2II), 920 * WORK(KFD2IJ),WORK(KFD2JI), 921 * WORK(KFD2AI),WORK(KFD2IA), 922 * WORK(KCMOF),WORK(KLAMDH), 923 * WORK(KLAMDP),WORK(KEND4),LWRK4, 924 * ISYMPQ) 925C 926 CALL CC_D2GAF(WORK(KD2GIJ),WORK(KD2GAB), 927 * WORK(KD2GAI),WORK(KD2GIA), 928 * WORK(KONEIJ),WORK(KONEAB), 929 * WORK(KONEAI),WORK(KONEIA), 930 * WORK(KCMOF),IDEL,ISYMD,G,ISYMG) 931C 932 ENDIF 933C 934C------------------------------------------------------- 935C Square up AO-integral distribution. 936C------------------------------------------------------- 937C 938 KOFFIN = KXINT + IDSAOG(ISYMG,ISYDIS) 939 * + NNBST(ISYMPQ)*(G - 1) 940C 941 CALL CCSD_SYMSQ(WORK(KOFFIN),ISYMPQ, 942 * WORK(KINTAO)) 943C 944C------------------------------------------------------------ 945C Backtransform density fully to AO basis. 946C------------------------------------------------------------ 947C 948 IF (IOPT .GE. 2) THEN 949C 950 CALL CC_DENAO(WORK(KD2AOB),ISYMPQ, 951 * WORK(KD2GAI),WORK(KD2GAB), 952 * WORK(KD2GIJ),WORK(KD2GIA),ISYMPQ, 953 * WORK(KLAMDP),1,WORK(KLAMDH),1, 954 * WORK(KEND4),LWRK4) 955C 956C--------------------------------------------------------------------- 957C Add relaxation terms to set up effective density. 958C--------------------------------------------------------------------- 959C 960 IF (IOPT .EQ. 3) THEN 961C 962 ICON = 1 963 CALL CC_D2EFF(WORK(KD2AOB),G,ISYMG,IDEL, 964 * ISYMD,WORK(KKABAO), 965 * WORK(KDHFAO),ICON) 966C 967 ENDIF 968C 969C---------------------------------------------------------------------- 970C Calcalate the 2 e- density contribution to E-ccsd. 971C---------------------------------------------------------------------- 972C 973 ECCSD2 = ECCSD2 + HALF*DDOT(N2BST(ISYMPQ), 974 * WORK(KD2AOB),1,WORK(KINTAO),1) 975C 976 ELSE 977C 978C----------------------------------------------- 979C Work space allocation five. 980C----------------------------------------------- 981C 982 KIJINT = KEND4 983 KAIINT = KIJINT + NMATIJ(ISYMPQ) 984 KIAINT = KAIINT + NT1AM(ISYMPQ) 985 KABINT = KIAINT + NT1AM(ISYMPQ) 986 KABTIN = KABINT + NMATAB(ISYMPQ) 987 KIJTIN = KABTIN + NMATAB(ISYMPQ) 988 KD2TAB = KIJTIN + NMATIJ(ISYMPQ) 989 KD2TIJ = KD2TAB + NMATAB(ISYMPQ) 990 KEND5 = KD2TIJ + NMATIJ(ISYMPQ) 991 LWRK5 = LWORK - KEND5 992 IF (FROIMP) THEN 993 KIIFRO = KEND5 994 KIJFRO = KIIFRO + NFROFR(ISYMPQ) 995 KJIFRO = KIJFRO + NCOFRO(ISYMPQ) 996 KAIFRO = KJIFRO + NCOFRO(ISYMPQ) 997 KIAFRO = KAIFRO + NT1FRO(ISYMPQ) 998 KEND5 = KIAFRO + NT1FRO(ISYMPQ) 999 LWRK5 = LWORK - KEND5 1000 ENDIF 1001C 1002 IF (LWRK5 .LT. 0) THEN 1003 WRITE(LUPRI,*) 'Available:', LWORK 1004 WRITE(LUPRI,*) 'Needed:', KEND5 1005 CALL QUIT('Insufficient work space '// 1006 & 'in CC_DEN') 1007 ENDIF 1008C 1009C---------------------------------------------------------------- 1010C Transform 2 integral indices to MO-basis. 1011C---------------------------------------------------------------- 1012C 1013 ISYM = ISYMPQ 1014 CALL CCDINTMO(WORK(KIJINT),WORK(KIAINT), 1015 * WORK(KABINT),WORK(KAIINT), 1016 * WORK(KINTAO),WORK(KLAMDP), 1017 * WORK(KLAMDH),WORK(KEND5), 1018 * LWRK5,ISYM) 1019C 1020 IF (FROIMP) THEN 1021C 1022 ISYM = ISYMPQ 1023 CALL CCIFROMO(WORK(KIJFRO),WORK(KJIFRO), 1024 * WORK(KAIFRO),WORK(KIAFRO), 1025 * WORK(KIIFRO),WORK(KINTAO), 1026 * WORK(KLAMDP),WORK(KLAMDH), 1027 * WORK(KCMOF),WORK(KEND5), 1028 * LWRK5,ISYM) 1029C 1030 ENDIF 1031C 1032C----------------------------------------------------------------- 1033C Set up transposed integrals and densities. 1034C----------------------------------------------------------------- 1035C 1036 CALL DCOPY(NMATAB(ISYMPQ),WORK(KABINT),1, 1037 * WORK(KABTIN),1) 1038 CALL DCOPY(NMATIJ(ISYMPQ),WORK(KIJINT),1, 1039 * WORK(KIJTIN),1) 1040 CALL DCOPY(NMATAB(ISYMPQ),WORK(KD2GAB),1, 1041 * WORK(KD2TAB),1) 1042 CALL DCOPY(NMATIJ(ISYMPQ),WORK(KD2GIJ),1, 1043 * WORK(KD2TIJ),1) 1044C 1045 CALL CC_EITR(WORK(KABTIN),WORK(KIJTIN), 1046 * WORK(KEND5),LWRK5,ISYMPQ) 1047 CALL CC_EITR(WORK(KD2TAB),WORK(KD2TIJ), 1048 * WORK(KEND5),LWRK5,ISYMPQ) 1049C 1050C------------------------------------------------------------------- 1051C Calculate 2 e- contribution to Zeta-Kappa-0. 1052C------------------------------------------------------------------- 1053C 1054 ISYM = ISYMPQ 1055 CALL CCDENZK0(AODEN,WORK(KIJINT),WORK(KAIINT), 1056 * WORK(KIAINT),WORK(KABINT), 1057 * WORK(KD2GIJ),WORK(KD2GAI), 1058 * WORK(KD2GIA),WORK(KD2GAB), 1059 * WORK(KIJTIN),WORK(KABTIN), 1060 * WORK(KD2TIJ),WORK(KD2TAB), 1061 * WORK(KT1AM),WORK(KEND5),LWRK5, 1062 * ISYM) 1063C 1064 IF (FROIMP) THEN 1065C 1066 ISYM = ISYMPQ 1067C 1068 CALL CCFRETAI(AODEN,WORK(KIJFRO), 1069 * WORK(KJIFRO),WORK(KAIFRO), 1070 * WORK(KIAFRO),WORK(KFD2IJ), 1071 * WORK(KFD2JI),WORK(KFD2AI), 1072 * WORK(KFD2IA),WORK(KT1AM), 1073 * WORK(KEND5),LWRK5,ISYM) 1074C 1075C----------------------------------------------------------------------- 1076C Calculate two-electron contribution to right- 1077C hand-side of correlated-frozen multipliers. 1078C----------------------------------------------------------------------- 1079C 1080 ICON = 2 1081 KOFRES = 2*NT1AMX + NMATIJ(1) + NMATAB(1) 1082 * + 2*NT1FRO(1) + 1 1083C 1084 CALL CC_ETIJF(ZKDIA(KOFRES),WORK(KD2GIJ), 1085 * WORK(KD2GAB),WORK(KD2GAI), 1086 * WORK(KD2GIA),WORK(KD2TIJ), 1087 * WORK(KFD2II),WORK(KFD2IJ), 1088 * WORK(KFD2JI),WORK(KFD2AI), 1089 * WORK(KFD2IA),WORK(KIJINT), 1090 * WORK(KAIINT),WORK(KIAINT), 1091 * WORK(KIJTIN),WORK(KABTIN), 1092 * WORK(KIJFRO),WORK(KJIFRO), 1093 * WORK(KAIFRO),WORK(KIAFRO), 1094 * WORK(KIIFRO),WORK(KT1AM), 1095 * WORK(KEND5),LWRK5,ISYM,ICON) 1096C 1097C----------------------------------------------------------------------- 1098C Calculate two-electron contribution to right- 1099C hand-side of virtual-frozen multipliers. 1100C----------------------------------------------------------------------- 1101C 1102 ICON = 2 1103 KOFRE = 2*NT1AMX + NMATIJ(1) + NMATAB(1) + 1 1104C 1105 CALL CC_ETAIF(ZKDIA(KOFRE),WORK(KD2GAB), 1106 * WORK(KD2GAI),WORK(KD2GIA), 1107 * WORK(KD2TIJ),WORK(KD2TAB), 1108 * WORK(KFD2II),WORK(KFD2IJ), 1109 * WORK(KFD2JI),WORK(KFD2AI), 1110 * WORK(KFD2IA),WORK(KIJINT), 1111 * WORK(KABINT),WORK(KAIINT), 1112 * WORK(KIAINT),WORK(KABTIN), 1113 * WORK(KIJFRO),WORK(KJIFRO), 1114 * WORK(KAIFRO),WORK(KIAFRO), 1115 * WORK(KIIFRO),WORK(KT1AM), 1116 * WORK(KEND5),LWRK5,ISYM,ICON) 1117C 1118C----------------------------------------------------------------------- 1119C Calculate two-electron contribution to right- 1120C hand-side of frozen-frozen multipliers. 1121C----------------------------------------------------------------------- 1122C 1123 ICON = 2 1124 KOFRES = 2*NT1AMX + NMATIJ(1) + NMATAB(1) 1125 * + 2*NT1FRO(1) + 2*NCOFRO(1) + 1 1126C 1127c CALL CC_ETIFJF(ZKDIA(KOFRES),WORK(KFD2II), 1128c * WORK(KFD2IJ),WORK(KFD2JI), 1129c * WORK(KFD2AI),WORK(KFD2IA), 1130c * WORK(KIJFRO),WORK(KJIFRO), 1131c * WORK(KAIFRO),WORK(KIAFRO), 1132c * WORK(KIIFRO),ISYM,ICON) 1133C 1134c CALL CCDIAZK0(ZKDIA,WORK(KIJINT), 1135c * WORK(KAIINT), 1136c * WORK(KIAINT),WORK(KABINT), 1137c * WORK(KD2GIJ),WORK(KD2GAI), 1138c * WORK(KD2GIA),WORK(KD2GAB), 1139c * WORK(KIJTIN),WORK(KABTIN), 1140c * WORK(KD2TIJ),WORK(KD2TAB), 1141c * WORK(KT1AM),WORK(KEND5), 1142c * LWRK5,ISYM) 1143C 1144c CALL CCFRETAB(ZKDIA,WORK(KIJFRO), 1145c * WORK(KJIFRO),WORK(KAIFRO), 1146c * WORK(KIAFRO),WORK(KFD2IJ), 1147c * WORK(KFD2JI),WORK(KFD2AI), 1148c * WORK(KFD2IA),WORK(KT1AM), 1149c * WORK(KEND5),LWRK5,ISYM) 1150C 1151c CALL CCFRETIJ(ZKDIA,WORK(KIJFRO), 1152c * WORK(KJIFRO),WORK(KAIFRO), 1153c * WORK(KIAFRO),WORK(KFD2IJ), 1154c * WORK(KFD2JI),WORK(KFD2AI), 1155c * WORK(KFD2IA),WORK(KT1AM), 1156c * WORK(KEND5),LWRK5,ISYM) 1157C 1158 ENDIF 1159 ENDIF 1160C 1161 140 CONTINUE 1162 130 CONTINUE 1163C 1164 AUTIME = SECOND() - AUTIME 1165 TIMRES = TIMRES + AUTIME 1166C 1167 120 CONTINUE 1168 110 CONTINUE 1169 100 CONTINUE 1170C 1171C----------------------- 1172C Regain work space. 1173C----------------------- 1174C 1175 KEND1 = KENDS2 1176 LWRK1 = LWRKS2 1177C 1178C--------------------------------------------------------------- 1179C Add nuclear nuclear repulsion energy and write out E-ccsd. 1180C--------------------------------------------------------------- 1181C 1182 IF (IOPT .GT. 1) THEN 1183C 1184 ECCSD = ECCSD1 + ECCSD2 + POTNUC 1185C 1186 WRITE(LUPRI,*) ' ' 1187 WRITE(LUPRI,*) 'Coupled Cluster energy constructed' 1188 WRITE(LUPRI,*) 'from density matrices:' 1189 IF (CCD) WRITE(LUPRI,*) 'CCD-energy:', ECCSD 1190 IF (CCSD) WRITE(LUPRI,*) 'CCSD-energy:', ECCSD 1191 WRITE(LUPRI,*) 'H1 energy, ECCSD1 = ',ECCSD1 1192 WRITE(LUPRI,*) 'H2 energy, ECCSD2 = ',ECCSD2 1193 WRITE(LUPRI,*) 'Nuc. Pot. energy = ',POTNUC 1194C 1195 ENDIF 1196C 1197C------------------------------------------------ 1198C Write out frozen block of eta-kappa-bar-0. 1199C------------------------------------------------ 1200C 1201 IF ((IPRINT .GT. 9) .AND. (FROIMP) .AND. (IOPT .EQ. 1)) THEN 1202 CALL AROUND('Eta-kappa-bar-0 correlated-frozen block') 1203 KOFRES = NMATIJ(1) + NMATAB(1) + 2*NT1AMX + 2*NT1FRO(1) + 1 1204 ZKAPF1 = DDOT(NCOFRO(1),ZKDIA(KOFRES),1, 1205 * ZKDIA(KOFRES),1) 1206 ZKAPFR = ZKAPF1**0.5 1207 WRITE(LUPRI,*) 'Norm of correlated-frozen block:', ZKAPFR 1208 ENDIF 1209C 1210C-------------------------------------------------- 1211C Calculate the diagonal blocks of kappa-bar-0. 1212C-------------------------------------------------- 1213C 1214c IF ((FROIMP) .AND. (IOPT .EQ. 1)) THEN 1215c CALL CCSD_ZKBLO(ZKDIA,WORK(KEND1),LWRK1) 1216c ENDIF 1217C 1218C----------------------------------------------------------- 1219C Calculate the correlated-frozen blocks of kappa-bar-0. 1220C----------------------------------------------------------- 1221C 1222 IF ((FROIMP) .AND. (IOPT .EQ. 1)) THEN 1223 KOFFIJ = 2*NT1AMX + NMATIJ(1) + NMATAB(1) + 2*NT1FRO(1) + 1 1224 CALL CC_ZKFCB(ZKDIA(KOFFIJ),WORK(KEND1),LWRK1) 1225C 1226C---------------------------------------------------------------- 1227C Calculate correction terms from correlated-frozen block. 1228C---------------------------------------------------------------- 1229C 1230 KOFFAI = 2*NT1AMX + NMATIJ(1) + NMATAB(1) + 1 1231 CALL CC_ETACOR(AODEN,ZKDIA(KOFFAI),ZKDIA(KOFFIJ), 1232 * WORK(KCMOF),LUN1,LUN2,WORK(KEND1),LWRK1) 1233C 1234 ENDIF 1235C 1236C--------------------- 1237C Reorder results. 1238C--------------------- 1239C 1240 KOFFAI = 2*NT1AMX + NMATIJ(1) + NMATAB(1) + 1 1241 CALL CC_ETARE(AODEN,ZKDIA(KOFFAI),WORK(KEND1),LWRK1) 1242C 1243C--------------------------------- 1244C Write out eta-ai and eta-aI. 1245C--------------------------------- 1246C 1247 IF (IPRINT .GT. 20) THEN 1248C 1249 CALL AROUND('Eta-kappa-bar-0-ai vector exiting CC_DEN') 1250C 1251 DO 20 ISYM = 1,NSYM 1252C 1253 WRITE(LUPRI,*) ' ' 1254 WRITE(LUPRI,444) 'Sub-symmetry block number:', ISYM 1255 WRITE(LUPRI,555) '--------------------------' 1256 444 FORMAT(3X,A26,2X,I1) 1257 555 FORMAT(3X,A25) 1258C 1259 KOFF = IALLAI(ISYM,ISYM) + 1 1260 CALL OUTPUT(AODEN(KOFF),1,NVIR(ISYM),1,NRHFS(ISYM), 1261 * NVIR(ISYM),NRHFS(ISYM),1,LUPRI) 1262C 1263 IF ((NVIR(ISYM) .EQ. 0) .OR. (NRHFS(ISYM) .EQ. 0)) THEN 1264 WRITE(LUPRI,*) 'This sub-symmetry is empty' 1265 ENDIF 1266C 1267 20 CONTINUE 1268 ENDIF 1269C 1270 IF (IPRINT .GT. 9) THEN 1271 ETAKAN = DDOT(NALLAI(1),AODEN,1,AODEN,1) 1272 WRITE(LUPRI,*) ' ' 1273 WRITE(LUPRI,*) 'Norm of occupied-virtual block:', ETAKAN 1274 ENDIF 1275C 1276C------------------------------------------- 1277C Write out frozen block of kappa-bar-0. 1278C------------------------------------------- 1279C 1280 IF ((IPRINT .GT. 9) .AND. (FROIMP) .AND. (IOPT .EQ. 1)) THEN 1281 CALL AROUND('Zeta-kappa-0 correlated-frozen block') 1282 KOFRES = 2*NT1AMX + NMATIJ(1) + NMATAB(1) + 2*NT1FRO(1) + 1 1283 ZKAPF1 = DDOT(NCOFRO(1),ZKDIA(KOFRES),1, 1284 * ZKDIA(KOFRES),1) 1285 ZKAPFR = ZKAPF1**0.5 1286 WRITE(LUPRI,*) 'Norm of correlated-frozen block:', ZKAPFR 1287 ENDIF 1288C 1289C------------------------------ 1290C Close intermediate files. 1291C------------------------------ 1292C 1293 IF ((FROIMP) .AND. (IOPT .EQ. 1)) THEN 1294C 1295 CALL WCLOSE2(LUN1,NAME1,'KEEP') 1296 CALL WCLOSE2(LUN2,NAME2,'KEEP') 1297 ENDIF 1298C 1299C----------------------- 1300C Write out timings. 1301C----------------------- 1302C 1303 99 TIMTOT = SECOND() - TIMTOT 1304C 1305 IF (IPRINT .GT. 3) THEN 1306 WRITE (LUPRI,*) ' ' 1307 WRITE (LUPRI,*) 'Requested density dependent '// 1308 & 'quantities calculated' 1309 WRITE (LUPRI,*) 'Total time used in CC_DEN:', TIMTOT 1310 ENDIF 1311 IF (IPRINT .GT. 9) THEN 1312 WRITE (LUPRI,*) 'Time used for setting up 2 e- density:',TIMDEN 1313 WRITE (LUPRI,*) 'Time used for contraction with integrals:', 1314 & TIMRES 1315 WRITE (LUPRI,*) 'Time used for reading 2 e- AO-integrals:', 1316 & TIRDAO 1317 WRITE (LUPRI,*) 'Time used for calculating 2 e- AO-integrals:', 1318 * TIMHE2 1319 WRITE (LUPRI,*) 'Time used for 1 e- density & intermediates:', 1320 * TIMONE 1321 ENDIF 1322C 1323 CALL QEXIT('CC_DEN') 1324 RETURN 1325 END 1326C /* Deck cc_den2 */ 1327 SUBROUTINE CC_DEN2(D2IJG,D2AIG,D2IAG,D2ABG,Z2AM,T2AM,T2AMT,XMINT, 1328 * XMAT,YMAT,DAB,DAI,DIA,XMIRES, 1329 * XLAMDH,ISYMLH,XLAMDP,ISYMLP, 1330 * WORK,LWORK,IDEL,ISYMD) 1331C 1332C Written by Asger Halkier 19/2 - 1996 1333C 1334C Version: 1.0 1335C 1336C Purpose: Directs the calculation of the 2 electron CC density 1337C d(pq,gam;del) for a given del (IDEL). The 4 blocks pq 1338C of the result is returned in D2IJG through D2ABG! 1339C 1340C Modifications for CC2 by A. Halkier and S. Coriani 13/01-2000. 1341C 1342#include "implicit.h" 1343 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 1344 INTEGER ISYMLP, ISYMLH 1345 DIMENSION D2IJG(*), D2AIG(*), D2IAG(*), D2ABG(*), Z2AM(*) 1346 DIMENSION T2AM(*), T2AMT(*), XMINT(*), XMAT(*), YMAT(*) 1347 DIMENSION DAB(*), DAI(*), DIA(*), XMIRES(*), XLAMDH(*) 1348 DIMENSION XLAMDP(*), WORK(LWORK) 1349#include "priunit.h" 1350#include "ccorb.h" 1351#include "ccsdsym.h" 1352#include "ccsdinp.h" 1353C 1354 CALL QENTER('CC_DEN2') 1355 IF (DEBUG) WRITE(LUPRI,*) 'Entered CC_DEN2' 1356C 1357C------------------------------- 1358C set some symmetries: 1359C------------------------------- 1360C 1361 ISYD1 = 1 ! one-particle density 1362 ISYMTR = 1 ! Z2AM 1363 ISYD2H = MULD2H(ISYMD,ISYMLH) ! lambdah backtransformed density 1364 ISYDEN = MULD2H(ISYD2H,ISYMLP) ! lambdah + lambdap transformed 1365 ISYMTI = MULD2H(ISYMLH,ISYMD) 1366C 1367C------------------------------- 1368C Work space allocation one. 1369C------------------------------- 1370C 1371 KT2DEL = 1 1372 KEND1 = KT2DEL + NT2BCD(ISYMTI) 1373 LWRK1 = LWORK - KEND1 1374C 1375 IF (LWRK1 .LT. 0) THEN 1376 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 1377 CALL QUIT('Insufficient space for first allocation in CC_DEN2') 1378 ENDIF 1379C 1380C------------------------------------------------------ 1381C Calculate the delta backtransformed t2-amplitude. 1382C------------------------------------------------------ 1383C 1384 CALL CC_TI(WORK(KT2DEL),ISYMTI,T2AM,ISYMOP,XLAMDH,ISYMLH, 1385 * WORK(KEND1),LWRK1,IDEL,ISYMD) 1386C 1387C--------------------------------------------------- 1388C Calculate the (occ.occ,vir;del) density block. 1389C--------------------------------------------------- 1390C 1391 KD2IJA = KEND1 1392 KEND2 = KD2IJA + NMAIJA(ISYD2H) 1393 LWRK2 = LWORK - KEND2 1394C 1395 IF (LWRK2 .LT. 0) THEN 1396 WRITE (LUPRI,*) 'Available:', LWORK, 'Needed:', KEND2 1397 CALL QUIT( 1398 * 'Insufficient space for second allocation in CC_DEN2') 1399 ENDIF 1400C 1401 CALL DZERO(WORK(KD2IJA),NMAIJA(ISYD2H)) 1402C 1403C------------------------------------------------ 1404C Contributions from projection against <u1|. 1405C Remember that D_ai = zeta_ai. 1406C------------------------------------------------ 1407C 1408 CALL DIJA11(WORK(KD2IJA),ISYD2H,DAI,ISYD1,XLAMDH,ISYMLH, 1409 & WORK(KEND2),LWRK2,IDEL,ISYMD) 1410C 1411C------------------------------------------------ 1412C Contributions from projection against <u2|. 1413C------------------------------------------------ 1414C 1415 IF (CCSD) THEN 1416 CALL DIJA21(WORK(KD2IJA),ISYD2H,DAB,ISYD1,XLAMDH,ISYMLH, 1417 * WORK(KEND2),LWRK2,IDEL,ISYMD) 1418 CALL DIJA22(WORK(KD2IJA),ISYD2H,Z2AM,WORK(KT2DEL),ISYMTI, 1419 * WORK(KEND2),LWRK2) 1420 CALL DIJA23(WORK(KD2IJA),ISYD2H,Z2AM,WORK(KT2DEL),ISYMTI, 1421 * WORK(KEND2),LWRK2) 1422 ENDIF 1423C 1424C----------------------------------------------------------------- 1425C Backtransform third virtual index to AO and store in result. 1426C----------------------------------------------------------------- 1427C 1428 ICON = 4 1429 CALL CCD2_PQAO(D2IJG,ISYDEN,WORK(KD2IJA),ISYD2H, 1430 & XLAMDP,ISYMLP,ICON) 1431C 1432 IF (CCSD) THEN 1433C 1434C------------------------------------------- 1435C Calculate the (vir.vir,vir;del) blocks 1436C contribution to D2ABG directly. 1437C------------------------------------------- 1438C 1439 DO 100 ISYMG = 1,NSYM 1440C 1441 ISYZ2I = MULD2H(MULD2H(ISYMG,ISYMTR),ISYMLP) 1442C 1443 KZ2GAM = KEND1 1444 KEND2 = KZ2GAM + NT2BCD(ISYZ2I) 1445 LWRK2 = LWORK - KEND2 1446C 1447 IF (LWRK2 .LT. 0) THEN 1448 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND2 1449 CALL QUIT 1450 & ('Insufficient space for 3. allocation in CC_DEN2') 1451 ENDIF 1452C 1453 DO 110 G = 1,NBAS(ISYMG) 1454C 1455C--------------------------------------------------------------- 1456C Calculate the gamma backtransformed zeta2-amplitude. 1457C--------------------------------------------------------------- 1458C 1459 ICON = 2 1460 CALL CC_T2AO(Z2AM,XLAMDP,ISYMLP,WORK(KZ2GAM), 1461 * WORK(KEND2),LWRK2, 1462 * G,ISYMG,ISYMTR,ICON) 1463C 1464C-------------------------------------- 1465C Calculate the contribution. 1466C-------------------------------------- 1467C 1468 CALL DABGAO(D2ABG,ISYDEN,WORK(KT2DEL),ISYMTI, 1469 * WORK(KZ2GAM),ISYZ2I,WORK(KEND2),LWRK2, 1470 * G,ISYMG) 1471C 1472 110 CONTINUE 1473 100 CONTINUE 1474C 1475 ENDIF 1476C 1477C------------------------------------------------------------ 1478C Calculate terms of (occ.vir,occ;del) block with t2-del. 1479C Note that the Z-intermediate is overwritten by the 1480C last calculation! 1481C------------------------------------------------------------ 1482C 1483 ISYMZI = MULD2H(ISYMTI,ISYMTR) 1484C 1485 KD2IAJ = KEND1 1486 KZINT = KD2IAJ + NT2BCD(ISYD2H) 1487 KEND2 = KZINT + NT2BCD(ISYMZI) 1488 LWRK2 = LWORK - KEND2 1489C 1490 IF (LWRK2 .LT. 0) THEN 1491 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND2 1492 CALL QUIT( 1493 * 'Insufficient space for fourth allocation in CC_DEN2') 1494 ENDIF 1495C 1496 CALL DZERO(WORK(KD2IAJ),NT2BCD(ISYD2H)) 1497C 1498 IF (CCSD) THEN 1499C 1500 ICON = 1 1501 CALL CC_ZWVI(WORK(KZINT),Z2AM,ISYMTR,WORK(KT2DEL),ISYMTI, 1502 * WORK(KEND2),LWRK2,ICON) 1503C 1504 CALL DIAJ25(WORK(KD2IAJ),ISYD2H,WORK(KT2DEL),ISYMTI,XMIRES, 1505 * WORK(KEND2),LWRK2) 1506 CALL DIAJ27(WORK(KD2IAJ),ISYD2H,T2AM,WORK(KZINT),ISYMZI, 1507 * WORK(KEND2),LWRK2) 1508C 1509 ICON = 2 1510 CALL CCSD_T2TP(Z2AM,WORK(KEND2),LWRK2,ISYMTR) 1511 CALL CC_ZWVI(WORK(KZINT),Z2AM,ISYMTR,WORK(KT2DEL),ISYMTI, 1512 * WORK(KEND2),LWRK2,ICON) 1513 CALL DIAJ27(WORK(KD2IAJ),ISYD2H,T2AM,WORK(KZINT),ISYMZI, 1514 * WORK(KEND2),LWRK2) 1515 CALL CCSD_T2TP(Z2AM,WORK(KEND2),LWRK2,ISYMTR) 1516C 1517 ENDIF 1518C 1519C-------------------------------------------------------------- 1520C Calculate delta backtransformed T2AMT & concomitant ZINT. 1521C Note that these overwrite the old intermediates, since 1522C all terms in need of the old ones have been calculated! 1523C-------------------------------------------------------------- 1524C 1525 CALL CC_TI(WORK(KT2DEL),ISYMTI,T2AMT,ISYMOP,XLAMDH,ISYMLH, 1526 * WORK(KEND2),LWRK2,IDEL,ISYMD) 1527C 1528 IF (CCSD) THEN 1529C 1530 ICON = 1 1531 CALL CC_ZWVI(WORK(KZINT),Z2AM,ISYMTR,WORK(KT2DEL),ISYMTI, 1532 * WORK(KEND2),LWRK2,ICON) 1533C 1534C-------------------------------------------------------------------- 1535C Calculate terms of (occ.vir,occ;del) block with Zbar (in ZINT). 1536C-------------------------------------------------------------------- 1537C 1538 CALL DIAJ26(WORK(KD2IAJ),ISYD2H,T2AM,WORK(KZINT),ISYMZI, 1539 * WORK(KEND2),LWRK2) 1540 CALL DIAJ28(WORK(KD2IAJ),ISYD2H,T2AMT,WORK(KZINT),ISYMZI, 1541 * WORK(KEND2),LWRK2) 1542C 1543 ENDIF 1544C 1545C----------------------------------------------------------------- 1546C Calculate the remaining terms of (occ.vir,occ;del) block 1547C and calculate the remainder of the (vir.occ,occ;del) block. 1548C Note that Zbar is the first contribution to this last block. 1549C----------------------------------------------------------------- 1550C 1551 KD2AIJ = KZINT 1552 KEND2 = KD2AIJ + NT2BCD(ISYD2H) 1553 LWRK2 = LWORK - KEND2 1554C 1555 IF (LWRK2 .LT. 0) THEN 1556 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND2 1557 CALL QUIT( 1558 * 'Insufficient space for fourth allocation in CC_DEN2') 1559 ENDIF 1560C 1561 IF (CC2) CALL DZERO(WORK(KD2AIJ),NT2BCD(ISYD2H)) 1562C 1563C------------------------------------------------ 1564C Contributions from projection against <HF|. 1565C------------------------------------------------ 1566C 1567 IF (ISYMTI.NE.ISYD2H) CALL QUIT('Symmetry mismatch in CC_DEN2') 1568 CALL DAXPY(NT2BCD(ISYMTI),TWO,WORK(KT2DEL),1,WORK(KD2IAJ),1) 1569C 1570C------------------------------------------------ 1571C Contributions from projection against <u1|. 1572C------------------------------------------------ 1573C 1574 CALL DAIJ11(WORK(KD2AIJ),ISYD2H,DAI,ISYD1,XLAMDH,ISYMLH, 1575 * IDEL,ISYMD) 1576 CALL DAIJ12(WORK(KD2AIJ),ISYD2H,DAI,ISYD1,XLAMDH,ISYMLH, 1577 * WORK(KEND2),LWRK2,IDEL,ISYMD) 1578 CALL DIAJ11(WORK(KD2IAJ),ISYD2H,DIA,ISYD1,XLAMDH,ISYMLH, 1579 * IDEL,ISYMD) 1580 CALL DIAJ13(WORK(KD2IAJ),ISYD2H,DAI,ISYD1,T2AMT,XLAMDH,ISYMLH, 1581 * WORK(KEND2),LWRK2,IDEL,ISYMD) 1582C 1583C------------------------------------------------ 1584C Contributions from projection against <u2|. 1585C------------------------------------------------ 1586C 1587 IF (CCSD) THEN 1588C 1589 CALL DAIJ22(WORK(KD2AIJ),ISYD2H,DAB,ISYD1,XLAMDH,ISYMLH, 1590 * WORK(KEND2),LWRK2,IDEL,ISYMD) 1591 CALL DIAJ21(WORK(KD2IAJ),ISYD2H,YMAT,WORK(KT2DEL),ISYMTI) 1592 CALL DIAJ22(WORK(KD2IAJ),ISYD2H,XMAT,WORK(KT2DEL),ISYMTI) 1593 CALL DIAJ23(WORK(KD2IAJ),ISYD2H,XMAT,WORK(KT2DEL),ISYMTI) 1594 CALL DIAJ24(WORK(KD2IAJ),ISYD2H,T2AMT,DAB,ISYD1,XLAMDH,ISYMLH, 1595 * WORK(KEND2),LWRK2,IDEL,ISYMD) 1596C 1597 ENDIF 1598C 1599C------------------------------------------------------------- 1600C Backtransform occupied index to AO and store in results. 1601C------------------------------------------------------------- 1602C 1603 ICON = 2 1604 CALL CCD2_PQAO(D2IAG,ISYDEN,WORK(KD2IAJ),ISYD2H, 1605 * XLAMDP,ISYMLP,ICON) 1606 CALL CCD2_PQAO(D2AIG,ISYDEN,WORK(KD2AIJ),ISYD2H, 1607 * XLAMDP,ISYMLP,ICON) 1608C 1609C-------------------------------------------------- 1610C Calculate the three blocks: (occ.vir,vir;del) 1611C (vir.occ,vir;del) & (vir.vir,occ;del). 1612C-------------------------------------------------- 1613C 1614 KD2IAB = KEND1 1615 KD2ABI = KD2IAB + NCKATR(ISYD2H) 1616 KD2AIB = KD2ABI + NCKASR(ISYD2H) 1617 KEND2 = KD2AIB + NCKATR(ISYD2H) 1618 LWRK2 = LWORK - KEND2 1619C 1620 IF (LWRK2 .LT. 0) THEN 1621 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND2 1622 CALL QUIT('Insufficient space for fifth allocation in CC_DEN2') 1623 ENDIF 1624C 1625 CALL DZERO(WORK(KD2IAB),NCKATR(ISYD2H)) 1626 CALL DZERO(WORK(KD2ABI),NCKASR(ISYD2H)) 1627 CALL DZERO(WORK(KD2AIB),NCKATR(ISYD2H)) 1628C 1629C-------------------------------------------------------------- 1630C Calculate backtransformed zeta2-amplitude zeta(ai,b;del). 1631C Note that this equals the d(ai,b;del) density block. 1632C-------------------------------------------------------------- 1633C 1634 CALL DAIB21(WORK(KD2AIB),ISYD2H,Z2AM,XLAMDH,ISYMLH, 1635 * WORK(KEND2),LWRK2,IDEL,ISYMD) 1636C 1637C------------------------------------------------------- 1638C Backtransform third index of the (vir.occ,vir;del) 1639C block to AO-basis and store in result. 1640C------------------------------------------------------- 1641C 1642 ICON = 5 1643 CALL CCD2_PQAO(D2AIG,ISYDEN,WORK(KD2AIB),ISYD2H, 1644 * XLAMDP,ISYMLP,ICON) 1645C 1646C----------------------------------------------- 1647C Calculate the zeta(ai,b;del)-containing 1648C contributions to the remaining two blocks. 1649C----------------------------------------------- 1650C 1651 KZ2AO = KD2AIB 1652 ISYZ2AO = ISYD2H 1653C 1654 IF (CCSD) THEN 1655C 1656 CALL DIAC22(WORK(KD2IAB),ISYD2H,T2AMT,WORK(KZ2AO),ISYZ2AO, 1657 * WORK(KEND2),LWRK2) 1658 CALL DABI22(WORK(KD2ABI),ISYD2H,T2AM,WORK(KZ2AO),ISYZ2AO, 1659 * WORK(KEND2),LWRK2) 1660 CALL DABI23(WORK(KD2ABI),ISYD2H,T2AM,WORK(KZ2AO),ISYZ2AO, 1661 * WORK(KEND2),LWRK2) 1662C 1663 ENDIF 1664C 1665C-------------------------------------------------------------------- 1666C Calculate remaining contributions from projection against <u1|. 1667C-------------------------------------------------------------------- 1668C 1669 CALL DABI11(WORK(KD2ABI),ISYD2H,DAI,ISYD1,WORK(KT2DEL),ISYMTI) 1670 CALL DIAC11(WORK(KD2IAB),ISYD2H,DAI,ISYD1,WORK(KT2DEL),ISYMTI) 1671C 1672C-------------------------------------------------------------------- 1673C Calculate remaining contributions from projection against <u2|. 1674C-------------------------------------------------------------------- 1675C 1676 IF (CCSD) THEN 1677C 1678 CALL DABI21(WORK(KD2ABI),ISYD2H,DAB,ISYD1,XLAMDH,ISYMLH, 1679 & IDEL,ISYMD) 1680 CALL DIAC21(WORK(KD2IAB),ISYD2H,YMAT,ISYD1,XLAMDH,ISYMLH, 1681 & IDEL,ISYMD) 1682C 1683 ENDIF 1684C 1685C--------------------------------------------------- 1686C Backtransform third index of the two remaining 1687C blocks to AO and store in results. 1688C--------------------------------------------------- 1689C 1690 ICON = 5 1691 CALL CCD2_PQAO(D2IAG,ISYDEN,WORK(KD2IAB),ISYD2H, 1692 * XLAMDP,ISYMLP,ICON) 1693 ICON = 3 1694 CALL CCD2_PQAO(D2ABG,ISYDEN,WORK(KD2ABI),ISYD2H, 1695 * XLAMDP,ISYMLP,ICON) 1696C 1697C--------------------------------------------------- 1698C Calculate the (occ.occ,occ;del) density block. 1699C--------------------------------------------------- 1700C 1701 KD2IJK = KEND1 1702 KEND2 = KD2IJK + NMAIJK(ISYD2H) 1703 LWRK2 = LWORK - KEND2 1704C 1705 IF (LWRK2 .LT. 0) THEN 1706 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND2 1707 CALL QUIT('Insufficient space for seventh allocation '// 1708 & 'in CC_DEN2') 1709 ENDIF 1710C 1711 CALL DZERO(WORK(KD2IJK),NMAIJK(ISYD2H)) 1712C 1713C------------------------------------------------ 1714C Contributions from projection against <HF|. 1715C------------------------------------------------ 1716C 1717 CALL DIJK01(WORK(KD2IJK),ISYD2H,XLAMDH,ISYMLH,IDEL,ISYMD) 1718 CALL DIJK02(WORK(KD2IJK),ISYD2H,XLAMDH,ISYMLH,IDEL,ISYMD) 1719C 1720C------------------------------------------------ 1721C Contributions from projection against <u1|. 1722C------------------------------------------------ 1723C 1724 CALL DIJK11(WORK(KD2IJK),ISYD2H,XLAMDH,ISYMLH,DIA,ISYD1, 1725 & WORK(KEND2),LWRK2,IDEL,ISYMD) 1726 CALL DIJK13(WORK(KD2IJK),ISYD2H,DAI,ISYD1,WORK(KT2DEL),ISYMTI) 1727C 1728C------------------------------------------------ 1729C Contributions from projection against <u2|. 1730C------------------------------------------------ 1731C 1732 IF (CCSD) THEN 1733C 1734 CALL DIJK21(WORK(KD2IJK),ISYD2H,XMAT,XLAMDH,ISYMLH, 1735 * WORK(KEND2),LWRK2,IDEL,ISYMD) 1736 CALL DIJK23(WORK(KD2IJK),ISYD2H,XMAT,XLAMDH,ISYMLH,IDEL,ISYMD) 1737 CALL DIJK24(WORK(KD2IJK),ISYD2H,XMAT,XLAMDH,ISYMLH,IDEL,ISYMD) 1738 CALL DIJK25(WORK(KD2IJK),ISYD2H,XMINT,XLAMDH,ISYMLH,IDEL,ISYMD) 1739C 1740 ENDIF 1741C 1742C------------------------------------------------------------------ 1743C Backtransform third occupied index to AO and store in result. 1744C------------------------------------------------------------------ 1745C 1746 ICON = 1 1747 CALL CCD2_PQAO(D2IJG,ISYDEN,WORK(KD2IJK),ISYD2H, 1748 * XLAMDP,ISYMLP,ICON) 1749C 1750 IF (DEBUG) WRITE(LUPRI,*) 'Leave CC_DEN2' 1751 CALL QEXIT('CC_DEN2') 1752C 1753 RETURN 1754 END 1755C /* Deck ccd2_pqao */ 1756 SUBROUTINE CCD2_PQAO(D2PQG,ISYPQG,D2PQMO,ISYDMO, 1757 * XLAMDP,ISYMLP,ICON) 1758C 1759C Written by Asger Halkier 19/2 - 1996 1760C 1761C Version: 1.0 1762C 1763C Purpose: To calculate the backtransformation of the CC density 1764C d(pq,r;del) to d(pq,gam;del), where p, q, & r are 1765C passed through the variable ICON! 1766C 1767#if defined (IMPLICIT_NONE) 1768 IMPLICIT NONE 1769#else 1770# include "implicit.h" 1771#endif 1772#include "priunit.h" 1773#include "ccorb.h" 1774#include "ccsdsym.h" 1775 1776#if defined (SYS_CRAY) 1777 REAL ZERO, ONE 1778#else 1779 DOUBLE PRECISION ZERO, ONE 1780#endif 1781 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 1782 1783 INTEGER ISYPQG, ISYDMO, ISYMLP, ICON 1784 1785#if defined (SYS_CRAY) 1786 REAL D2PQG(*), D2PQMO(*), XLAMDP(*) 1787#else 1788 DOUBLE PRECISION D2PQG(*), D2PQMO(*), XLAMDP(*) 1789#endif 1790 1791 INTEGER ISYMK, ISYMG, ISYMA, ISYMB, KOFF1, KOFF2, KOFF3, 1792 * NTOTIJ, NTOTAI, NTOTAB, NTOTG, ISYMIJ, ISYMAB, ISYMAI 1793C 1794 CALL QENTER('CCD2_PQAO') 1795C 1796 IF (MULD2H(ISYDMO,ISYMLP).NE.ISYPQG) THEN 1797 CALL QUIT('Symmetry mismatch in CCD2_PQAO') 1798 END IF 1799C 1800C----------------------------------------------------------------- 1801C Calculate the transformation of the (occ.occ,occ;del) block. 1802C----------------------------------------------------------------- 1803C 1804 IF (ICON .EQ. 1) THEN 1805C 1806 DO 100 ISYMIJ = 1,NSYM 1807C 1808 ISYMK = MULD2H(ISYMIJ,ISYDMO) 1809 ISYMG = MULD2H(ISYMK,ISYMLP) 1810C 1811 KOFF1 = IMAIJK(ISYMIJ,ISYMK) + 1 1812 KOFF2 = IGLMRH(ISYMG,ISYMK) + 1 1813 KOFF3 = ID2IJG(ISYMIJ,ISYMG) + 1 1814C 1815 NTOTIJ = MAX(NMATIJ(ISYMIJ),1) 1816 NTOTG = MAX(NBAS(ISYMG),1) 1817C 1818 CALL DGEMM('N','T',NMATIJ(ISYMIJ),NBAS(ISYMG),NRHF(ISYMK), 1819 * ONE,D2PQMO(KOFF1),NTOTIJ,XLAMDP(KOFF2),NTOTG, 1820 * ONE,D2PQG(KOFF3),NTOTIJ) 1821C 1822 100 CONTINUE 1823C 1824C--------------------------------------------------------------- 1825C Calculate the transformation of both the (occ.vir,occ;del) 1826C block and the (vir.occ,occ;del) block (stored equally). 1827C--------------------------------------------------------------- 1828C 1829 ELSE IF (ICON .EQ. 2) THEN 1830C 1831 DO 110 ISYMAI = 1,NSYM 1832C 1833 ISYMK = MULD2H(ISYMAI,ISYDMO) 1834 ISYMG = MULD2H(ISYMK,ISYMLP) 1835C 1836 KOFF1 = IT2BCD(ISYMAI,ISYMK) + 1 1837 KOFF2 = IGLMRH(ISYMG,ISYMK) + 1 1838 KOFF3 = ID2AIG(ISYMAI,ISYMG) + 1 1839C 1840 NTOTAI = MAX(NT1AM(ISYMAI),1) 1841 NTOTG = MAX(NBAS(ISYMG),1) 1842C 1843 CALL DGEMM('N','T',NT1AM(ISYMAI),NBAS(ISYMG),NRHF(ISYMK), 1844 * ONE,D2PQMO(KOFF1),NTOTAI,XLAMDP(KOFF2),NTOTG, 1845 * ONE,D2PQG(KOFF3),NTOTAI) 1846C 1847 110 CONTINUE 1848C 1849C----------------------------------------------------------------- 1850C Calculate the transformation of the (vir.vir,occ;del) block. 1851C----------------------------------------------------------------- 1852C 1853 ELSE IF (ICON .EQ. 3) THEN 1854C 1855 DO 120 ISYMAB = 1,NSYM 1856C 1857 ISYMK = MULD2H(ISYMAB,ISYDMO) 1858 ISYMG = MULD2H(ISYMK,ISYMLP) 1859C 1860 KOFF1 = ICKASR(ISYMAB,ISYMK) + 1 1861 KOFF2 = IGLMRH(ISYMG,ISYMK) + 1 1862 KOFF3 = ID2ABG(ISYMAB,ISYMG) + 1 1863C 1864 NTOTAB = MAX(NMATAB(ISYMAB),1) 1865 NTOTG = MAX(NBAS(ISYMG),1) 1866C 1867 CALL DGEMM('N','T',NMATAB(ISYMAB),NBAS(ISYMG),NRHF(ISYMK), 1868 * ONE,D2PQMO(KOFF1),NTOTAB,XLAMDP(KOFF2),NTOTG, 1869 * ONE,D2PQG(KOFF3),NTOTAB) 1870C 1871 120 CONTINUE 1872C 1873C----------------------------------------------------------------- 1874C Calculate the transformation of the (occ.occ,vir;del) block. 1875C----------------------------------------------------------------- 1876C 1877 ELSE IF (ICON .EQ. 4) THEN 1878C 1879 DO 130 ISYMIJ = 1,NSYM 1880C 1881 ISYMA = MULD2H(ISYMIJ,ISYDMO) 1882 ISYMG = MULD2H(ISYMA,ISYMLP) 1883C 1884 KOFF1 = IMAIJA(ISYMIJ,ISYMA) + 1 1885 KOFF2 = IGLMVI(ISYMG,ISYMA) + 1 1886 KOFF3 = ID2IJG(ISYMIJ,ISYMG) + 1 1887C 1888 NTOTIJ = MAX(NMATIJ(ISYMIJ),1) 1889 NTOTG = MAX(NBAS(ISYMG),1) 1890C 1891 CALL DGEMM('N','T',NMATIJ(ISYMIJ),NBAS(ISYMG),NVIR(ISYMA), 1892 * ONE,D2PQMO(KOFF1),NTOTIJ,XLAMDP(KOFF2),NTOTG, 1893 * ONE,D2PQG(KOFF3),NTOTIJ) 1894C 1895 130 CONTINUE 1896C 1897C----------------------------------------------------------------- 1898C Calculate the transformation of the (occ.vir,vir;del) block. 1899C----------------------------------------------------------------- 1900C 1901 ELSE IF (ICON .EQ. 5) THEN 1902C 1903 DO 140 ISYMAI = 1,NSYM 1904C 1905 ISYMB = MULD2H(ISYMAI,ISYDMO) 1906 ISYMG = MULD2H(ISYMB,ISYMLP) 1907C 1908 KOFF1 = ICKATR(ISYMAI,ISYMB) + 1 1909 KOFF2 = IGLMVI(ISYMG,ISYMB) + 1 1910 KOFF3 = ID2AIG(ISYMAI,ISYMG) + 1 1911C 1912 NTOTAI = MAX(NT1AM(ISYMAI),1) 1913 NTOTG = MAX(NBAS(ISYMG),1) 1914C 1915 CALL DGEMM('N','T',NT1AM(ISYMAI),NBAS(ISYMG),NVIR(ISYMB), 1916 * ONE,D2PQMO(KOFF1),NTOTAI,XLAMDP(KOFF2),NTOTG, 1917 * ONE,D2PQG(KOFF3),NTOTAI) 1918C 1919 140 CONTINUE 1920C 1921 ENDIF 1922C 1923 CALL QEXIT('CCD2_PQAO') 1924C 1925 RETURN 1926 END 1927C /* Deck diac11 */ 1928 SUBROUTINE DIAC11(D2IAB,ISYIAB,DAI,ISYD1,T2TIN,ISYMTI) 1929C 1930C Written by Asger Halkier 20/2 - 1996 1931C 1932C Version: 1.0 1933C 1934C Purpose: To calculate the one and only contribution to D2IAB 1935C from projection against singles! 1936C 1937#include "implicit.h" 1938 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 1939 DIMENSION D2IAB(*), DAI(*), T2TIN(*) 1940#include "priunit.h" 1941#include "ccorb.h" 1942#include "ccsdsym.h" 1943C 1944 CALL QENTER('DIAC11') 1945C 1946 IF (MULD2H(ISYMTI,ISYD1).NE.ISYIAB) 1947 & CALL QUIT('Symmetry mismatch in DIAC11') 1948C 1949C---------------------------- 1950C Calculate contribution. 1951C---------------------------- 1952C 1953 DO 100 ISYMAI = 1,NSYM 1954C 1955 ISYMM = MULD2H(ISYMAI,ISYMTI) 1956 ISYMC = MULD2H(ISYMM,ISYD1) 1957C 1958 KOFF1 = IT2BCD(ISYMAI,ISYMM) + 1 1959 KOFF2 = IT1AM(ISYMC,ISYMM) + 1 1960 KOFF3 = ICKATR(ISYMAI,ISYMC) + 1 1961C 1962 NTOTAI = MAX(NT1AM(ISYMAI),1) 1963 NTOTC = MAX(NVIR(ISYMC),1) 1964C 1965 CALL DGEMM('N','T',NT1AM(ISYMAI),NVIR(ISYMC),NRHF(ISYMM),ONE, 1966 * T2TIN(KOFF1),NTOTAI,DAI(KOFF2),NTOTC,ONE, 1967 * D2IAB(KOFF3),NTOTAI) 1968C 1969 100 CONTINUE 1970C 1971 CALL QEXIT('DIAC11') 1972C 1973 RETURN 1974 END 1975C /* Deck dija11 */ 1976 SUBROUTINE DIJA11(D2IJA,ISYDEN,DAI,ISYDAI,XLAMDH,ISYMLH, 1977 & WORK,LWORK,IDEL,ISYMD) 1978C 1979C Written by Asger Halkier 20/2 - 1996 1980C 1981C Version: 1.0 1982C 1983C Purpose: To calculate the first and second contribution to D2IJA 1984C from projection against singles! 1985C 1986#if defined (IMPLICIT_NONE) 1987 IMPLICIT NONE 1988#else 1989# include "implicit.h" 1990#endif 1991#include "priunit.h" 1992#include "ccorb.h" 1993#include "ccsdsym.h" 1994 1995#if defined (SYS_CRAY) 1996 REAL ZERO, ONE, TWO 1997#else 1998 DOUBLE PRECISION ZERO, ONE, TWO 1999#endif 2000 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 2001 2002 INTEGER ISYDEN, ISYDAI, ISYMLH, LWORK, IDEL, ISYMD 2003 INTEGER ISYMK, ISYMA, ISYMIJ, ISYMI, ISYMJ, ISYMAJ, 2004 & KOFF1, KOFF2, NTOTA, NIJA, NDEI, NAJ 2005 2006#if defined (SYS_CRAY) 2007 REAL D2IJA(*), DAI(*), XLAMDH(*), WORK(LWORK) 2008#else 2009 DOUBLE PRECISION D2IJA(*), DAI(*), XLAMDH(*), WORK(LWORK) 2010#endif 2011 2012 CALL QENTER('DIJA11') 2013C 2014C 2015 IF (ISYDAI.NE.1) CALL QUIT('Illegal ISYDAI in DIJA11') 2016C 2017 ISYMK = MULD2H(ISYMLH,ISYMD) 2018 ISYMA = MULD2H(ISYDAI,ISYMK) 2019C 2020 IF (ISYMA.NE.ISYDEN) THEN 2021 CALL QUIT('Symmetry mismatch in DIJA11. (2)') 2022 END IF 2023C 2024C-------------------------------------------------------- 2025C Calculate contraction of D(ai) & lamda hole matrix. 2026C (back transform occupied index of D(ai) with lambda: 2027c d(a d) = sum_k D(ak) XL(dk) 2028C-------------------------------------------------------- 2029C 2030 IF (LWORK .LT. NVIR(ISYMA)) THEN 2031 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', NVIR(ISYMA) 2032 CALL QUIT('Insufficient space for allocation in DIJA11') 2033 ENDIF 2034C 2035 CALL DZERO(WORK,NVIR(ISYMA)) 2036C 2037 KOFF1 = IT1AM(ISYMA,ISYMK) + 1 2038 KOFF2 = IGLMRH(ISYMD,ISYMK) + IDEL - IBAS(ISYMD) 2039C 2040 NTOTA = MAX(NVIR(ISYMA),1) 2041C 2042 CALL DGEMV('N',NVIR(ISYMA),NRHF(ISYMK),ONE,DAI(KOFF1),NTOTA, 2043 * XLAMDH(KOFF2),NBAS(ISYMD),ZERO,WORK,1) 2044C 2045C-------------------------------------- 2046C Calculate the first contribution. 2047C-------------------------------------- 2048C 2049 ISYMIJ = 1 2050C 2051 DO 100 A = 1,NVIR(ISYMA) 2052C 2053 DO 110 ISYMJ = 1,NSYM 2054C 2055 ISYMI = ISYMJ 2056C 2057 DO 120 J = 1,NRHF(ISYMJ) 2058C 2059 I = J 2060 NIJA = IMAIJA(ISYMIJ,ISYMA) + NMATIJ(ISYMIJ)*(A - 1) 2061 * + IMATIJ(ISYMI,ISYMJ) + NRHF(ISYMI)*(J - 1) + I 2062C 2063 D2IJA(NIJA) = D2IJA(NIJA) + TWO*WORK(A) 2064C 2065 120 CONTINUE 2066 110 CONTINUE 2067 100 CONTINUE 2068C 2069C--------------------------------------- 2070C Calculate the second contribution. 2071C--------------------------------------- 2072C 2073 ISYMAJ = ISYDAI 2074 ISYMI = MULD2H(ISYMLH,ISYMD) 2075C 2076 IF (ISYMI.NE.ISYDEN) THEN 2077 CALL QUIT('Symmetry mismatch in DIJA11. (1)') 2078 END IF 2079C 2080 DO 130 I = 1,NRHF(ISYMI) 2081C 2082 NDEI = IGLMRH(ISYMD,ISYMI) + NBAS(ISYMD)*(I - 1) 2083 * + IDEL - IBAS(ISYMD) 2084C 2085 DO 140 ISYMA = 1,NSYM 2086C 2087 ISYMJ = MULD2H(ISYMA,ISYMAJ) 2088 ISYMIJ = MULD2H(ISYMI,ISYMJ) 2089C 2090 DO 150 J = 1,NRHF(ISYMJ) 2091C 2092 DO 160 A = 1,NVIR(ISYMA) 2093C 2094 NAJ = IT1AM(ISYMA,ISYMJ) + NVIR(ISYMA)*(J - 1) + A 2095 NIJA = IMAIJA(ISYMIJ,ISYMA) + NMATIJ(ISYMIJ)*(A - 1) 2096 * + IMATIJ(ISYMI,ISYMJ) + NRHF(ISYMI)*(J - 1) + I 2097C 2098 D2IJA(NIJA) = D2IJA(NIJA) - DAI(NAJ)*XLAMDH(NDEI) 2099C 2100 160 CONTINUE 2101 150 CONTINUE 2102 140 CONTINUE 2103 130 CONTINUE 2104C 2105 CALL QEXIT('DIJA11') 2106C 2107 RETURN 2108 END 2109C /* Deck diac21 */ 2110 SUBROUTINE DIAC21(D2IAC,ISYIAC,YMAT,ISYMY,XLAMDH,ISYMLH, 2111 * IDEL,ISYMD) 2112C 2113C Written by Asger Halkier 21/2 - 1996 2114C 2115C Version: 1.0 2116C 2117C Purpose: To calculate the first contribution to D2IAC 2118C from projection against doubles! 2119C 2120#include "implicit.h" 2121 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 2122 DIMENSION D2IAC(*), YMAT(*), XLAMDH(*) 2123#include "priunit.h" 2124#include "ccorb.h" 2125#include "ccsdsym.h" 2126C 2127 CALL QENTER('DIAC21') 2128C 2129 ISYMI = MULD2H(ISYMD,ISYMLH) 2130 ISYMAC = ISYMY 2131C 2132 IF (MULD2H(ISYMI,ISYMAC).NE.ISYIAC) 2133 * CALL QUIT('Symmetry mismatch in DIAC21') 2134C 2135C-------------------------------- 2136C Calculate the contribution. 2137C-------------------------------- 2138C 2139 DO 100 I = 1,NRHF(ISYMI) 2140C 2141 NDEI = IGLMRH(ISYMD,ISYMI) + NBAS(ISYMD)*(I - 1) 2142 * + IDEL - IBAS(ISYMD) 2143C 2144 FACT = -XLAMDH(NDEI) 2145C 2146 DO 110 ISYMA = 1,NSYM 2147C 2148 ISYMC = MULD2H(ISYMA,ISYMAC) 2149 ISYMAI = MULD2H(ISYMA,ISYMI) 2150C 2151 DO 120 C = 1,NVIR(ISYMC) 2152C 2153 NAC = IMATAB(ISYMA,ISYMC) + NVIR(ISYMA)*(C - 1) + 1 2154 NAIC = ICKATR(ISYMAI,ISYMC) + NT1AM(ISYMAI)*(C - 1) 2155 * + IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I - 1) + 1 2156C 2157 CALL DAXPY(NVIR(ISYMA),FACT,YMAT(NAC),1, 2158 * D2IAC(NAIC),1) 2159C 2160 120 CONTINUE 2161 110 CONTINUE 2162 100 CONTINUE 2163C 2164 CALL QEXIT('DIAC21') 2165C 2166 RETURN 2167 END 2168C /* Deck dija21 */ 2169 SUBROUTINE DIJA21(D2IJA,ISYDEN,DAB,ISYDAB,XLAMDH,ISYMLH, 2170 & WORK,LWORK,IDEL,ISYMD) 2171C 2172C Written by Asger Halkier 21/2 - 1996 2173C 2174C Version: 1.0 2175C 2176C Purpose: To calculate the first contribution to D2IJA 2177C from projection against doubles! 2178C 2179#if defined (IMPLICIT_NONE) 2180 IMPLICIT NONE 2181#else 2182# include "implicit.h" 2183#endif 2184#include "priunit.h" 2185#include "ccorb.h" 2186#include "ccsdsym.h" 2187C 2188#if defined (SYS_CRAY) 2189 REAL ZERO, ONE, TWO 2190#else 2191 DOUBLE PRECISION ZERO, ONE, TWO 2192#endif 2193 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 2194 2195 INTEGER ISYDEN, ISYDAB, ISYMLH, IDEL, ISYMD, LWORK 2196 INTEGER ISYMB, ISYMA, ISYMIJ, ISYMI, ISYMJ, 2197 & KOFF1, KOFF2, NTOTA, NIJA 2198 2199#if defined (SYS_CRAY) 2200 REAL D2IJA(*), DAB(*), XLAMDH(*), WORK(LWORK) 2201#else 2202 DOUBLE PRECISION D2IJA(*), DAB(*), XLAMDH(*), WORK(LWORK) 2203#endif 2204C 2205 CALL QENTER('DIJA21') 2206C 2207 IF (ISYDAB.NE.1) CALL QUIT('Illegal ISYDAI in DIJA21') 2208C 2209 ISYMB = MULD2H(ISYMLH,ISYMD) 2210 ISYMA = MULD2H(ISYDAB,ISYMB) 2211 ISYMIJ = 1 2212C 2213 IF (ISYMA.NE.ISYDEN) THEN 2214 CALL QUIT('Symmetry mismatch in DIJA21. (2)') 2215 END IF 2216C 2217C------------------------------------------------------------------- 2218C Calculate contraction of transposed Y-matrix and lamda matrix. 2219C------------------------------------------------------------------- 2220C 2221 IF (LWORK .LT. NVIR(ISYMA)) THEN 2222 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', NVIR(ISYMA) 2223 CALL QUIT('Insufficient space for allocation in DIJA21') 2224 ENDIF 2225C 2226 CALL DZERO(WORK,NVIR(ISYMA)) 2227C 2228 KOFF1 = IMATAB(ISYMA,ISYMB) + 1 2229 KOFF2 = IGLMVI(ISYMD,ISYMB) + IDEL - IBAS(ISYMD) 2230C 2231 NTOTA = MAX(NVIR(ISYMA),1) 2232C 2233 CALL DGEMV('N',NVIR(ISYMA),NVIR(ISYMB),ONE,DAB(KOFF1),NTOTA, 2234 * XLAMDH(KOFF2),NBAS(ISYMD),ZERO,WORK,1) 2235C 2236 DO 100 A = 1,NVIR(ISYMA) 2237C 2238 DO 110 ISYMI = 1,NSYM 2239C 2240 ISYMJ = ISYMI 2241C 2242 DO 120 I = 1,NRHF(ISYMI) 2243C 2244 J = I 2245 NIJA = IMAIJA(ISYMIJ,ISYMA) + NMATIJ(ISYMIJ)*(A - 1) 2246 * + IMATIJ(ISYMI,ISYMJ) + NRHF(ISYMI)*(J - 1) + I 2247C 2248 D2IJA(NIJA) = D2IJA(NIJA) + TWO*WORK(A) 2249C 2250 120 CONTINUE 2251 110 CONTINUE 2252 100 CONTINUE 2253C 2254 CALL QEXIT('DIJA21') 2255C 2256 RETURN 2257 END 2258C /* Deck dija22 */ 2259 SUBROUTINE DIJA22(D2IJA,ISYDEN,Z2AM,T2INT,ISYMTI,WORK,LWORK) 2260C 2261C Written by Asger Halkier 21/2 - 1996 2262C 2263C Version: 1.0 2264C 2265C Purpose: To calculate the second contribution to D2IJA 2266C from projection against doubles! 2267C 2268#include "implicit.h" 2269 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 2270 DIMENSION D2IJA(*), Z2AM(*), T2INT(*), WORK(LWORK) 2271#include "priunit.h" 2272#include "ccorb.h" 2273#include "ccsdsym.h" 2274C 2275 CALL QENTER('DIJA22') 2276C 2277 IF (ISYDEN.NE.ISYMTI) CALL QUIT('Symmetry mismatch in DIJA22') 2278C 2279 ISYCKI = ISYDEN 2280 ISYIJA = ISYDEN 2281C 2282 DO 100 ISYMI = 1,NSYM 2283C 2284 ISYMCK = MULD2H(ISYMI,ISYCKI) 2285 ISYMAJ = ISYMCK 2286C 2287C------------------------------ 2288C Work space allocation. 2289C------------------------------ 2290C 2291 KSCR = 1 2292 KEND = KSCR + NT1AM(ISYMAJ)*NRHF(ISYMI) 2293 LWRK = LWORK - KEND 2294C 2295 IF (LWRK .LT. 0) THEN 2296 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND 2297 CALL QUIT('Insufficient work space for '// 2298 & 'allocation in DIJA22') 2299 ENDIF 2300C 2301 CALL DZERO(WORK(KSCR),NT1AM(ISYMAJ)*NRHF(ISYMI)) 2302C 2303C-------------------------------------------- 2304C Calculate contraction of zeta and t. 2305C-------------------------------------------- 2306C 2307 KOFF1 = IT2SQ(ISYMAJ,ISYMCK) + 1 2308 KOFF2 = IT2BCD(ISYMCK,ISYMI) + 1 2309C 2310 NTOTAJ = MAX(NT1AM(ISYMAJ),1) 2311 NTOTCK = MAX(NT1AM(ISYMCK),1) 2312C 2313 CALL DGEMM('N','N',NT1AM(ISYMAJ),NRHF(ISYMI),NT1AM(ISYMCK), 2314 * ONE,Z2AM(KOFF1),NTOTAJ,T2INT(KOFF2),NTOTCK,ZERO, 2315 * WORK(KSCR),NTOTAJ) 2316C 2317C--------------------------------------------------------- 2318C Store properly reordered scratch-array in result. 2319C--------------------------------------------------------- 2320C 2321 DO 110 ISYMA = 1,NSYM 2322C 2323 ISYMJ = MULD2H(ISYMA,ISYMAJ) 2324 ISYMIJ = MULD2H(ISYMJ,ISYMI) 2325C 2326 DO 120 I = 1,NRHF(ISYMI) 2327C 2328 DO 130 J = 1,NRHF(ISYMJ) 2329C 2330 DO 140 A = 1,NVIR(ISYMA) 2331C 2332 NAJI = NT1AM(ISYMAJ)*(I - 1) + IT1AM(ISYMA,ISYMJ) 2333 * + NVIR(ISYMA)*(J - 1) + A 2334 NIJA = IMAIJA(ISYMIJ,ISYMA) 2335 * + NMATIJ(ISYMIJ)*(A - 1) 2336 * + IMATIJ(ISYMI,ISYMJ) 2337 * + NRHF(ISYMI)*(J - 1) + I 2338C 2339 D2IJA(NIJA) = D2IJA(NIJA) - WORK(KSCR + NAJI - 1) 2340C 2341 140 CONTINUE 2342 130 CONTINUE 2343 120 CONTINUE 2344 110 CONTINUE 2345 100 CONTINUE 2346C 2347 CALL QEXIT('DIJA22') 2348C 2349 RETURN 2350 END 2351C /* Deck dija23 */ 2352 SUBROUTINE DIJA23(D2IJA,ISYDEN,Z2AM,T2INT,ISYMTI,WORK,LWORK) 2353C 2354C Written by Asger Halkier 21/2 - 1996 2355C 2356C Version: 1.0 2357C 2358C Purpose: To calculate the third contribution to D2IJA 2359C from projection against doubles! 2360C 2361#include "implicit.h" 2362 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 2363 DIMENSION D2IJA(*), Z2AM(*), T2INT(*), WORK(LWORK) 2364#include "priunit.h" 2365#include "ccorb.h" 2366#include "ccsdsym.h" 2367C 2368 CALL QENTER('DIJA23') 2369C 2370 IF (ISYDEN.NE.ISYMTI) CALL QUIT('Symmetry mismatch in DIJA23') 2371C 2372 ISYCIK = ISYDEN 2373 ISYIJA = ISYDEN 2374C 2375 DO 100 ISYMI = 1,NSYM 2376C 2377 ISYMCK = MULD2H(ISYMI,ISYCIK) 2378 ISYMAJ = ISYMCK 2379C 2380C---------------------------------- 2381C Work space allocation one. 2382C---------------------------------- 2383C 2384 KSCR1 = 1 2385 KEND1 = KSCR1 + NRHF(ISYMI)*NT1AM(ISYMCK) 2386 LWRK1 = LWORK - KEND1 2387C 2388 IF (LWRK1 .LT. 0) THEN 2389 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 2390 CALL QUIT( 2391 * 'Insufficient work space for allocation in DIJA23') 2392 ENDIF 2393C 2394 CALL DZERO(WORK(KSCR1),NRHF(ISYMI)*NT1AM(ISYMCK)) 2395C 2396C------------------------------------------------------------- 2397C Reorder backtransformed t-amplitude in scratch-array. 2398C------------------------------------------------------------- 2399C 2400 DO 110 ISYMC = 1,NSYM 2401C 2402 ISYMK = MULD2H(ISYMC,ISYMCK) 2403 ISYMCI = MULD2H(ISYMC,ISYMI) 2404C 2405 DO 120 K = 1,NRHF(ISYMK) 2406C 2407 DO 130 I = 1,NRHF(ISYMI) 2408C 2409 DO 140 C = 1,NVIR(ISYMC) 2410C 2411 NCIK = IT2BCD(ISYMCI,ISYMK) 2412 * + NT1AM(ISYMCI)*(K - 1) + IT1AM(ISYMC,ISYMI) 2413 * + NVIR(ISYMC)*(I - 1) + C 2414 NCK = IT1AM(ISYMC,ISYMK) 2415 * + NVIR(ISYMC)*(K - 1) + C 2416 NICK = KSCR1 + NRHF(ISYMI)*(NCK - 1) + I - 1 2417C 2418 WORK(NICK) = T2INT(NCIK) 2419C 2420 140 CONTINUE 2421 130 CONTINUE 2422 120 CONTINUE 2423 110 CONTINUE 2424C 2425 DO 150 ISYMA = 1,NSYM 2426C 2427 ISYMJ = MULD2H(ISYMA,ISYMAJ) 2428 ISYMIJ = MULD2H(ISYMA,ISYDEN) 2429C 2430C------------------------------------- 2431C Work space allocation two. 2432C------------------------------------- 2433C 2434 KSCR2 = KEND1 2435 KEND2 = KSCR2 + NT1AM(ISYMCK)*NRHF(ISYMJ) 2436 LWRK2 = LWORK - KEND2 2437C 2438 IF (LWRK2 .LT. 0) THEN 2439 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND2 2440 CALL QUIT('Insufficient work space for '// 2441 & 'allocation in DIJA23') 2442 ENDIF 2443C 2444 DO 160 A = 1,NVIR(ISYMA) 2445C 2446 CALL DZERO(WORK(KSCR2),NT1AM(ISYMCK)*NRHF(ISYMJ)) 2447C 2448C--------------------------------------------------- 2449C Copy submatrix out of zeta-amplitude. 2450C--------------------------------------------------- 2451C 2452 DO 170 ISYMC = 1,NSYM 2453C 2454 ISYMCJ = MULD2H(ISYMC,ISYMJ) 2455 ISYMAK = ISYMCJ 2456 ISYMK = MULD2H(ISYMC,ISYMCK) 2457C 2458 DO 180 K = 1,NRHF(ISYMK) 2459C 2460 NAK = IT1AM(ISYMA,ISYMK) + NVIR(ISYMA)*(K - 1) + A 2461C 2462 DO 190 J = 1,NRHF(ISYMJ) 2463C 2464 NCJAK = IT2SQ(ISYMCJ,ISYMAK) 2465 * + NT1AM(ISYMCJ)*(NAK - 1) 2466 * + IT1AM(ISYMC,ISYMJ) 2467 * + NVIR(ISYMC)*(J - 1) + 1 2468 NCKJ = KSCR2 + NT1AM(ISYMCK)*(J - 1) 2469 * + IT1AM(ISYMC,ISYMK) 2470 * + NVIR(ISYMC)*(K - 1) 2471C 2472 CALL DCOPY(NVIR(ISYMC),Z2AM(NCJAK),1, 2473 * WORK(NCKJ),1) 2474C 2475 190 CONTINUE 2476 180 CONTINUE 2477 170 CONTINUE 2478C 2479C------------------------------------------------------ 2480C Final contraction and storage in result. 2481C------------------------------------------------------ 2482C 2483 KOFF1 = IMAIJA(ISYMIJ,ISYMA) + NMATIJ(ISYMIJ)*(A - 1) 2484 * + IMATIJ(ISYMI,ISYMJ) + 1 2485C 2486 NTOTI = MAX(NRHF(ISYMI),1) 2487 NTOTCK = MAX(NT1AM(ISYMCK),1) 2488C 2489 CALL DGEMM('N','N',NRHF(ISYMI),NRHF(ISYMJ), 2490 * NT1AM(ISYMCK),-ONE,WORK(KSCR1),NTOTI, 2491 * WORK(KSCR2),NTOTCK,ONE,D2IJA(KOFF1),NTOTI) 2492C 2493 160 CONTINUE 2494 150 CONTINUE 2495 100 CONTINUE 2496C 2497 CALL QEXIT('DIJA23') 2498C 2499 RETURN 2500 END 2501C /* Deck daib21 */ 2502 SUBROUTINE DAIB21(D2AIB,ISYAIB,Z2AM,XLAMDH,ISYMLH,WORK,LWORK, 2503 * IDEL,ISYMD) 2504C 2505C Written by Asger Halkier 22/2 - 1996 2506C 2507C Version: 1.0 2508C 2509C Purpose: To calculate the one and only contribution to D2AIB 2510C at all which arises from projection against doubles! 2511C 2512#include "implicit.h" 2513 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 2514 DIMENSION D2AIB(*), Z2AM(*), XLAMDH(*), WORK(LWORK) 2515#include "priunit.h" 2516#include "ccorb.h" 2517#include "ccsdsym.h" 2518C 2519 CALL QENTER('DAIB21') 2520C 2521 ISYMJ = MULD2H(ISYMD,ISYMLH) 2522C 2523 IF (MULD2H(ISYMLH,ISYMD).NE.ISYAIB) 2524 & CALL QUIT('Symmetry mismatch in DAIB21') 2525C 2526C------------------------------- 2527C Work space allocation one. 2528C------------------------------- 2529C 2530 KJVEC = 1 2531 KEND1 = KJVEC + NRHF(ISYMJ) 2532 LWRK1 = LWORK - KEND1 2533C 2534 IF (LWRK1 .LT. 0) THEN 2535 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 2536 CALL QUIT('Insufficient work space for allocation in DAIB21') 2537 ENDIF 2538C 2539 CALL DZERO(WORK(KJVEC),NRHF(ISYMJ)) 2540C 2541C------------------------------------------ 2542C Copy vector out of lamda hole matrix. 2543C------------------------------------------ 2544C 2545 KOFF1 = IGLMRH(ISYMD,ISYMJ) + IDEL - IBAS(ISYMD) 2546C 2547 CALL DCOPY(NRHF(ISYMJ),XLAMDH(KOFF1),NBAS(ISYMD),WORK(KJVEC),1) 2548C 2549 DO 100 ISYMB = 1,NSYM 2550C 2551 ISYMAI = MULD2H(ISYMB,ISYAIB) 2552 ISYMBJ = MULD2H(ISYMB,ISYMJ) 2553C 2554C---------------------------------- 2555C Work space allocation two. 2556C---------------------------------- 2557C 2558 KZSUB = KEND1 2559 KEND2 = KZSUB + NT1AM(ISYMAI)*NRHF(ISYMJ) 2560 LWRK2 = LWORK - KEND2 2561C 2562 IF (LWRK2 .LT. 0) THEN 2563 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND2 2564 CALL QUIT( 2565 * 'Insufficient work space for allocation in DAIB21') 2566 ENDIF 2567C 2568 CALL DZERO(WORK(KZSUB),NT1AM(ISYMAI)*NRHF(ISYMJ)) 2569C 2570 DO 110 B = 1,NVIR(ISYMB) 2571C 2572C------------------------------------------------ 2573C Copy submatrix out of zeta-amplitude. 2574C------------------------------------------------ 2575C 2576 DO 120 J = 1,NRHF(ISYMJ) 2577C 2578 NBJ = IT1AM(ISYMB,ISYMJ) + NVIR(ISYMB)*(J - 1) + B 2579 NAIBJ = IT2SQ(ISYMAI,ISYMBJ) 2580 * + NT1AM(ISYMAI)*(NBJ - 1) + 1 2581 NAIJ = KZSUB + NT1AM(ISYMAI)*(J - 1) 2582C 2583 CALL DCOPY(NT1AM(ISYMAI),Z2AM(NAIBJ),1,WORK(NAIJ),1) 2584C 2585 120 CONTINUE 2586C 2587C--------------------------------------------------- 2588C Final contraction and storage in result. 2589C--------------------------------------------------- 2590C 2591 KOFF2 = ICKATR(ISYMAI,ISYMB) + NT1AM(ISYMAI)*(B - 1) + 1 2592C 2593 NTOTAI = MAX(NT1AM(ISYMAI),1) 2594C 2595 CALL DGEMV('N',NT1AM(ISYMAI),NRHF(ISYMJ),ONE,WORK(KZSUB), 2596 * NTOTAI,WORK(KJVEC),1,ONE,D2AIB(KOFF2),1) 2597C 2598 110 CONTINUE 2599 100 CONTINUE 2600C 2601 CALL QEXIT('DAIB21') 2602C 2603 RETURN 2604 END 2605C /* Deck diac22 */ 2606 SUBROUTINE DIAC22(D2IAB,ISYIAC,T2AMT,Z2AO,ISYZ2AO,WORK,LWORK) 2607C 2608C Written by Asger Halkier 24/2 - 1996 2609C 2610C Version: 1.0 2611C 2612C Purpose: To calculate the second contribution to D2IAC 2613C from projection against doubles! 2614C 2615#include "implicit.h" 2616 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 2617 DIMENSION D2IAB(*), T2AMT(*), Z2AO(*), WORK(LWORK) 2618#include "priunit.h" 2619#include "ccorb.h" 2620#include "ccsdsym.h" 2621C 2622 INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 2623C 2624 CALL QENTER('DIAC22') 2625C 2626 DO 100 ISYMI = 1,NSYM 2627C 2628 ISYMAC = MULD2H(ISYMI,ISYIAC) 2629C 2630 DO 110 I = 1,NRHF(ISYMI) 2631C 2632 DO 120 ISYMA = 1,NSYM 2633C 2634 ISYMC = MULD2H(ISYMA,ISYMAC) 2635 ISYMDK = MULD2H(ISYMC,ISYZ2AO) 2636 ISYMAI = MULD2H(ISYMA,ISYMI) 2637C 2638C------------------------------------ 2639C Work space allocation. 2640C------------------------------------ 2641C 2642 KT2SUB = 1 2643 KACRES = KT2SUB + NVIR(ISYMA)*NT1AM(ISYMDK) 2644 KEND1 = KACRES + NVIR(ISYMA)*NVIR(ISYMC) 2645 LWRK1 = LWORK - KEND1 2646C 2647 IF (LWRK1 .LT. 0) THEN 2648 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 2649 CALL QUIT('Insufficient work space in routine DIAC22') 2650 ENDIF 2651C 2652 CALL DZERO(WORK(KT2SUB),NVIR(ISYMA)*NT1AM(ISYMDK)) 2653 CALL DZERO(WORK(KACRES),NVIR(ISYMA)*NVIR(ISYMC)) 2654C 2655C------------------------------------------ 2656C Copy submatrix out of t2amt. 2657C------------------------------------------ 2658C 2659 DO 130 NDK = 1,NT1AM(ISYMDK) 2660C 2661 DO 140 A = 1,NVIR(ISYMA) 2662C 2663 NAI = IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I - 1) + A 2664 NAIDK = IT2AM(ISYMAI,ISYMDK) + INDEX(NAI,NDK) 2665 NADK = KT2SUB + NVIR(ISYMA)*(NDK - 1) + A - 1 2666C 2667 WORK(NADK) = T2AMT(NAIDK) 2668C 2669 140 CONTINUE 2670 130 CONTINUE 2671C 2672C----------------------------------------------------------- 2673C Calculate contraction of t2amt(sub) and z2ao. 2674C----------------------------------------------------------- 2675C 2676 KOFF1 = ICKATR(ISYMDK,ISYMC) + 1 2677C 2678 NTOTA = MAX(NVIR(ISYMA),1) 2679 NTOTDK = MAX(NT1AM(ISYMDK),1) 2680C 2681 CALL DGEMM('N','N',NVIR(ISYMA),NVIR(ISYMC), 2682 * NT1AM(ISYMDK),ONE,WORK(KT2SUB),NTOTA, 2683 * Z2AO(KOFF1),NTOTDK,ZERO,WORK(KACRES),NTOTA) 2684C 2685C-------------------------------------- 2686C Final storage in result. 2687C-------------------------------------- 2688C 2689 DO 150 C = 1,NVIR(ISYMC) 2690C 2691 NAC = KACRES + NVIR(ISYMA)*(C - 1) 2692 NAIC = ICKATR(ISYMAI,ISYMC) + NT1AM(ISYMAI)*(C - 1) 2693 * + IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I - 1) + 1 2694C 2695 CALL DAXPY(NVIR(ISYMA),ONE,WORK(NAC),1,D2IAB(NAIC),1) 2696C 2697 150 CONTINUE 2698 120 CONTINUE 2699 110 CONTINUE 2700 100 CONTINUE 2701C 2702 CALL QEXIT('DIAC22') 2703C 2704 RETURN 2705 END 2706C /* Deck dabi22 */ 2707 SUBROUTINE DABI22(D2ABI,ISYABI,T2AM,Z2AO,ISYZ2AO,WORK,LWORK) 2708C 2709C Written by Asger Halkier 24/2 - 1996 2710C 2711C Version: 1.0 2712C 2713C Purpose: To calculate the second contribution to D2AIB 2714C from projection against doubles! 2715C 2716#include "implicit.h" 2717 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 2718 DIMENSION D2ABI(*), T2AM(*), Z2AO(*), WORK(LWORK) 2719#include "priunit.h" 2720#include "ccorb.h" 2721#include "ccsdsym.h" 2722C 2723 INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 2724C 2725 CALL QENTER('DABI22') 2726C 2727 DO 100 ISYMI = 1,NSYM 2728C 2729 ISYMAB = MULD2H(ISYMI,ISYABI) 2730C 2731 DO 110 I = 1,NRHF(ISYMI) 2732C 2733 DO 120 ISYMA = 1,NSYM 2734C 2735 ISYMB = MULD2H(ISYMA,ISYMAB) 2736 ISYMCK = MULD2H(ISYMA,ISYZ2AO) 2737 ISYMBI = MULD2H(ISYMA,ISYABI) 2738C 2739C------------------------------------ 2740C Work space allocation. 2741C------------------------------------ 2742C 2743 KT2SUB = 1 2744 KBARES = KT2SUB + NVIR(ISYMB)*NT1AM(ISYMCK) 2745 KEND1 = KBARES + NVIR(ISYMB)*NVIR(ISYMA) 2746 LWRK1 = LWORK - KEND1 2747C 2748 IF (LWRK1 .LT. 0) THEN 2749 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 2750 CALL QUIT('Insufficient work space in routine DABI22') 2751 ENDIF 2752C 2753 CALL DZERO(WORK(KT2SUB),NVIR(ISYMB)*NT1AM(ISYMCK)) 2754 CALL DZERO(WORK(KBARES),NVIR(ISYMB)*NVIR(ISYMA)) 2755C 2756C------------------------------------------ 2757C Copy submatrix out of t2amt. 2758C------------------------------------------ 2759C 2760 DO 130 NCK = 1,NT1AM(ISYMCK) 2761C 2762 DO 140 B = 1,NVIR(ISYMB) 2763C 2764 NBI = IT1AM(ISYMB,ISYMI) + NVIR(ISYMB)*(I - 1) + B 2765 NBICK = IT2AM(ISYMBI,ISYMCK) + INDEX(NBI,NCK) 2766 NBCK = KT2SUB + NVIR(ISYMB)*(NCK - 1) + B - 1 2767C 2768 WORK(NBCK) = T2AM(NBICK) 2769C 2770 140 CONTINUE 2771 130 CONTINUE 2772C 2773C---------------------------------------------------------- 2774C Calculate contraction of t2am(sub) and z2ao. 2775C---------------------------------------------------------- 2776C 2777 KOFF1 = ICKATR(ISYMCK,ISYMA) + 1 2778C 2779 NTOTB = MAX(NVIR(ISYMB),1) 2780 NTOTCK = MAX(NT1AM(ISYMCK),1) 2781C 2782 CALL DGEMM('N','N',NVIR(ISYMB),NVIR(ISYMA), 2783 * NT1AM(ISYMCK),ONE,WORK(KT2SUB),NTOTB, 2784 * Z2AO(KOFF1),NTOTCK,ZERO,WORK(KBARES),NTOTB) 2785C 2786C-------------------------------------- 2787C Final storage in result. 2788C-------------------------------------- 2789C 2790 DO 150 A = 1,NVIR(ISYMA) 2791C 2792 NBA = KBARES + NVIR(ISYMB)*(A - 1) 2793 NAB = ICKASR(ISYMAB,ISYMI) + NMATAB(ISYMAB)*(I - 1) 2794 * + IMATAB(ISYMA,ISYMB) + A 2795C 2796 CALL DAXPY(NVIR(ISYMB),-ONE,WORK(NBA),1,D2ABI(NAB), 2797 * NVIR(ISYMA)) 2798C 2799 150 CONTINUE 2800 120 CONTINUE 2801 110 CONTINUE 2802 100 CONTINUE 2803C 2804 CALL QEXIT('DABI22') 2805C 2806 RETURN 2807 END 2808C /* Deck dabi23 */ 2809 SUBROUTINE DABI23(D2ABI,ISYABI,T2AM,Z2AO,ISYZ2AO,WORK,LWORK) 2810C 2811C Written by Asger Halkier 24/2 - 1996 2812C 2813C Version: 1.0 2814C 2815C Purpose: To calculate the third contribution to D2AIB 2816C from projection against doubles! 2817C 2818#include "implicit.h" 2819 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 2820 DIMENSION D2ABI(*), T2AM(*), Z2AO(*), WORK(LWORK) 2821#include "priunit.h" 2822#include "ccorb.h" 2823#include "ccsdsym.h" 2824C 2825 INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 2826C 2827 CALL QENTER('DABI23') 2828C 2829 DO 100 ISYMA = 1,NSYM 2830C 2831 ISYMCK = MULD2H(ISYMA,ISYZ2AO) 2832 ISYMBI = ISYMCK 2833C 2834C---------------------------------- 2835C Work space allocation one. 2836C---------------------------------- 2837C 2838 KZREOR = 1 2839 KEND1 = KZREOR + NVIR(ISYMA)*NT1AM(ISYMCK) 2840 LWRK1 = LWORK - KEND1 2841C 2842 IF (LWRK1 .LT. 0) THEN 2843 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 2844 CALL QUIT('Insufficient work space in routine DABI23') 2845 ENDIF 2846C 2847 CALL DZERO(WORK(KZREOR),NVIR(ISYMA)*NT1AM(ISYMCK)) 2848C 2849C------------------------------------------------- 2850C Reorder z2ao(ak,c;del) to z2ao(a,ck;del). 2851C------------------------------------------------- 2852C 2853 DO 110 ISYMC = 1,NSYM 2854C 2855 ISYMK = MULD2H(ISYMC,ISYMCK) 2856 ISYMAK = MULD2H(ISYMC,ISYZ2AO) 2857C 2858 DO 120 C = 1,NVIR(ISYMC) 2859C 2860 DO 130 K = 1,NRHF(ISYMK) 2861C 2862 NAKC = ICKATR(ISYMAK,ISYMC) + NT1AM(ISYMAK)*(C - 1) 2863 * + IT1AM(ISYMA,ISYMK) + NVIR(ISYMA)*(K - 1) + 1 2864 NCK = IT1AM(ISYMC,ISYMK) + NVIR(ISYMC)*(K - 1) + C 2865 NACK = KZREOR + NVIR(ISYMA)*(NCK - 1) 2866C 2867 CALL DCOPY(NVIR(ISYMA),Z2AO(NAKC),1,WORK(NACK),1) 2868C 2869 130 CONTINUE 2870 120 CONTINUE 2871 110 CONTINUE 2872C 2873 DO 140 ISYMI = 1,NSYM 2874C 2875 ISYMB = MULD2H(ISYMI,ISYMBI) 2876 ISYMAB = MULD2H(ISYMI,ISYABI) 2877C 2878C------------------------------------- 2879C Work space allocation two. 2880C------------------------------------- 2881C 2882 KT2SUB = KEND1 2883 KEND2 = KT2SUB + NT1AM(ISYMCK)*NVIR(ISYMB) 2884 LWRK2 = LWORK - KEND2 2885C 2886 IF (LWRK2 .LT. 0) THEN 2887 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND2 2888 CALL QUIT('Insufficient work space in routine DABI23') 2889 ENDIF 2890C 2891 DO 150 I = 1,NRHF(ISYMI) 2892C 2893 CALL DZERO(WORK(KT2SUB),NT1AM(ISYMCK)*NVIR(ISYMB)) 2894C 2895C----------------------------------------- 2896C Copy submatrix out of t2am. 2897C----------------------------------------- 2898C 2899 DO 160 ISYMK = 1,NSYM 2900C 2901 ISYMC = MULD2H(ISYMK,ISYMCK) 2902 ISYMCI = MULD2H(ISYMC,ISYMI) 2903 ISYMBK = ISYMCI 2904C 2905 DO 170 K = 1,NRHF(ISYMK) 2906C 2907 DO 180 B = 1,NVIR(ISYMB) 2908C 2909 NBK = IT1AM(ISYMB,ISYMK) + NVIR(ISYMB)*(K - 1) 2910 * + B 2911C 2912 DO 190 C = 1,NVIR(ISYMC) 2913C 2914 NCI = IT1AM(ISYMC,ISYMI) 2915 * + NVIR(ISYMC)*(I - 1) + C 2916 NCIBK = IT2AM(ISYMCI,ISYMBK) 2917 * + INDEX(NCI,NBK) 2918 NCKB = KT2SUB + NT1AM(ISYMCK)*(B - 1) 2919 * + IT1AM(ISYMC,ISYMK) 2920 * + NVIR(ISYMC)*(K - 1) + C - 1 2921C 2922 WORK(NCKB) = T2AM(NCIBK) 2923C 2924 190 CONTINUE 2925 180 CONTINUE 2926 170 CONTINUE 2927 160 CONTINUE 2928C 2929C------------------------------------------------------ 2930C Final contraction and storage in result. 2931C------------------------------------------------------ 2932C 2933 KOFF1 = ICKASR(ISYMAB,ISYMI) + NMATAB(ISYMAB)*(I - 1) 2934 * + IMATAB(ISYMA,ISYMB) + 1 2935C 2936 NTOTA = MAX(NVIR(ISYMA),1) 2937 NTOTCK = MAX(NT1AM(ISYMCK),1) 2938C 2939 CALL DGEMM('N','N',NVIR(ISYMA),NVIR(ISYMB), 2940 * NT1AM(ISYMCK),-ONE,WORK(KZREOR),NTOTA, 2941 * WORK(KT2SUB),NTOTCK,ONE,D2ABI(KOFF1),NTOTA) 2942C 2943 150 CONTINUE 2944 140 CONTINUE 2945 100 CONTINUE 2946C 2947 CALL QEXIT('DABI23') 2948C 2949 RETURN 2950 END 2951C /* Deck diaj26 */ 2952 SUBROUTINE DIAJ26(D2IAJ,ISYIAJ,T2AM,ZINT,ISYMZI,WORK,LWORK) 2953C 2954C Written by Asger Halkier 25/2 - 1996 2955C 2956C Version: 1.0 2957C 2958C Purpose: To calculate the sixth contribution to D2IAJ 2959C from projection against doubles! 2960C 2961#include "implicit.h" 2962 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 2963 DIMENSION D2IAJ(*), T2AM(*), ZINT(*), WORK(LWORK) 2964#include "priunit.h" 2965#include "ccorb.h" 2966#include "ccsdsym.h" 2967C 2968 INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 2969C 2970 CALL QENTER('DIAJ26') 2971C 2972 DO 100 ISYMJ = 1,NSYM 2973C 2974 ISYMAI = MULD2H(ISYMJ,ISYIAJ) 2975C 2976 DO 110 J = 1,NRHF(ISYMJ) 2977C 2978 DO 120 ISYMA = 1,NSYM 2979C 2980 ISYMI = MULD2H(ISYMA,ISYMAI) 2981 ISYMEM = MULD2H(ISYMI,ISYMZI) 2982 ISYMAJ = ISYMEM 2983C 2984C------------------------------------ 2985C Work space allocation. 2986C------------------------------------ 2987C 2988 KT2SUB = 1 2989 KEND1 = KT2SUB + NVIR(ISYMA)*NT1AM(ISYMEM) 2990 LWRK1 = LWORK - KEND1 2991C 2992 IF (LWRK1 .LT. 0) THEN 2993 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 2994 CALL QUIT('Insufficient work space in routine DIAJ26') 2995 ENDIF 2996C 2997 CALL DZERO(WORK(KT2SUB),NVIR(ISYMA)*NT1AM(ISYMEM)) 2998C 2999C------------------------------------------------- 3000C Copy submatrix out of t2-amplitude. 3001C------------------------------------------------- 3002C 3003 DO 130 NEM = 1,NT1AM(ISYMEM) 3004C 3005 DO 140 A = 1,NVIR(ISYMA) 3006C 3007 NAJ = IT1AM(ISYMA,ISYMJ) + NVIR(ISYMA)*(J - 1) + A 3008 NAJEM = IT2AM(ISYMAJ,ISYMEM) + INDEX(NAJ,NEM) 3009 NAEM = KT2SUB + NVIR(ISYMA)*(NEM - 1) + A - 1 3010C 3011 WORK(NAEM) = T2AM(NAJEM) 3012C 3013 140 CONTINUE 3014 130 CONTINUE 3015C 3016C------------------------------------------------------ 3017C Final contraction and storage in result. 3018C------------------------------------------------------ 3019C 3020 KOFF1 = IT2BCD(ISYMEM,ISYMI) + 1 3021 KOFF2 = IT2BCD(ISYMAI,ISYMJ) + NT1AM(ISYMAI)*(J - 1) 3022 * + IT1AM(ISYMA,ISYMI) + 1 3023C 3024 NTOTA = MAX(NVIR(ISYMA),1) 3025 NTOTEM = MAX(NT1AM(ISYMEM),1) 3026C 3027 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI), 3028 * NT1AM(ISYMEM),-ONE,WORK(KT2SUB),NTOTA, 3029 * ZINT(KOFF1),NTOTEM,ONE,D2IAJ(KOFF2),NTOTA) 3030C 3031 120 CONTINUE 3032 110 CONTINUE 3033 100 CONTINUE 3034C 3035 CALL QEXIT('DIAJ26') 3036C 3037 RETURN 3038 END 3039C /* Deck diaj27 */ 3040 SUBROUTINE DIAJ27(D2IAJ,ISYIAJ,T2AM,ZINT,ISYMZI,WORK,LWORK) 3041C 3042C Written by Asger Halkier 25/2 - 1996 3043C 3044C Version: 1.0 3045C 3046C Purpose: To calculate the seventh contribution to D2IAJ 3047C from projection against doubles! 3048C 3049#include "implicit.h" 3050 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 3051 DIMENSION D2IAJ(*), T2AM(*), ZINT(*), WORK(LWORK) 3052#include "priunit.h" 3053#include "ccorb.h" 3054#include "ccsdsym.h" 3055C 3056 INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 3057C 3058 CALL QENTER('DIAJ27') 3059C 3060 IF (ISYIAJ.NE.ISYMZI) CALL QUIT('Symmetry mismatch in DIAJ27') 3061C 3062 DO 100 ISYMJ = 1,NSYM 3063C 3064 ISYMAI = MULD2H(ISYMJ,ISYIAJ) 3065C 3066 DO 110 J = 1,NRHF(ISYMJ) 3067C 3068 DO 120 ISYMA = 1,NSYM 3069C 3070 ISYMI = MULD2H(ISYMA,ISYMAI) 3071 ISYMEM = MULD2H(ISYMI,ISYMZI) 3072 ISYMAJ = ISYMEM 3073C 3074C------------------------------------ 3075C Work space allocation. 3076C------------------------------------ 3077C 3078 KT2SUB = 1 3079 KEND1 = KT2SUB + NVIR(ISYMA)*NT1AM(ISYMEM) 3080 LWRK1 = LWORK - KEND1 3081C 3082 IF (LWRK1 .LT. 0) THEN 3083 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 3084 CALL QUIT('Insufficient work space in routine DIAJ27') 3085 ENDIF 3086C 3087 CALL DZERO(WORK(KT2SUB),NVIR(ISYMA)*NT1AM(ISYMEM)) 3088C 3089C------------------------------------------------- 3090C Copy submatrix out of t2-amplitude. 3091C------------------------------------------------- 3092C 3093 DO 130 ISYME = 1,NSYM 3094C 3095 ISYMM = MULD2H(ISYME,ISYMEM) 3096 ISYMAM = MULD2H(ISYMM,ISYMA) 3097 ISYMEJ = ISYMAM 3098C 3099 DO 140 E = 1,NVIR(ISYME) 3100C 3101 NEJ = IT1AM(ISYME,ISYMJ) + NVIR(ISYME)*(J - 1) + E 3102C 3103 DO 150 M = 1,NRHF(ISYMM) 3104C 3105 NEM = IT1AM(ISYME,ISYMM) 3106 * + NVIR(ISYME)*(M - 1) + E 3107C 3108 DO 160 A = 1,NVIR(ISYMA) 3109C 3110 NAM = IT1AM(ISYMA,ISYMM) 3111 * + NVIR(ISYMA)*(M - 1) + A 3112 NAMEJ = IT2AM(ISYMAM,ISYMEJ) 3113 * + INDEX(NAM,NEJ) 3114 NAEM = KT2SUB + NVIR(ISYMA)*(NEM - 1) 3115 * + A - 1 3116C 3117 WORK(NAEM) = T2AM(NAMEJ) 3118C 3119 160 CONTINUE 3120 150 CONTINUE 3121 140 CONTINUE 3122 130 CONTINUE 3123C 3124C------------------------------------------------------ 3125C Final contraction and storage in result. 3126C------------------------------------------------------ 3127C 3128 KOFF1 = IT2BCD(ISYMEM,ISYMI) + 1 3129 KOFF2 = IT2BCD(ISYMAI,ISYMJ) + NT1AM(ISYMAI)*(J - 1) 3130 * + IT1AM(ISYMA,ISYMI) + 1 3131C 3132 NTOTA = MAX(NVIR(ISYMA),1) 3133 NTOTEM = MAX(NT1AM(ISYMEM),1) 3134C 3135 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI), 3136 * NT1AM(ISYMEM),ONE,WORK(KT2SUB),NTOTA, 3137 * ZINT(KOFF1),NTOTEM,ONE,D2IAJ(KOFF2),NTOTA) 3138C 3139 120 CONTINUE 3140 110 CONTINUE 3141 100 CONTINUE 3142C 3143 CALL QEXIT('DIAJ27') 3144C 3145 RETURN 3146 END 3147C /* Deck diaj28 */ 3148 SUBROUTINE DIAJ28(D2IAJ,ISYIAJ,T2AMT,ZINT,ISYMZI,WORK,LWORK) 3149C 3150C Written by Asger Halkier 26/2 - 1996 3151C 3152C Version: 1.0 3153C 3154C Purpose: To calculate the eighth contribution to D2IAJ 3155C from projection against doubles! 3156C 3157#include "implicit.h" 3158 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 3159 DIMENSION D2IAJ(*), T2AMT(*), ZINT(*), WORK(LWORK) 3160#include "priunit.h" 3161#include "ccorb.h" 3162#include "ccsdsym.h" 3163C 3164 INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 3165C 3166 CALL QENTER('DIAJ28') 3167C 3168 DO 100 ISYMI = 1,NSYM 3169C 3170 ISYMAJ = MULD2H(ISYMI,ISYIAJ) 3171C 3172 DO 110 I = 1,NRHF(ISYMI) 3173C 3174 DO 120 ISYMA = 1,NSYM 3175C 3176 ISYMJ = MULD2H(ISYMA,ISYMAJ) 3177 ISYMEM = MULD2H(ISYMJ,ISYMZI) 3178 ISYMAI = ISYMEM 3179C 3180C------------------------------------ 3181C Work space allocation. 3182C------------------------------------ 3183C 3184 KT2SUB = 1 3185 KEND1 = KT2SUB + NVIR(ISYMA)*NT1AM(ISYMEM) 3186 LWRK1 = LWORK - KEND1 3187C 3188 IF (LWRK1 .LT. 0) THEN 3189 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 3190 CALL QUIT('Insufficient work space in routine DIAJ28') 3191 ENDIF 3192C 3193 CALL DZERO(WORK(KT2SUB),NVIR(ISYMA)*NT1AM(ISYMEM)) 3194C 3195C------------------------------------------ 3196C Copy submatrix out of t2amt. 3197C------------------------------------------ 3198C 3199 DO 130 NEM = 1,NT1AM(ISYMEM) 3200C 3201 DO 140 A = 1,NVIR(ISYMA) 3202C 3203 NAI = IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I - 1) + A 3204 NAIEM = IT2AM(ISYMAI,ISYMEM) + INDEX(NAI,NEM) 3205 NAEM = KT2SUB + NVIR(ISYMA)*(NEM - 1) + A - 1 3206C 3207 WORK(NAEM) = T2AMT(NAIEM) 3208C 3209 140 CONTINUE 3210 130 CONTINUE 3211C 3212C--------------------------------------------------------- 3213C Calculate final contraction in loop over j. 3214C--------------------------------------------------------- 3215C 3216 DO 150 J = 1,NRHF(ISYMJ) 3217C 3218 KOFF1 = IT2BCD(ISYMEM,ISYMJ) 3219 * + NT1AM(ISYMEM)*(J - 1) + 1 3220 KOFF2 = IT2BCD(ISYMAI,ISYMJ) 3221 * + NT1AM(ISYMAI)*(J - 1) 3222 * + IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I - 1) + 1 3223C 3224 NTOTA = MAX(NVIR(ISYMA),1) 3225C 3226 CALL DGEMV('N',NVIR(ISYMA),NT1AM(ISYMEM),ONE, 3227 * WORK(KT2SUB),NTOTA,ZINT(KOFF1),1,ONE, 3228 * D2IAJ(KOFF2),1) 3229C 3230 150 CONTINUE 3231C 3232 120 CONTINUE 3233 110 CONTINUE 3234 100 CONTINUE 3235C 3236 CALL QEXIT('DIAJ28') 3237C 3238 RETURN 3239 END 3240C /* Deck cc_mirs */ 3241 SUBROUTINE CC_MIRS(XMIRES,XMINT) 3242C 3243C Written by Asger Halkier 26/2 - 1996 3244C 3245C Version: 1.0 3246C 3247C Purpose: To resort the M-intermediate XMINT (im,j;k) 3248C to (mk,i;j) and store the result in XMIRES! 3249C 3250#include "implicit.h" 3251 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 3252 DIMENSION XMIRES(*), XMINT(*) 3253#include "priunit.h" 3254#include "ccorb.h" 3255#include "ccsdsym.h" 3256C 3257 CALL QENTER('CC_MIRS') 3258C 3259 CALL DZERO(XMIRES,N3ORHF(1)) 3260C 3261C---------------------------------------------- 3262C Do the reordering through loads of loops. 3263C---------------------------------------------- 3264C 3265 DO 100 ISYMK = 1,NSYM 3266C 3267 ISYIMJ = ISYMK 3268C 3269 DO 110 K = 1,NRHF(ISYMK) 3270C 3271 DO 120 ISYMJ = 1,NSYM 3272C 3273 ISYMIM = MULD2H(ISYMJ,ISYIMJ) 3274 ISYMKI = ISYMJ 3275C 3276 DO 130 J = 1,NRHF(ISYMJ) 3277C 3278 DO 140 ISYMM = 1,NSYM 3279C 3280 ISYMI = MULD2H(ISYMM,ISYMIM) 3281 ISYMMK = MULD2H(ISYMI,ISYMKI) 3282C 3283 DO 150 M = 1,NRHF(ISYMM) 3284C 3285 NIMJK = I3ORHF(ISYIMJ,ISYMK) 3286 * + NMAIJK(ISYIMJ)*(K - 1) 3287 * + IMAIJK(ISYMIM,ISYMJ) 3288 * + NMATIJ(ISYMIM)*(J - 1) 3289 * + IMATIJ(ISYMI,ISYMM) 3290 * + NRHF(ISYMI)*(M - 1) + 1 3291 NMK = IMATIJ(ISYMM,ISYMK) 3292 * + NRHF(ISYMM)*(K - 1) + M 3293 NMKIJ = I3ORHF(ISYMKI,ISYMJ) 3294 * + NMAIJK(ISYMKI)*(J - 1) 3295 * + IMAIJK(ISYMMK,ISYMI) + NMK 3296C 3297 CALL DCOPY(NRHF(ISYMI),XMINT(NIMJK),1, 3298 * XMIRES(NMKIJ),NMATIJ(ISYMMK)) 3299C 3300 150 CONTINUE 3301 140 CONTINUE 3302 130 CONTINUE 3303 120 CONTINUE 3304 110 CONTINUE 3305 100 CONTINUE 3306C 3307 CALL QEXIT('CC_MIRS') 3308C 3309 RETURN 3310 END 3311C /* Deck diaj25 */ 3312 SUBROUTINE DIAJ25(D2IAJ,ISYIAJ,T2INT,ISYMTI,XMIRES,WORK,LWORK) 3313C 3314C Written by Asger Halkier 26/2 - 1996 3315C 3316C Version: 1.0 3317C 3318C Purpose: To calculate the fifth contribution to D2IAJ 3319C from projection against doubles! 3320C 3321#include "implicit.h" 3322 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 3323 DIMENSION D2IAJ(*), T2INT(*), XMIRES(*), WORK(LWORK) 3324#include "priunit.h" 3325#include "ccorb.h" 3326#include "ccsdsym.h" 3327C 3328 CALL QENTER('DIAJ25') 3329C 3330 IF (ISYMTI.NE.ISYIAJ) CALL QUIT('Symmetry mismatch in DIAJ25') 3331C 3332 DO 100 ISYMA = 1,NSYM 3333C 3334 ISYMMK = MULD2H(ISYMA,ISYMTI) 3335 ISYMIJ = MULD2H(ISYMA,ISYIAJ) 3336C 3337C------------------------------ 3338C Work space allocation. 3339C------------------------------ 3340C 3341 KTREOR = 1 3342 KEND1 = KTREOR + NVIR(ISYMA)*NMATIJ(ISYMMK) 3343 LWRK1 = LWORK - KEND1 3344C 3345 IF (LWRK1 .LT. 0) THEN 3346 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 3347 CALL QUIT('Insufficient work space in routine DIAJ25') 3348 ENDIF 3349C 3350 CALL DZERO(WORK(KTREOR),NVIR(ISYMA)*NMATIJ(ISYMMK)) 3351C 3352C-------------------------------------------- 3353C Reorder t2int in scratch-work-array. 3354C-------------------------------------------- 3355C 3356 DO 110 ISYMM = 1,NSYM 3357C 3358 ISYMK = MULD2H(ISYMM,ISYMMK) 3359 ISYMAM = MULD2H(ISYMK,ISYMTI) 3360C 3361 DO 120 K = 1,NRHF(ISYMK) 3362C 3363 DO 130 M = 1,NRHF(ISYMM) 3364C 3365 NMK = IMATIJ(ISYMM,ISYMK) + NRHF(ISYMM)*(K - 1) + M 3366 NAMKN = IT2BCD(ISYMAM,ISYMK) + NT1AM(ISYMAM)*(K - 1) 3367 * + IT1AM(ISYMA,ISYMM) + NVIR(ISYMA)*(M - 1) + 1 3368 NAMKR = KTREOR + NVIR(ISYMA)*(NMK - 1) 3369C 3370 CALL DCOPY(NVIR(ISYMA),T2INT(NAMKN),1,WORK(NAMKR),1) 3371C 3372 130 CONTINUE 3373 120 CONTINUE 3374 110 CONTINUE 3375C 3376C--------------------------------------------------------- 3377C Calculate contraction with xmires in loop over j. 3378C--------------------------------------------------------- 3379C 3380 DO 140 ISYMJ = 1,NSYM 3381C 3382 ISYMI = MULD2H(ISYMJ,ISYMIJ) 3383 ISYMAI = MULD2H(ISYMJ,ISYIAJ) 3384 ISYMKI = ISYMJ 3385C 3386 DO 150 J = 1,NRHF(ISYMJ) 3387C 3388 KOFF1 = I3ORHF(ISYMKI,ISYMJ) 3389 * + NMAIJK(ISYMKI)*(J - 1) 3390 * + IMAIJK(ISYMMK,ISYMI) + 1 3391 KOFF2 = IT2BCD(ISYMAI,ISYMJ) + NT1AM(ISYMAI)*(J - 1) 3392 * + IT1AM(ISYMA,ISYMI) + 1 3393C 3394 NTOTA = MAX(NVIR(ISYMA),1) 3395 NTOTMK = MAX(NMATIJ(ISYMMK),1) 3396C 3397 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI), 3398 * NMATIJ(ISYMMK),ONE,WORK(KTREOR),NTOTA, 3399 * XMIRES(KOFF1),NTOTMK,ONE,D2IAJ(KOFF2),NTOTA) 3400C 3401 150 CONTINUE 3402 140 CONTINUE 3403 100 CONTINUE 3404C 3405 CALL QEXIT('DIAJ25') 3406C 3407 RETURN 3408 END 3409C /* Deck dabgao */ 3410 SUBROUTINE DABGAO(D2ABG,ISYDEN,T2INT,ISYMTI,Z2INT,ISYMZI, 3411 * WORK,LWORK,G,ISYMG) 3412C 3413C Written by Asger Halkier 27/2 - 1996 3414C 3415C Version: 1.0 3416C 3417C Purpose: To calculate the only contribution from D2ABC 3418C (originating from projection against doubles) 3419C to D2ABG directly! 3420C 3421#if defined (IMPLICIT_NONE) 3422 IMPLICIT NONE 3423#else 3424# include "implicit.h" 3425#endif 3426#include "priunit.h" 3427#include "ccorb.h" 3428#include "ccsdsym.h" 3429C 3430#if defined (SYS_CRAY) 3431 REAL ZERO, ONE 3432#else 3433 DOUBLE PRECISION ZERO, ONE 3434#endif 3435 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 3436 3437 INTEGER ISYDEN, ISYMTI, ISYMZI, ISYMG, LWORK 3438 3439#if defined (SYS_CRAY) 3440 REAL D2ABG(*), T2INT(*), Z2INT(*), WORK(LWORK) 3441#else 3442 DOUBLE PRECISION D2ABG(*), T2INT(*), Z2INT(*), WORK(LWORK) 3443#endif 3444 3445 INTEGER ISYMA,ISYMAB,ISYMB,ISYMM, ISYMAM, ISYMBM, ISYMMN, ISYMN, 3446 & KZREOR, KTREOR, KEND1, LWRK1, kOFF1, NTOTA, NTOTMN, 3447 & NMN, NAMNN, NAMNR, NBMN, NMNB 3448C 3449 CALL QENTER('DABGAO') 3450C 3451 IF (MULD2H(ISYMTI,ISYMZI).NE.MULD2H(ISYDEN,ISYMG)) THEN 3452 WRITE(LUPRI,*) 'ISYMTI,ISYMZI,ISYMG,ISYDEN:', 3453 & ISYMTI,ISYMZI,ISYMG,ISYDEN 3454 CALL QUIT('Symmetry mismatch in DABGAO. (1)') 3455 END IF 3456C 3457 ISYMAB = MULD2H(ISYDEN,ISYMG) 3458C 3459 DO 100 ISYMA = 1,NSYM 3460C 3461 ISYMB = MULD2H(ISYMA,ISYMAB) 3462 ISYMMN = MULD2H(ISYMA,ISYMZI) 3463 IF (MULD2H(ISYMB,ISYMTI).NE.ISYMMN) THEN 3464 CALL QUIT('Symmetry mismatch in DABGAO. (2)') 3465 END IF 3466C 3467C------------------------------ 3468C Work space allocation. 3469C------------------------------ 3470C 3471 KZREOR = 1 3472 KTREOR = KZREOR + NVIR(ISYMA)*NMATIJ(ISYMMN) 3473 KEND1 = KTREOR + NMATIJ(ISYMMN)*NVIR(ISYMB) 3474 LWRK1 = LWORK - KEND1 3475C 3476 IF (LWRK1 .LT. 0) THEN 3477 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 3478 CALL QUIT('Insufficient work space in routine DABGAO') 3479 ENDIF 3480C 3481 CALL DZERO(WORK(KZREOR),NVIR(ISYMA)*NMATIJ(ISYMMN)) 3482 CALL DZERO(WORK(KTREOR),NMATIJ(ISYMMN)*NVIR(ISYMB)) 3483C 3484 DO 110 ISYMM = 1,NSYM 3485C 3486 ISYMN = MULD2H(ISYMM,ISYMMN) 3487C 3488C----------------------------------------------- 3489C Reorder z2int in scratch-work-array. 3490C----------------------------------------------- 3491C 3492 ISYMAM = MULD2H(ISYMN,ISYMZI) 3493 3494 DO 120 N = 1,NRHF(ISYMN) 3495C 3496 DO 130 M = 1,NRHF(ISYMM) 3497C 3498 NMN = IMATIJ(ISYMM,ISYMN) + NRHF(ISYMM)*(N - 1) + M 3499 NAMNN = IT2BCD(ISYMAM,ISYMN) + NT1AM(ISYMAM)*(N - 1) 3500 * + IT1AM(ISYMA,ISYMM) + NVIR(ISYMA)*(M - 1) + 1 3501 NAMNR = KZREOR + NVIR(ISYMA)*(NMN - 1) 3502C 3503 CALL DCOPY(NVIR(ISYMA),Z2INT(NAMNN),1,WORK(NAMNR),1) 3504C 3505 130 CONTINUE 3506 120 CONTINUE 3507C 3508C----------------------------------------------- 3509C Reorder t2int in scratch-work-array. 3510C----------------------------------------------- 3511C 3512 ISYMBM = MULD2H(ISYMN,ISYMTI) 3513 3514 DO 140 N = 1,NRHF(ISYMN) 3515C 3516 DO 150 M = 1,NRHF(ISYMM) 3517C 3518 NMN = IMATIJ(ISYMM,ISYMN) + NRHF(ISYMM)*(N - 1) + M 3519 NBMN = IT2BCD(ISYMBM,ISYMN) + NT1AM(ISYMBM)*(N - 1) 3520 * + IT1AM(ISYMB,ISYMM) + NVIR(ISYMB)*(M - 1) + 1 3521 NMNB = KTREOR + NMN - 1 3522C 3523 CALL DCOPY(NVIR(ISYMB),T2INT(NBMN),1,WORK(NMNB), 3524 * NMATIJ(ISYMMN)) 3525C 3526 150 CONTINUE 3527 140 CONTINUE 3528 110 CONTINUE 3529C 3530C------------------------------------------------------------------ 3531C Contraction of two scratch-matrices and storage in result. 3532C------------------------------------------------------------------ 3533C 3534 KOFF1 = ID2ABG(ISYMAB,ISYMG) + NMATAB(ISYMAB)*(G - 1) 3535 * + IMATAB(ISYMA,ISYMB) + 1 3536C 3537 NTOTA = MAX(NVIR(ISYMA),1) 3538 NTOTMN = MAX(NMATIJ(ISYMMN),1) 3539C 3540 CALL DGEMM('N','N',NVIR(ISYMA),NVIR(ISYMB),NMATIJ(ISYMMN), 3541 * ONE,WORK(KZREOR),NTOTA,WORK(KTREOR),NTOTMN,ONE, 3542 * D2ABG(KOFF1),NTOTA) 3543C 3544 100 CONTINUE 3545C 3546 CALL QEXIT('DABGAO') 3547C 3548 RETURN 3549 END 3550C /* Deck cc1intmo */ 3551 SUBROUTINE CCDINTMO(XINTIJ,XINTIA,XINTAB,XINTAI,XAOINT, 3552 * XLAMDP,XLAMDH,WORK,LWORK,ISYM) 3553C 3554C Written by Asger Halkier 19/3 - 1996 3555C 3556C Version: 1.0 3557C 3558C Purpose: To transform the AO integrals in (XAOINT) 3559C to MO-basis (stored in XINTIJ through XINTAI)! 3560C Note that the AO integrals either can be the one 3561C electron integrals or a submatrix alfa,beta of the two 3562C electron integrals (alfa beta|gamma delta)! 3563C ISYM is the integral symmetry! 3564C 3565#include "implicit.h" 3566 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 3567 DIMENSION XINTIJ(*), XINTIA(*), XINTAB(*), XINTAI(*) 3568 DIMENSION XAOINT(*), XLAMDP(*), XLAMDH(*), WORK(LWORK) 3569#include "priunit.h" 3570#include "ccorb.h" 3571#include "ccsdsym.h" 3572C 3573 CALL QENTER('CCDINTMO') 3574C 3575C------------------------------ 3576C Initialize result arrays. 3577C------------------------------ 3578C 3579 CALL DZERO(XINTIJ,NMATIJ(ISYM)) 3580 CALL DZERO(XINTIA,NT1AM(ISYM)) 3581 CALL DZERO(XINTAB,NMATAB(ISYM)) 3582 CALL DZERO(XINTAI,NT1AM(ISYM)) 3583C 3584 DO 100 ISYMAL = 1,NSYM 3585C 3586 ISYMBE = MULD2H(ISYMAL,ISYM) 3587 ISYMP = ISYMAL 3588 ISYMQ = ISYMBE 3589C 3590C------------------------------ 3591C Work space allocation. 3592C------------------------------ 3593C 3594 LESCR = MAX(NBAS(ISYMAL)*NRHF(ISYMQ),NBAS(ISYMAL)*NVIR(ISYMQ)) 3595C 3596 KSCR = 1 3597 KEND1 = KSCR + LESCR 3598 LWRK1 = LWORK - KEND1 3599C 3600 IF (LWRK1 .LT. 0) THEN 3601 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 3602 CALL QUIT('Insufficient memory for allocation in CCDINTMO') 3603 ENDIF 3604C 3605 CALL DZERO(WORK(KSCR),NBAS(ISYMAL)*NRHF(ISYMQ)) 3606C 3607C---------------------------------------------------------- 3608C Transform second integral index to occupied space. 3609C---------------------------------------------------------- 3610C 3611 ISYMI = ISYMQ 3612C 3613 KOFF1 = IAODIS(ISYMAL,ISYMBE) + 1 3614 KOFF2 = IGLMRH(ISYMBE,ISYMI) + 1 3615C 3616 NTOTAL = MAX(NBAS(ISYMAL),1) 3617 NTOTBE = MAX(NBAS(ISYMBE),1) 3618C 3619 CALL DGEMM('N','N',NBAS(ISYMAL),NRHF(ISYMI),NBAS(ISYMBE),ONE, 3620 * XAOINT(KOFF1),NTOTAL,XLAMDH(KOFF2),NTOTBE,ZERO, 3621 * WORK(KSCR),NTOTAL) 3622C 3623C------------------------------------------ 3624C Calculate the (occ.occ) integrals. 3625C------------------------------------------ 3626C 3627 ISYMJ = ISYMP 3628C 3629 KOFF3 = IGLMRH(ISYMAL,ISYMJ) + 1 3630 KOFF4 = IMATIJ(ISYMJ,ISYMI) + 1 3631C 3632 NTOTAL = MAX(NBAS(ISYMAL),1) 3633 NTOTJ = MAX(NRHF(ISYMJ),1) 3634C 3635 CALL DGEMM('T','N',NRHF(ISYMJ),NRHF(ISYMI),NBAS(ISYMAL),ONE, 3636 * XLAMDP(KOFF3),NTOTAL,WORK(KSCR),NTOTAL,ZERO, 3637 * XINTIJ(KOFF4),NTOTJ) 3638C 3639C------------------------------------------ 3640C Calculate the (vir.occ) integrals. 3641C------------------------------------------ 3642C 3643 ISYMA = ISYMP 3644C 3645 KOFF5 = IGLMVI(ISYMAL,ISYMA) + 1 3646 KOFF6 = IT1AM(ISYMA,ISYMI) + 1 3647C 3648 NTOTAL = MAX(NBAS(ISYMAL),1) 3649 NTOTA = MAX(NVIR(ISYMA),1) 3650C 3651 CALL DGEMM('T','N',NVIR(ISYMA),NRHF(ISYMI),NBAS(ISYMAL),ONE, 3652 * XLAMDP(KOFF5),NTOTAL,WORK(KSCR),NTOTAL,ZERO, 3653 * XINTAI(KOFF6),NTOTA) 3654C 3655C--------------------------------------------------------- 3656C Transform second integral index to virtual space. 3657C--------------------------------------------------------- 3658C 3659 CALL DZERO(WORK(KSCR),NBAS(ISYMAL)*NVIR(ISYMQ)) 3660C 3661 ISYMA = ISYMQ 3662C 3663 KOFF7 = IAODIS(ISYMAL,ISYMBE) + 1 3664 KOFF8 = IGLMVI(ISYMBE,ISYMA) + 1 3665C 3666 NTOTAL = MAX(NBAS(ISYMAL),1) 3667 NTOTBE = MAX(NBAS(ISYMBE),1) 3668C 3669 CALL DGEMM('N','N',NBAS(ISYMAL),NVIR(ISYMA),NBAS(ISYMBE),ONE, 3670 * XAOINT(KOFF7),NTOTAL,XLAMDH(KOFF8),NTOTBE,ZERO, 3671 * WORK(KSCR),NTOTAL) 3672C 3673C------------------------------------------ 3674C Calculate the (occ.vir) integrals. 3675C Note that these are stored trans- 3676C posed, i.e. (vir.occ) like a t1! 3677C------------------------------------------ 3678C 3679 ISYMI = ISYMP 3680C 3681 KOFF9 = IGLMRH(ISYMAL,ISYMI) + 1 3682 KOFF10 = IT1AM(ISYMA,ISYMI) + 1 3683C 3684 NTOTAL = MAX(NBAS(ISYMAL),1) 3685 NTOTA = MAX(NVIR(ISYMA),1) 3686C 3687 CALL DGEMM('T','N',NVIR(ISYMA),NRHF(ISYMI),NBAS(ISYMAL),ONE, 3688 * WORK(KSCR),NTOTAL,XLAMDP(KOFF9),NTOTAL,ZERO, 3689 * XINTIA(KOFF10),NTOTA) 3690C 3691C------------------------------------------ 3692C Calculate the (vir.vir) integrals. 3693C------------------------------------------ 3694C 3695 ISYMB = ISYMP 3696C 3697 KOFF11 = IGLMVI(ISYMAL,ISYMB) + 1 3698 KOFF12 = IMATAB(ISYMB,ISYMA) + 1 3699C 3700 NTOTAL = MAX(NBAS(ISYMAL),1) 3701 NTOTB = MAX(NVIR(ISYMB),1) 3702C 3703 CALL DGEMM('T','N',NVIR(ISYMB),NVIR(ISYMA),NBAS(ISYMAL),ONE, 3704 * XLAMDP(KOFF11),NTOTAL,WORK(KSCR),NTOTAL,ZERO, 3705 * XINTAB(KOFF12),NTOTB) 3706C 3707 100 CONTINUE 3708C 3709 CALL QEXIT('CCDINTMO') 3710C 3711 RETURN 3712 END 3713C /* Deck ccdenzk0 */ 3714 SUBROUTINE CCDENZK0(ETAKA,XINTIJ,XINTAI,XINTIA,XINTAB,DIJ,DAI, 3715 * DIA,DAB,XIJT,XABT,DIJT,DABT,T1AM, 3716 * WORK,LWORK,ISYM) 3717C 3718C Written by Asger Halkier 20/3 - 1996 3719C 3720C Version: 1.0 3721C 3722C Purpose: To set up the right hand side of the equation for 3723C zeta-kappa-0 (ETAKA) from MO-integrals (XI*), Coupled 3724C Cluster densities (D*) and t1-amplitudes (T1AM)! 3725C Note that due to the symmetry in the formulas, this 3726C routine is able to handle both the one- and the two- 3727C electron contributions! 3728C ISYM is the symmetry of both the density and the 3729C integrals! 3730C 3731#include "implicit.h" 3732 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 3733 DIMENSION ETAKA(*), XINTIJ(*), XINTAI(*), XINTIA(*), XINTAB(*) 3734 DIMENSION DIJ(*), DAI(*), DIA(*), DAB(*), XIJT(*), XABT(*) 3735 DIMENSION DIJT(*), DABT(*), T1AM(*), WORK(LWORK) 3736#include "priunit.h" 3737#include "ccorb.h" 3738#include "ccsdsym.h" 3739C 3740 CALL QENTER('CCDENZK0') 3741C 3742 DO 100 ISYMA = 1,NSYM 3743C 3744C------------------------------------------------------- 3745C Calculate terms originating from [H(t1),E(ai)]. 3746C------------------------------------------------------- 3747C 3748 ISYMI = ISYMA 3749 ISYMB = MULD2H(ISYMA,ISYM) 3750 ISYMJ = MULD2H(ISYMA,ISYM) 3751C 3752 KOFFRE = IT1AM(ISYMA,ISYMI) + 1 3753C 3754 NTOTRE = MAX(NVIR(ISYMA),1) 3755 NTOTI = MAX(NRHF(ISYMI),1) 3756 NTOTB = MAX(NVIR(ISYMB),1) 3757 NTOTJ = MAX(NRHF(ISYMJ),1) 3758C 3759 KOFF1 = IMATAB(ISYMA,ISYMB) + 1 3760 KOFF2 = IT1AM(ISYMB,ISYMI) + 1 3761C 3762 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI),NVIR(ISYMB), 3763 * ONE,XABT(KOFF1),NTOTRE,DAI(KOFF2),NTOTB,ONE, 3764 * ETAKA(KOFFRE),NTOTRE) 3765C 3766 KOFF3 = IMATAB(ISYMA,ISYMB) + 1 3767 KOFF4 = IT1AM(ISYMB,ISYMI) + 1 3768C 3769 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI),NVIR(ISYMB), 3770 * -ONE,DAB(KOFF3),NTOTRE,XINTIA(KOFF4),NTOTB,ONE, 3771 * ETAKA(KOFFRE),NTOTRE) 3772C 3773 KOFF5 = IT1AM(ISYMA,ISYMJ) + 1 3774 KOFF6 = IMATIJ(ISYMJ,ISYMI) + 1 3775C 3776 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI),NRHF(ISYMJ), 3777 * ONE,XINTIA(KOFF5),NTOTRE,DIJ(KOFF6),NTOTJ,ONE, 3778 * ETAKA(KOFFRE),NTOTRE) 3779C 3780 KOFF7 = IT1AM(ISYMA,ISYMJ) + 1 3781 KOFF8 = IMATIJ(ISYMJ,ISYMI) + 1 3782C 3783 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI),NRHF(ISYMJ), 3784 * -ONE,DAI(KOFF7),NTOTRE,XIJT(KOFF8),NTOTJ,ONE, 3785 * ETAKA(KOFFRE),NTOTRE) 3786C 3787C------------------------------------------------------- 3788C Calculate terms originating from [H(t1),E(ia)]. 3789C------------------------------------------------------- 3790C 3791 KOFF9 = IMATAB(ISYMA,ISYMB) + 1 3792 KOFF10 = IT1AM(ISYMB,ISYMI) + 1 3793C 3794 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI),NVIR(ISYMB), 3795 * -ONE,DABT(KOFF9),NTOTRE,XINTAI(KOFF10),NTOTB, 3796 * ONE,ETAKA(KOFFRE),NTOTRE) 3797C 3798 KOFF11 = IMATAB(ISYMA,ISYMB) + 1 3799 KOFF12 = IT1AM(ISYMB,ISYMI) + 1 3800C 3801 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI),NVIR(ISYMB), 3802 * ONE,XINTAB(KOFF11),NTOTRE,DIA(KOFF12),NTOTB, 3803 * ONE,ETAKA(KOFFRE),NTOTRE) 3804C 3805 KOFF13 = IT1AM(ISYMA,ISYMJ) + 1 3806 KOFF14 = IMATIJ(ISYMJ,ISYMI) + 1 3807C 3808 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI),NRHF(ISYMJ), 3809 * -ONE,DIA(KOFF13),NTOTRE,XINTIJ(KOFF14),NTOTJ, 3810 * ONE,ETAKA(KOFFRE),NTOTRE) 3811C 3812 KOFF15 = IT1AM(ISYMA,ISYMJ) + 1 3813 KOFF16 = IMATIJ(ISYMJ,ISYMI) + 1 3814C 3815 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI),NRHF(ISYMJ), 3816 * ONE,XINTAI(KOFF15),NTOTRE,DIJT(KOFF16),NTOTJ, 3817 * ONE,ETAKA(KOFFRE),NTOTRE) 3818C 3819 100 CONTINUE 3820C 3821C---------------------------------------------------- 3822C Calculate terms originating from [H(t1),E(ik)]. 3823C---------------------------------------------------- 3824C 3825 DO 110 ISYMA = 1,NSYM 3826C 3827 ISYMI = ISYMA 3828 ISYMK = ISYMA 3829 ISYMB = MULD2H(ISYMI,ISYM) 3830 ISYMJ = MULD2H(ISYMI,ISYM) 3831C 3832 KOFFR = IT1AM(ISYMA,ISYMI) + 1 3833 KOFFT = IT1AM(ISYMA,ISYMK) + 1 3834C 3835 NTOTR = MAX(NVIR(ISYMA),1) 3836 NTOTB = MAX(NVIR(ISYMB),1) 3837 NTOTK = MAX(NRHF(ISYMK),1) 3838 NTOTI = MAX(NRHF(ISYMI),1) 3839 NTOTJ = MAX(NRHF(ISYMJ),1) 3840C 3841C---------------------------------- 3842C Work space allocation one. 3843C---------------------------------- 3844C 3845 KSCR1 = 1 3846 KEND1 = KSCR1 + NRHF(ISYMK)*NRHF(ISYMI) 3847 LWRK1 = LWORK - KEND1 3848C 3849 IF (LWRK1 .LT. 0) THEN 3850 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 3851 CALL QUIT('Insufficient memory for work '// 3852 & 'allocation in CCDENZK0') 3853 ENDIF 3854C 3855 CALL DZERO(WORK(KSCR1),NRHF(ISYMK)*NRHF(ISYMI)) 3856C 3857 KOFF17 = IT1AM(ISYMB,ISYMK) + 1 3858 KOFF18 = IT1AM(ISYMB,ISYMI) + 1 3859C 3860 CALL DGEMM('T','N',NRHF(ISYMK),NRHF(ISYMI),NVIR(ISYMB),ONE, 3861 * XINTIA(KOFF17),NTOTB,DIA(KOFF18),NTOTB,ONE, 3862 * WORK(KSCR1),NTOTK) 3863C 3864 KOFF19 = IT1AM(ISYMB,ISYMK) + 1 3865 KOFF20 = IT1AM(ISYMB,ISYMI) + 1 3866C 3867 CALL DGEMM('T','N',NRHF(ISYMK),NRHF(ISYMI),NVIR(ISYMB),-ONE, 3868 * DAI(KOFF19),NTOTB,XINTAI(KOFF20),NTOTB,ONE, 3869 * WORK(KSCR1),NTOTK) 3870C 3871 KOFF21 = IMATIJ(ISYMK,ISYMJ) + 1 3872 KOFF22 = IMATIJ(ISYMJ,ISYMI) + 1 3873C 3874 CALL DGEMM('N','N',NRHF(ISYMK),NRHF(ISYMI),NRHF(ISYMJ),ONE, 3875 * XINTIJ(KOFF21),NTOTK,DIJT(KOFF22),NTOTJ,ONE, 3876 * WORK(KSCR1),NTOTK) 3877C 3878 KOFF23 = IMATIJ(ISYMK,ISYMJ) + 1 3879 KOFF24 = IMATIJ(ISYMJ,ISYMI) + 1 3880C 3881 CALL DGEMM('N','N',NRHF(ISYMK),NRHF(ISYMI),NRHF(ISYMJ),-ONE, 3882 * DIJT(KOFF23),NTOTK,XINTIJ(KOFF24),NTOTJ,ONE, 3883 * WORK(KSCR1),NTOTK) 3884C 3885 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI),NRHF(ISYMK),ONE, 3886 * T1AM(KOFFT),NTOTR,WORK(KSCR1),NTOTK,ONE, 3887 * ETAKA(KOFFR),NTOTR) 3888C 3889 110 CONTINUE 3890C 3891C---------------------------------------------------- 3892C Calculate terms originating from [H(t1),E(ca)]. 3893C---------------------------------------------------- 3894C 3895 DO 120 ISYMA = 1,NSYM 3896C 3897 ISYMI = ISYMA 3898 ISYMC = ISYMI 3899 ISYMB = MULD2H(ISYMA,ISYM) 3900 ISYMJ = MULD2H(ISYMA,ISYM) 3901C 3902 KOFFR = IT1AM(ISYMA,ISYMI) + 1 3903 KOFFT = IT1AM(ISYMC,ISYMI) + 1 3904C 3905 NTOTR = MAX(NVIR(ISYMA),1) 3906 NTOTC = MAX(NVIR(ISYMC),1) 3907 NTOTB = MAX(NVIR(ISYMB),1) 3908C 3909C---------------------------------- 3910C Work space allocation two. 3911C---------------------------------- 3912C 3913 KSCR2 = 1 3914 KEND2 = KSCR2 + NVIR(ISYMA)*NVIR(ISYMC) 3915 LWRK2 = LWORK - KEND2 3916C 3917 IF (LWRK2 .LT. 0) THEN 3918 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND2 3919 CALL QUIT('Insufficient memory for work '// 3920 & 'allocation in CCDENZK0') 3921 ENDIF 3922C 3923 CALL DZERO(WORK(KSCR2),NVIR(ISYMA)*NVIR(ISYMC)) 3924C 3925 KOFF25 = IT1AM(ISYMA,ISYMJ) + 1 3926 KOFF26 = IT1AM(ISYMC,ISYMJ) + 1 3927C 3928 CALL DGEMM('N','T',NVIR(ISYMA),NVIR(ISYMC),NRHF(ISYMJ),ONE, 3929 * DIA(KOFF25),NTOTR,XINTIA(KOFF26),NTOTC,ONE, 3930 * WORK(KSCR2),NTOTR) 3931C 3932 KOFF27 = IT1AM(ISYMA,ISYMJ) + 1 3933 KOFF28 = IT1AM(ISYMC,ISYMJ) + 1 3934C 3935 CALL DGEMM('N','T',NVIR(ISYMA),NVIR(ISYMC),NRHF(ISYMJ),-ONE, 3936 * XINTAI(KOFF27),NTOTR,DAI(KOFF28),NTOTC,ONE, 3937 * WORK(KSCR2),NTOTR) 3938C 3939 KOFF29 = IMATAB(ISYMA,ISYMB) + 1 3940 KOFF30 = IMATAB(ISYMB,ISYMC) + 1 3941C 3942 CALL DGEMM('N','N',NVIR(ISYMA),NVIR(ISYMC),NVIR(ISYMB),ONE, 3943 * DABT(KOFF29),NTOTR,XINTAB(KOFF30),NTOTB,ONE, 3944 * WORK(KSCR2),NTOTR) 3945C 3946 KOFF31 = IMATAB(ISYMA,ISYMB) + 1 3947 KOFF32 = IMATAB(ISYMB,ISYMC) + 1 3948C 3949 CALL DGEMM('N','N',NVIR(ISYMA),NVIR(ISYMC),NVIR(ISYMB),-ONE, 3950 * XINTAB(KOFF31),NTOTR,DABT(KOFF32),NTOTB,ONE, 3951 * WORK(KSCR2),NTOTR) 3952C 3953 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI),NVIR(ISYMC),ONE, 3954 * WORK(KSCR2),NTOTR,T1AM(KOFFT),NTOTC,ONE, 3955 * ETAKA(KOFFR),NTOTR) 3956C 3957 120 CONTINUE 3958C 3959C---------------------------------------------------- 3960C Calculate terms originating from [H(t1),E(ck)]. 3961C---------------------------------------------------- 3962C 3963 DO 130 ISYMA = 1,NSYM 3964C 3965 ISYMI = ISYMA 3966 ISYMK = ISYMA 3967 ISYMC = ISYMI 3968 ISYMB = MULD2H(ISYMC,ISYM) 3969 ISYMJ = MULD2H(ISYMC,ISYM) 3970C 3971 KOFFR = IT1AM(ISYMA,ISYMI) + 1 3972 KOFTO = IT1AM(ISYMA,ISYMK) + 1 3973 KOFTI = IT1AM(ISYMC,ISYMI) + 1 3974C 3975 NTOTR = MAX(NVIR(ISYMA),1) 3976 NTOTC = MAX(NVIR(ISYMC),1) 3977 NTOTB = MAX(NVIR(ISYMB),1) 3978 NTOTK = MAX(NRHF(ISYMK),1) 3979 NTOTJ = MAX(NRHF(ISYMJ),1) 3980C 3981C------------------------------------ 3982C Work space allocation three. 3983C------------------------------------ 3984C 3985 KSCR3 = 1 3986 KSCR4 = KSCR3 + NVIR(ISYMC)*NRHF(ISYMK) 3987 KEND3 = KSCR4 + NRHF(ISYMK)*NRHF(ISYMI) 3988 LWRK3 = LWORK - KEND3 3989C 3990 IF (LWRK3 .LT. 0) THEN 3991 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND3 3992 CALL QUIT('Insufficient memory for work allocation '// 3993 & 'in CCDENZK0') 3994 ENDIF 3995C 3996 CALL DZERO(WORK(KSCR3),NVIR(ISYMC)*NRHF(ISYMK)) 3997 CALL DZERO(WORK(KSCR4),NRHF(ISYMK)*NRHF(ISYMI)) 3998C 3999 KOFF33 = IMATAB(ISYMC,ISYMB) + 1 4000 KOFF34 = IT1AM(ISYMB,ISYMK) + 1 4001C 4002 CALL DGEMM('N','N',NVIR(ISYMC),NRHF(ISYMK),NVIR(ISYMB),ONE, 4003 * XABT(KOFF33),NTOTC,DAI(KOFF34),NTOTB,ONE, 4004 * WORK(KSCR3),NTOTC) 4005C 4006 KOFF35 = IMATAB(ISYMC,ISYMB) + 1 4007 KOFF36 = IT1AM(ISYMB,ISYMK) + 1 4008C 4009 CALL DGEMM('N','N',NVIR(ISYMC),NRHF(ISYMK),NVIR(ISYMB),-ONE, 4010 * DAB(KOFF35),NTOTC,XINTIA(KOFF36),NTOTB,ONE, 4011 * WORK(KSCR3),NTOTC) 4012C 4013 KOFF37 = IT1AM(ISYMC,ISYMJ) + 1 4014 KOFF38 = IMATIJ(ISYMJ,ISYMK) + 1 4015C 4016 CALL DGEMM('N','N',NVIR(ISYMC),NRHF(ISYMK),NRHF(ISYMJ),ONE, 4017 * XINTIA(KOFF37),NTOTC,DIJ(KOFF38),NTOTJ,ONE, 4018 * WORK(KSCR3),NTOTC) 4019C 4020 KOFF39 = IT1AM(ISYMC,ISYMJ) + 1 4021 KOFF40 = IMATIJ(ISYMJ,ISYMK) + 1 4022C 4023 CALL DGEMM('N','N',NVIR(ISYMC),NRHF(ISYMK),NRHF(ISYMJ),-ONE, 4024 * DAI(KOFF39),NTOTC,XIJT(KOFF40),NTOTJ,ONE, 4025 * WORK(KSCR3),NTOTC) 4026C 4027 CALL DGEMM('T','N',NRHF(ISYMK),NRHF(ISYMI),NVIR(ISYMC),ONE, 4028 * WORK(KSCR3),NTOTC,T1AM(KOFTI),NTOTC,ZERO, 4029 * WORK(KSCR4),NTOTK) 4030C 4031 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI),NRHF(ISYMK),ONE, 4032 * T1AM(KOFTO),NTOTR,WORK(KSCR4),NTOTK,ONE, 4033 * ETAKA(KOFFR),NTOTR) 4034C 4035 130 CONTINUE 4036C 4037 CALL QEXIT('CCDENZK0') 4038C 4039 RETURN 4040 END 4041C /* Deck cc_d2eff */ 4042 SUBROUTINE CC_D2EFF(AODEN,G,ISYMG,IDEL,ISYMD,ZKABAO,DHFAO,ICON) 4043C 4044C Written by Asger Halkier 12/5 - 1998 4045C 4046C Version: 1.0 4047C 4048C Purpose: To add the extra terms that add to the 2-elctron 4049C coupled cluster lamda density to form the effective 4050C 2-electron density. 4051C 4052#include "implicit.h" 4053 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 4054 DIMENSION AODEN(*), ZKABAO(*), DHFAO(*) 4055#include "priunit.h" 4056#include "ccorb.h" 4057#include "ccsdsym.h" 4058C 4059 CALL QENTER('CC_D2EFF') 4060C 4061 FACI = ONE 4062 IF (ICON .EQ. 2) FACI = ONE/TWO 4063C 4064C----------------------- 4065C Add coulomb terms. 4066C----------------------- 4067C 4068 ISALBE = MULD2H(ISYMG,ISYMD) 4069 D = IDEL - IBAS(ISYMD) 4070C 4071 IF (ISYMG .EQ. ISYMD) THEN 4072 KOFFGD = IAODIS(ISYMG,ISYMD) + NBAS(ISYMG)*(D - 1) + G 4073 FAC1 = TWO*DHFAO(KOFFGD)*FACI 4074 CALL DAXPY(N2BST(ISALBE),FAC1,ZKABAO,1,AODEN,1) 4075 ENDIF 4076C 4077C------------------------ 4078C Add exchange terms. 4079C------------------------ 4080C 4081 ISYMB = ISYMG 4082 ISYMA = ISYMD 4083C 4084 DO 100 B = 1,NBAS(ISYMB) 4085C 4086 KOFFGB = IAODIS(ISYMG,ISYMB) + NBAS(ISYMG)*(B - 1) + G 4087 KOFFAD = IAODIS(ISYMA,ISYMD) + NBAS(ISYMA)*(D - 1) + 1 4088 KOFFAB = IAODIS(ISYMA,ISYMB) + NBAS(ISYMA)*(B - 1) + 1 4089C 4090 FAC2 = -DHFAO(KOFFGB)*FACI 4091 CALL DAXPY(NBAS(ISYMA),FAC2,ZKABAO(KOFFAD),1, 4092 * AODEN(KOFFAB),1) 4093C 4094 100 CONTINUE 4095C 4096 CALL QEXIT('CC_D2EFF') 4097C 4098 RETURN 4099 END 4100C /* Deck ccdiazk0 */ 4101 SUBROUTINE CCDIAZK0(ZKDIA,XINTIJ,XINTAI,XINTIA,XINTAB,DIJ,DAI, 4102 * DIA,DAB,XIJT,XABT,DIJT,DABT,T1AM, 4103 * WORK,LWORK,ISYM) 4104C 4105C Written by Asger Halkier 17/6 - 1998 4106C 4107C Purpose: To set up the right hand side for the diagonal blocks 4108C of kappa-bar-0 (ZKDIA) from MO-integrals (XI*), Coupled 4109C Cluster densities (D*) and t1-amplitudes (T1AM). 4110C 4111C Note that due to the symmetry in the formulas, this 4112C routine is able to handle both the one- and the two- 4113C electron contributions. 4114C 4115#include "implicit.h" 4116 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 4117 DIMENSION ZKDIA(*), XINTIJ(*), XINTAI(*), XINTIA(*), XINTAB(*) 4118 DIMENSION DIJ(*), DAI(*), DIA(*), DAB(*), XIJT(*), XABT(*) 4119 DIMENSION DIJT(*), DABT(*), T1AM(*), WORK(LWORK) 4120#include "priunit.h" 4121#include "ccorb.h" 4122#include "ccsdsym.h" 4123C 4124 CALL QENTER('CCDIAZK0') 4125C 4126C=================================== 4127C First we do the virtual block. 4128C=================================== 4129C 4130 DO 100 ISYMA = 1,NSYM 4131C 4132 ISYMB = ISYMA 4133 ISYMC = MULD2H(ISYMA,ISYM) 4134 ISYMI = MULD2H(ISYMA,ISYM) 4135C 4136 KOFFR = NMATIJ(1) + IMATAB(ISYMA,ISYMB) + 1 4137C 4138 NTOTA = MAX(NVIR(ISYMA),1) 4139 NTOTB = MAX(NVIR(ISYMB),1) 4140 NTOTC = MAX(NVIR(ISYMC),1) 4141C 4142C----------------------------------------------------------------- 4143C Direct contributions with intermediate loop over virtual. 4144C----------------------------------------------------------------- 4145C 4146 KOFF1 = IMATAB(ISYMA,ISYMC) + 1 4147 KOFF2 = IMATAB(ISYMC,ISYMB) + 1 4148C 4149 CALL DGEMM('N','N',NVIR(ISYMA),NVIR(ISYMB),NVIR(ISYMC),ONE, 4150 * XABT(KOFF1),NTOTA,DAB(KOFF2),NTOTC,ONE, 4151 * ZKDIA(KOFFR),NTOTA) 4152C 4153 KOFF3 = IMATAB(ISYMA,ISYMC) + 1 4154 KOFF4 = IMATAB(ISYMC,ISYMB) + 1 4155C 4156 CALL DGEMM('N','N',NVIR(ISYMA),NVIR(ISYMB),NVIR(ISYMC),-ONE, 4157 * DAB(KOFF3),NTOTA,XABT(KOFF4),NTOTC,ONE, 4158 * ZKDIA(KOFFR),NTOTA) 4159C 4160 KOFF5 = IMATAB(ISYMA,ISYMC) + 1 4161 KOFF6 = IMATAB(ISYMC,ISYMB) + 1 4162C 4163 CALL DGEMM('N','N',NVIR(ISYMA),NVIR(ISYMB),NVIR(ISYMC),-ONE, 4164 * DABT(KOFF5),NTOTA,XINTAB(KOFF6),NTOTC,ONE, 4165 * ZKDIA(KOFFR),NTOTA) 4166C 4167 KOFF7 = IMATAB(ISYMA,ISYMC) + 1 4168 KOFF8 = IMATAB(ISYMC,ISYMB) + 1 4169C 4170 CALL DGEMM('N','N',NVIR(ISYMA),NVIR(ISYMB),NVIR(ISYMC),ONE, 4171 * XINTAB(KOFF7),NTOTA,DABT(KOFF8),NTOTC,ONE, 4172 * ZKDIA(KOFFR),NTOTA) 4173C 4174C------------------------------------------------------------------ 4175C Direct contributions with intermediate loop over occupied. 4176C------------------------------------------------------------------ 4177C 4178 KOFF9 = IT1AM(ISYMA,ISYMI) + 1 4179 KOFF10 = IT1AM(ISYMB,ISYMI) + 1 4180C 4181 CALL DGEMM('N','T',NVIR(ISYMA),NVIR(ISYMB),NRHF(ISYMI),ONE, 4182 * XINTIA(KOFF9),NTOTA,DIA(KOFF10),NTOTB,ONE, 4183 * ZKDIA(KOFFR),NTOTA) 4184C 4185 KOFF11 = IT1AM(ISYMA,ISYMI) + 1 4186 KOFF12 = IT1AM(ISYMB,ISYMI) + 1 4187C 4188 CALL DGEMM('N','T',NVIR(ISYMA),NVIR(ISYMB),NRHF(ISYMI),-ONE, 4189 * DAI(KOFF11),NTOTA,XINTAI(KOFF12),NTOTB,ONE, 4190 * ZKDIA(KOFFR),NTOTA) 4191C 4192 KOFF13 = IT1AM(ISYMA,ISYMI) + 1 4193 KOFF14 = IT1AM(ISYMB,ISYMI) + 1 4194C 4195 CALL DGEMM('N','T',NVIR(ISYMA),NVIR(ISYMB),NRHF(ISYMI),-ONE, 4196 * DIA(KOFF13),NTOTA,XINTIA(KOFF14),NTOTB,ONE, 4197 * ZKDIA(KOFFR),NTOTA) 4198C 4199 KOFF15 = IT1AM(ISYMA,ISYMI) + 1 4200 KOFF16 = IT1AM(ISYMB,ISYMI) + 1 4201C 4202 CALL DGEMM('N','T',NVIR(ISYMA),NVIR(ISYMB),NRHF(ISYMI),ONE, 4203 * XINTAI(KOFF15),NTOTA,DAI(KOFF16),NTOTB,ONE, 4204 * ZKDIA(KOFFR),NTOTA) 4205C 4206 100 CONTINUE 4207C 4208 DO 110 ISYMA = 1,NSYM 4209C 4210 ISYMB = ISYMA 4211 ISYMK = ISYMA 4212 ISYMC = MULD2H(ISYMA,ISYM) 4213 ISYMI = MULD2H(ISYMA,ISYM) 4214C 4215 KOFFR = NMATIJ(1) + IMATAB(ISYMA,ISYMB) + 1 4216 KOFFT = IT1AM(ISYMA,ISYMK) + 1 4217C 4218 NTOTA = MAX(NVIR(ISYMA),1) 4219 NTOTB = MAX(NVIR(ISYMB),1) 4220 NTOTC = MAX(NVIR(ISYMC),1) 4221 NTOTI = MAX(NRHF(ISYMI),1) 4222C 4223C---------------------------------- 4224C Work space allocation one. 4225C---------------------------------- 4226C 4227 KSCR1 = 1 4228 KEND1 = KSCR1 + NVIR(ISYMA)*NRHF(ISYMK) 4229 LWRK1 = LWORK - KEND1 4230C 4231 IF (LWRK1 .LT. 0) THEN 4232 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 4233 CALL QUIT('Insufficient memory for work allocation '// 4234 & 'in CCDIAZK0') 4235 ENDIF 4236C 4237 LEN = NVIR(ISYMA)*NRHF(ISYMK) 4238C 4239 CALL DZERO(WORK(KSCR1),LEN) 4240C 4241C------------------------------------------------------------------- 4242C Indirect contributions with intermediate loop over virtual. 4243C------------------------------------------------------------------- 4244C 4245 KOFF1 = IMATAB(ISYMA,ISYMC) + 1 4246 KOFF2 = IT1AM(ISYMC,ISYMK) + 1 4247C 4248 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMK),NVIR(ISYMC),ONE, 4249 * XABT(KOFF1),NTOTA,DAI(KOFF2),NTOTC,ZERO, 4250 * WORK(KSCR1),NTOTA) 4251C 4252 CALL DGEMM('N','T',NVIR(ISYMA),NVIR(ISYMB),NRHF(ISYMK),ONE, 4253 * WORK(KSCR1),NTOTA,T1AM(KOFFT),NTOTB,ONE, 4254 * ZKDIA(KOFFR),NTOTA) 4255C 4256 CALL DZERO(WORK(KSCR1),LEN) 4257C 4258 KOFF3 = IMATAB(ISYMB,ISYMC) + 1 4259 KOFF4 = IT1AM(ISYMC,ISYMK) + 1 4260C 4261 CALL DGEMM('N','N',NVIR(ISYMB),NRHF(ISYMK),NVIR(ISYMC),ONE, 4262 * XABT(KOFF3),NTOTB,DAI(KOFF4),NTOTC,ZERO, 4263 * WORK(KSCR1),NTOTB) 4264C 4265 CALL DGEMM('N','T',NVIR(ISYMA),NVIR(ISYMB),NRHF(ISYMK),-ONE, 4266 * T1AM(KOFFT),NTOTA,WORK(KSCR1),NTOTB,ONE, 4267 * ZKDIA(KOFFR),NTOTA) 4268C 4269 CALL DZERO(WORK(KSCR1),LEN) 4270C 4271 KOFF5 = IMATAB(ISYMA,ISYMC) + 1 4272 KOFF6 = IT1AM(ISYMC,ISYMK) + 1 4273C 4274 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMK),NVIR(ISYMC),ONE, 4275 * DAB(KOFF5),NTOTA,XINTIA(KOFF6),NTOTC,ZERO, 4276 * WORK(KSCR1),NTOTA) 4277C 4278 CALL DGEMM('N','T',NVIR(ISYMA),NVIR(ISYMB),NRHF(ISYMK),-ONE, 4279 * WORK(KSCR1),NTOTA,T1AM(KOFFT),NTOTB,ONE, 4280 * ZKDIA(KOFFR),NTOTA) 4281C 4282 CALL DZERO(WORK(KSCR1),LEN) 4283C 4284 KOFF7 = IMATAB(ISYMB,ISYMC) + 1 4285 KOFF8 = IT1AM(ISYMC,ISYMK) + 1 4286C 4287 CALL DGEMM('N','N',NVIR(ISYMB),NRHF(ISYMK),NVIR(ISYMC),ONE, 4288 * DAB(KOFF7),NTOTB,XINTIA(KOFF8),NTOTC,ZERO, 4289 * WORK(KSCR1),NTOTB) 4290C 4291 CALL DGEMM('N','T',NVIR(ISYMA),NVIR(ISYMB),NRHF(ISYMK),ONE, 4292 * T1AM(KOFFT),NTOTA,WORK(KSCR1),NTOTB,ONE, 4293 * ZKDIA(KOFFR),NTOTA) 4294C 4295 CALL DZERO(WORK(KSCR1),LEN) 4296C 4297C-------------------------------------------------------------------- 4298C Indirect contributions with intermediate loop over occupied. 4299C-------------------------------------------------------------------- 4300C 4301 KOFF9 = IT1AM(ISYMA,ISYMI) + 1 4302 KOFF10 = IMATIJ(ISYMI,ISYMK) + 1 4303C 4304 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMK),NRHF(ISYMI),ONE, 4305 * XINTIA(KOFF9),NTOTA,DIJ(KOFF10),NTOTI,ZERO, 4306 * WORK(KSCR1),NTOTA) 4307C 4308 CALL DGEMM('N','T',NVIR(ISYMA),NVIR(ISYMB),NRHF(ISYMK),ONE, 4309 * WORK(KSCR1),NTOTA,T1AM(KOFFT),NTOTB,ONE, 4310 * ZKDIA(KOFFR),NTOTA) 4311C 4312 CALL DZERO(WORK(KSCR1),LEN) 4313C 4314 KOFF11 = IT1AM(ISYMB,ISYMI) + 1 4315 KOFF12 = IMATIJ(ISYMI,ISYMK) + 1 4316C 4317 CALL DGEMM('N','N',NVIR(ISYMB),NRHF(ISYMK),NRHF(ISYMI),ONE, 4318 * XINTIA(KOFF11),NTOTB,DIJ(KOFF12),NTOTI,ZERO, 4319 * WORK(KSCR1),NTOTB) 4320C 4321 CALL DGEMM('N','T',NVIR(ISYMA),NVIR(ISYMB),NRHF(ISYMK),-ONE, 4322 * T1AM(KOFFT),NTOTA,WORK(KSCR1),NTOTB,ONE, 4323 * ZKDIA(KOFFR),NTOTA) 4324C 4325 CALL DZERO(WORK(KSCR1),LEN) 4326C 4327 KOFF13 = IT1AM(ISYMA,ISYMI) + 1 4328 KOFF14 = IMATIJ(ISYMI,ISYMK) + 1 4329C 4330 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMK),NRHF(ISYMI),ONE, 4331 * DAI(KOFF13),NTOTA,XIJT(KOFF14),NTOTI,ZERO, 4332 * WORK(KSCR1),NTOTA) 4333C 4334 CALL DGEMM('N','T',NVIR(ISYMA),NVIR(ISYMB),NRHF(ISYMK),-ONE, 4335 * WORK(KSCR1),NTOTA,T1AM(KOFFT),NTOTB,ONE, 4336 * ZKDIA(KOFFR),NTOTA) 4337C 4338 CALL DZERO(WORK(KSCR1),LEN) 4339C 4340 KOFF15 = IT1AM(ISYMB,ISYMI) + 1 4341 KOFF16 = IMATIJ(ISYMI,ISYMK) + 1 4342C 4343 CALL DGEMM('N','N',NVIR(ISYMB),NRHF(ISYMK),NRHF(ISYMI),ONE, 4344 * DAI(KOFF15),NTOTB,XIJT(KOFF16),NTOTI,ZERO, 4345 * WORK(KSCR1),NTOTB) 4346C 4347 CALL DGEMM('N','T',NVIR(ISYMA),NVIR(ISYMB),NRHF(ISYMK),ONE, 4348 * T1AM(KOFFT),NTOTA,WORK(KSCR1),NTOTB,ONE, 4349 * ZKDIA(KOFFR),NTOTA) 4350C 4351 110 CONTINUE 4352C 4353C=================================== 4354C Then we do the occupied block. 4355C=================================== 4356C 4357 DO 120 ISYMI = 1,NSYM 4358C 4359 ISYMJ = ISYMI 4360 ISYMC = MULD2H(ISYMI,ISYM) 4361 ISYMK = MULD2H(ISYMI,ISYM) 4362C 4363 KOFFR = IMATIJ(ISYMI,ISYMJ) + 1 4364C 4365 NTOTI = MAX(NRHF(ISYMI),1) 4366 NTOTC = MAX(NVIR(ISYMC),1) 4367 NTOTK = MAX(NRHF(ISYMK),1) 4368C 4369C----------------------------------------------------------------- 4370C Direct contributions with intermediate loop over virtual. 4371C----------------------------------------------------------------- 4372C 4373 KOFF17 = IT1AM(ISYMC,ISYMI) + 1 4374 KOFF18 = IT1AM(ISYMC,ISYMJ) + 1 4375C 4376 CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),NVIR(ISYMC),ONE, 4377 * XINTAI(KOFF17),NTOTC,DAI(KOFF18),NTOTC,ONE, 4378 * ZKDIA(KOFFR),NTOTI) 4379C 4380 KOFF19 = IT1AM(ISYMC,ISYMI) + 1 4381 KOFF20 = IT1AM(ISYMC,ISYMJ) + 1 4382C 4383 CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),NVIR(ISYMC),-ONE, 4384 * DAI(KOFF19),NTOTC,XINTAI(KOFF20),NTOTC,ONE, 4385 * ZKDIA(KOFFR),NTOTI) 4386C 4387 KOFF21 = IT1AM(ISYMC,ISYMI) + 1 4388 KOFF22 = IT1AM(ISYMC,ISYMJ) + 1 4389C 4390 CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),NVIR(ISYMC),-ONE, 4391 * DIA(KOFF21),NTOTC,XINTIA(KOFF22),NTOTC,ONE, 4392 * ZKDIA(KOFFR),NTOTI) 4393C 4394 KOFF23 = IT1AM(ISYMC,ISYMI) + 1 4395 KOFF24 = IT1AM(ISYMC,ISYMJ) + 1 4396C 4397 CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),NVIR(ISYMC),ONE, 4398 * XINTIA(KOFF23),NTOTC,DIA(KOFF24),NTOTC,ONE, 4399 * ZKDIA(KOFFR),NTOTI) 4400C 4401C------------------------------------------------------------------ 4402C Direct contributions with intermediate loop over occupied. 4403C------------------------------------------------------------------ 4404C 4405 KOFF25 = IMATIJ(ISYMI,ISYMK) + 1 4406 KOFF26 = IMATIJ(ISYMK,ISYMJ) + 1 4407C 4408 CALL DGEMM('N','N',NRHF(ISYMI),NRHF(ISYMJ),NRHF(ISYMK),ONE, 4409 * XIJT(KOFF25),NTOTI,DIJ(KOFF26),NTOTK,ONE, 4410 * ZKDIA(KOFFR),NTOTI) 4411C 4412 KOFF27 = IMATIJ(ISYMI,ISYMK) + 1 4413 KOFF28 = IMATIJ(ISYMK,ISYMJ) + 1 4414C 4415 CALL DGEMM('N','N',NRHF(ISYMI),NRHF(ISYMJ),NRHF(ISYMK),-ONE, 4416 * DIJT(KOFF27),NTOTI,XINTIJ(KOFF28),NTOTK,ONE, 4417 * ZKDIA(KOFFR),NTOTI) 4418C 4419 KOFF29 = IMATIJ(ISYMI,ISYMK) + 1 4420 KOFF30 = IMATIJ(ISYMK,ISYMJ) + 1 4421C 4422 CALL DGEMM('N','N',NRHF(ISYMI),NRHF(ISYMJ),NRHF(ISYMK),-ONE, 4423 * DIJ(KOFF29),NTOTI,XIJT(KOFF30),NTOTK,ONE, 4424 * ZKDIA(KOFFR),NTOTI) 4425C 4426 KOFF31 = IMATIJ(ISYMI,ISYMK) + 1 4427 KOFF32 = IMATIJ(ISYMK,ISYMJ) + 1 4428C 4429 CALL DGEMM('N','N',NRHF(ISYMI),NRHF(ISYMJ),NRHF(ISYMK),ONE, 4430 * XINTIJ(KOFF31),NTOTI,DIJT(KOFF32),NTOTK,ONE, 4431 * ZKDIA(KOFFR),NTOTI) 4432C 4433 120 CONTINUE 4434C 4435 DO 130 ISYMI = 1,NSYM 4436C 4437 ISYMJ = ISYMI 4438 ISYMC = ISYMI 4439 ISYMD = MULD2H(ISYMI,ISYM) 4440 ISYMK = MULD2H(ISYMI,ISYM) 4441C 4442 KOFFR = IMATIJ(ISYMI,ISYMJ) + 1 4443 KOFFT = IT1AM(ISYMC,ISYMI) + 1 4444C 4445 NTOTI = MAX(NRHF(ISYMI),1) 4446 NTOTK = MAX(NRHF(ISYMK),1) 4447 NTOTC = MAX(NVIR(ISYMC),1) 4448 NTOTD = MAX(NVIR(ISYMD),1) 4449C 4450C---------------------------------- 4451C Work space allocation two. 4452C---------------------------------- 4453C 4454 KSCR2 = 1 4455 KEND2 = KSCR2 + NVIR(ISYMC)*NRHF(ISYMI) 4456 LWRK2 = LWORK - KEND2 4457C 4458 IF (LWRK2 .LT. 0) THEN 4459 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND2 4460 CALL QUIT('Insufficient memory for work allocation '// 4461 & 'in CCDIAZK0') 4462 ENDIF 4463C 4464 LEN = NVIR(ISYMC)*NRHF(ISYMI) 4465C 4466 CALL DZERO(WORK(KSCR2),LEN) 4467C 4468C------------------------------------------------------------------- 4469C Indirect contributions with intermediate loop over virtual. 4470C------------------------------------------------------------------- 4471C 4472 KOFF33 = IMATAB(ISYMC,ISYMD) + 1 4473 KOFF34 = IT1AM(ISYMD,ISYMI) + 1 4474C 4475 CALL DGEMM('N','N',NVIR(ISYMC),NRHF(ISYMI),NVIR(ISYMD),ONE, 4476 * XABT(KOFF33),NTOTC,DAI(KOFF34),NTOTD,ZERO, 4477 * WORK(KSCR2),NTOTC) 4478C 4479 CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),NVIR(ISYMC),ONE, 4480 * WORK(KSCR2),NTOTC,T1AM(KOFFT),NTOTC,ONE, 4481 * ZKDIA(KOFFR),NTOTI) 4482C 4483 CALL DZERO(WORK(KSCR2),LEN) 4484C 4485 KOFF35 = IMATAB(ISYMC,ISYMD) + 1 4486 KOFF36 = IT1AM(ISYMD,ISYMJ) + 1 4487C 4488 CALL DGEMM('N','N',NVIR(ISYMC),NRHF(ISYMJ),NVIR(ISYMD),ONE, 4489 * XABT(KOFF35),NTOTC,DAI(KOFF36),NTOTD,ZERO, 4490 * WORK(KSCR2),NTOTC) 4491C 4492 CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),NVIR(ISYMC),-ONE, 4493 * T1AM(KOFFT),NTOTC,WORK(KSCR2),NTOTC,ONE, 4494 * ZKDIA(KOFFR),NTOTI) 4495C 4496 CALL DZERO(WORK(KSCR2),LEN) 4497C 4498 KOFF37 = IMATAB(ISYMC,ISYMD) + 1 4499 KOFF38 = IT1AM(ISYMD,ISYMI) + 1 4500C 4501 CALL DGEMM('N','N',NVIR(ISYMC),NRHF(ISYMI),NVIR(ISYMD),ONE, 4502 * DAB(KOFF37),NTOTC,XINTIA(KOFF38),NTOTD,ZERO, 4503 * WORK(KSCR2),NTOTC) 4504C 4505 CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),NVIR(ISYMC),-ONE, 4506 * WORK(KSCR2),NTOTC,T1AM(KOFFT),NTOTC,ONE, 4507 * ZKDIA(KOFFR),NTOTI) 4508C 4509 CALL DZERO(WORK(KSCR2),LEN) 4510C 4511 KOFF39 = IMATAB(ISYMC,ISYMD) + 1 4512 KOFF40 = IT1AM(ISYMD,ISYMJ) + 1 4513C 4514 CALL DGEMM('N','N',NVIR(ISYMC),NRHF(ISYMJ),NVIR(ISYMD),ONE, 4515 * DAB(KOFF39),NTOTC,XINTIA(KOFF40),NTOTD,ZERO, 4516 * WORK(KSCR2),NTOTC) 4517C 4518 CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),NVIR(ISYMC),ONE, 4519 * T1AM(KOFFT),NTOTC,WORK(KSCR2),NTOTC,ONE, 4520 * ZKDIA(KOFFR),NTOTI) 4521C 4522C------------------------------------------------------------------- 4523C Indirect contributions with intermediate loop over virtual. 4524C------------------------------------------------------------------- 4525C 4526 CALL DZERO(WORK(KSCR2),LEN) 4527C 4528 KOFF41 = IT1AM(ISYMC,ISYMK) + 1 4529 KOFF42 = IMATIJ(ISYMK,ISYMI) + 1 4530C 4531 CALL DGEMM('N','N',NVIR(ISYMC),NRHF(ISYMI),NRHF(ISYMK),ONE, 4532 * XINTIA(KOFF41),NTOTC,DIJ(KOFF42),NTOTK,ZERO, 4533 * WORK(KSCR2),NTOTC) 4534C 4535 CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),NVIR(ISYMC),ONE, 4536 * WORK(KSCR2),NTOTC,T1AM(KOFFT),NTOTC,ONE, 4537 * ZKDIA(KOFFR),NTOTI) 4538C 4539 CALL DZERO(WORK(KSCR2),LEN) 4540C 4541 KOFF43 = IT1AM(ISYMC,ISYMK) + 1 4542 KOFF44 = IMATIJ(ISYMK,ISYMJ) + 1 4543C 4544 CALL DGEMM('N','N',NVIR(ISYMC),NRHF(ISYMJ),NRHF(ISYMK),ONE, 4545 * XINTIA(KOFF43),NTOTC,DIJ(KOFF44),NTOTK,ZERO, 4546 * WORK(KSCR2),NTOTC) 4547C 4548 CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),NVIR(ISYMC),-ONE, 4549 * T1AM(KOFFT),NTOTC,WORK(KSCR2),NTOTC,ONE, 4550 * ZKDIA(KOFFR),NTOTI) 4551C 4552 CALL DZERO(WORK(KSCR2),LEN) 4553C 4554 KOFF45 = IT1AM(ISYMC,ISYMK) + 1 4555 KOFF46 = IMATIJ(ISYMK,ISYMI) + 1 4556C 4557 CALL DGEMM('N','N',NVIR(ISYMC),NRHF(ISYMI),NRHF(ISYMK),ONE, 4558 * DAI(KOFF45),NTOTC,XIJT(KOFF46),NTOTK,ZERO, 4559 * WORK(KSCR2),NTOTC) 4560C 4561 CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),NVIR(ISYMC),-ONE, 4562 * WORK(KSCR2),NTOTC,T1AM(KOFFT),NTOTC,ONE, 4563 * ZKDIA(KOFFR),NTOTI) 4564C 4565 CALL DZERO(WORK(KSCR2),LEN) 4566C 4567 KOFF47 = IT1AM(ISYMC,ISYMK) + 1 4568 KOFF48 = IMATIJ(ISYMK,ISYMJ) + 1 4569C 4570 CALL DGEMM('N','N',NVIR(ISYMC),NRHF(ISYMJ),NRHF(ISYMK),ONE, 4571 * DAI(KOFF47),NTOTC,XIJT(KOFF48),NTOTK,ZERO, 4572 * WORK(KSCR2),NTOTC) 4573C 4574 CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),NVIR(ISYMC),ONE, 4575 * T1AM(KOFFT),NTOTC,WORK(KSCR2),NTOTC,ONE, 4576 * ZKDIA(KOFFR),NTOTI) 4577C 4578 130 CONTINUE 4579C 4580 CALL QEXIT('CCDIAZK0') 4581C 4582 RETURN 4583 END 4584C /* Deck ccsd_zkblo */ 4585 SUBROUTINE CCSD_ZKBLO(ZKDIA,WORK,LWORK) 4586C 4587C Written by Asger Halkier 4/8 - 1998 4588C 4589C Version: 1.0 4590C 4591C Purpose: To calculate the ab & ij parts of the CCSD kappa-bar-0, 4592C from right-hand-sides (ZKDIA on input) and canonical 4593C orbital energies. 4594C Control numerical instabilities. Sonia, 2002 4595C Additional numerical stability, Thomas Bondo Pedersen, Jan. 2013. 4596C - if RHS (eta) is zero, then kappa-bar-0 is set to zero. 4597C - if RHS is non-zero and denominator is zero, the equation 4598C system is singular and we have to quit. 4599C - in addition, redundant copying and zeroing eliminated. 4600C 4601#include "implicit.h" 4602#include "dummy.h" 4603 PARAMETER(HALF = 0.5D0) 4604 PARAMETER (EPSN = 1.0D-12, EPSD=1.0d-12) 4605 DIMENSION ZKDIA(*), WORK(LWORK) 4606#include "priunit.h" 4607#include "maxorb.h" 4608#include "ccorb.h" 4609#include "iratdef.h" 4610#include "inftap.h" 4611#include "ccsdsym.h" 4612#include "ccsdio.h" 4613#include "ccsdinp.h" 4614C 4615 CALL QENTER('CCSD_ZKBLO') 4616C 4617C--------------------------- 4618C Work space allocation. 4619C--------------------------- 4620C 4621 KFOCKD = 1 4622 KEND1 = KFOCKD + NORBTS 4623 LWRK1 = LWORK - KEND1 4624C 4625 IF (LWRK1 .LT. 0) THEN 4626 WRITE(LUPRI,*) 'Need:', KEND1, 'Available:', LWORK 4627 CALL QUIT('Insufficient memory for allocation in CCSD_ZKBLO') 4628 ENDIF 4629C 4630C------------------------------------- 4631C Read canonical orbital energies. 4632C------------------------------------- 4633C 4634 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY, 4635 & .FALSE.) 4636 REWIND (LUSIFC) 4637C 4638 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 4639 READ (LUSIFC) 4640 READ (LUSIFC) (WORK(KFOCKD + I - 1), I = 1,NORBTS) 4641C 4642 CALL GPCLOSE(LUSIFC,'KEEP') 4643C 4644C---------------------------------------------------------------- 4645C Change symmetry ordering of the canonical orbital energies. 4646C---------------------------------------------------------------- 4647C 4648 IF (FROIMP) 4649 * CALL CCSD_DELFRO(WORK(KFOCKD),WORK(KEND1),LWRK1) 4650C 4651 CALL FOCK_REORDER(WORK(KFOCKD),WORK(KEND1),LWRK1) 4652C 4653C--------------------------- 4654C Calculate the results: 4655C Occupied block: 4656C--------------------------- 4657C 4658 DO ISYMI = 1,NSYM 4659 ISYMJ = ISYMI 4660 DO J = 1,NRHF(ISYMJ) 4661 NJJ = IMATIJ(ISYMJ,ISYMJ) + NRHF(ISYMJ)*(J - 1) + J 4662 ZKDIA(NJJ) = 0.0D0 4663 KOFFJ = KFOCKD + IRHF(ISYMJ) + J - 1 4664 DO I = J+1,NRHF(ISYMI) 4665 KOFFI = KFOCKD + IRHF(ISYMI) + I - 1 4666 DENOM = WORK(KOFFJ) - WORK(KOFFI) 4667 NIJ = IMATIJ(ISYMI,ISYMJ) + NRHF(ISYMI)*(J - 1) + I 4668 ETAIJ = HALF*ZKDIA(NIJ) 4669 ZKDIA(NIJ) = CC_PROTECTED_DIVISION(ETAIJ,DENOM,EPSN,EPSD) 4670 NJI = IMATIJ(ISYMJ,ISYMI) + NRHF(ISYMJ)*(I - 1) + J 4671 ZKDIA(NJI) = ZKDIA(NIJ) 4672 END DO 4673 END DO 4674 END DO 4675C 4676C------------------- 4677C Virtual block: 4678C------------------- 4679C 4680 DO ISYMA = 1,NSYM 4681 ISYMB = ISYMA 4682 DO B = 1,NVIR(ISYMB) 4683 NBB = NMATIJ(1)+IMATAB(ISYMB,ISYMB)+NVIR(ISYMB)*(B-1)+B 4684 ZKDIA(NBB) = 0.0D0 4685 KOFFB = KFOCKD + IVIR(ISYMB) + B - 1 4686 DO A = B+1,NVIR(ISYMA) 4687 KOFFA = KFOCKD + IVIR(ISYMA) + A - 1 4688 DENOM = WORK(KOFFB) - WORK(KOFFA) 4689 NAB = NMATIJ(1)+IMATAB(ISYMA,ISYMB)+NVIR(ISYMA)*(B-1)+A 4690 ETAAB = HALF*ZKDIA(NAB) 4691 ZKDIA(NAB) = CC_PROTECTED_DIVISION(ETAAB,DENOM,EPSN,EPSD) 4692 NBA = NMATIJ(1)+IMATAB(ISYMB,ISYMA)+NVIR(ISYMB)*(A-1)+B 4693 ZKDIA(NBA) = ZKDIA(NAB) 4694 END DO 4695 END DO 4696 END DO 4697C 4698 CALL QEXIT('CCSD_ZKBLO') 4699C 4700 RETURN 4701 END 4702C /* Deck cc_d2gaf */ 4703 SUBROUTINE CC_D2GAF(D2GIJ,D2GAB,D2GAI,D2GIA,DIJ,DAB,DAI,DIA, 4704 * CMO,IDEL,ISYMD,G,ISYMG) 4705C 4706C Written by Asger Halkier 12/8 - 1998 4707C 4708C Purpose: To calculate the contributions to d(pq,gam;del) where 4709C gamma has been backtransformed from a frozen index. 4710C 4711#include "implicit.h" 4712 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 4713 DIMENSION D2GIJ(*), D2GAB(*), D2GAI(*), D2GIA(*) 4714 DIMENSION DIJ(*), DAB(*), DAI(*), DIA(*), CMO(*) 4715#include "priunit.h" 4716#include "ccorb.h" 4717#include "ccsdsym.h" 4718C 4719 CALL QENTER('CC_D2GAF') 4720C 4721 IF (ISYMG .EQ. ISYMD) THEN 4722C 4723 ISYML = ISYMD 4724C 4725 ND = ILRHSI(ISYMD) + IDEL - IBAS(ISYMD) 4726 NG = ILRHSI(ISYMG) + G 4727C 4728 FACT = TWO*DDOT(NRHFFR(ISYML),CMO(ND),NBAS(ISYMD),CMO(NG), 4729 * NBAS(ISYMG)) 4730C 4731 CALL DAXPY(NMATIJ(1),FACT,DIJ,1,D2GIJ,1) 4732 CALL DAXPY(NMATAB(1),FACT,DAB,1,D2GAB,1) 4733 CALL DAXPY(NT1AMX,FACT,DAI,1,D2GAI,1) 4734 CALL DAXPY(NT1AMX,FACT,DIA,1,D2GIA,1) 4735C 4736 ENDIF 4737C 4738 CALL QEXIT('CC_D2GAF') 4739C 4740 RETURN 4741 END 4742C /* Deck cc_fd2ao */ 4743 SUBROUTINE CC_FD2AO(D2AO,D2II,D2IJ,D2JI,D2AI,D2IA,CMO,XLAMDH, 4744 * XLAMDP,WORK,LWORK,ISYMPQ) 4745C 4746C Written by Asger Halkier 12/8 - 1998 4747C 4748C Purpose: To calculate the contributions to D2AO from d(pq,gam;del) 4749C where at least one of the two indices p & q is frozen. 4750C 4751#include "implicit.h" 4752 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 4753 DIMENSION D2AO(*), D2II(*), D2IJ(*), D2JI(*), D2AI(*) 4754 DIMENSION D2IA(*), CMO(*), XLAMDH(*), XLAMDP(*), WORK(LWORK) 4755#include "priunit.h" 4756#include "ccorb.h" 4757#include "ccsdsym.h" 4758#include "ccfro.h" 4759C 4760 CALL QENTER('CC_FD2AO') 4761C 4762C-------------------------------- 4763C Do the frozen-frozen block. 4764C-------------------------------- 4765C 4766 DO 100 ISYMAL = 1,NSYM 4767C 4768 ISYMBE = MULD2H(ISYMPQ,ISYMAL) 4769 ISYMI = ISYMAL 4770 ISYMJ = ISYMBE 4771C 4772 IF ((NRHFFR(ISYMI) .EQ. 0).OR.(NRHFFR(ISYMJ) .EQ. 0)) GOTO 100 4773C 4774 KSCR1 = 1 4775 KEND1 = KSCR1 + NBAS(ISYMAL)*NRHFFR(ISYMJ) 4776 LWRK1 = LWORK - KEND1 4777C 4778 IF (LWRK1 .LT. 0) THEN 4779 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 4780 CALL QUIT('Insufficient work space in CC_FD2AO') 4781 ENDIF 4782C 4783 CALL DZERO(WORK(KSCR1),KEND1) 4784C 4785 KOFF1 = ILRHSI(ISYMAL) + 1 4786 KOFF2 = IFROFR(ISYMI,ISYMJ) + 1 4787C 4788 NTOTAL = MAX(NBAS(ISYMAL),1) 4789 NTOTI = MAX(NRHFFR(ISYMI),1) 4790C 4791 CALL DGEMM('N','N',NBAS(ISYMAL),NRHFFR(ISYMJ),NRHFFR(ISYMI), 4792 * ONE,CMO(KOFF1),NTOTAL,D2II(KOFF2),NTOTI,ZERO, 4793 * WORK(KSCR1),NTOTAL) 4794C 4795 KOFF3 = ILRHSI(ISYMBE) + 1 4796 KOFF4 = IAODIS(ISYMAL,ISYMBE) + 1 4797C 4798 NTOTAL = MAX(NBAS(ISYMAL),1) 4799 NTOTBE = MAX(NBAS(ISYMBE),1) 4800C 4801 CALL DGEMM('N','T',NBAS(ISYMAL),NBAS(ISYMBE),NRHFFR(ISYMJ), 4802 * ONE,WORK(KSCR1),NTOTAL,CMO(KOFF3),NTOTBE,ONE, 4803 * D2AO(KOFF4),NTOTAL) 4804C 4805 100 CONTINUE 4806C 4807C------------------------------------ 4808C Do the correlated-frozen block. 4809C------------------------------------ 4810C 4811 DO 110 ISYMAL = 1,NSYM 4812C 4813 ISYMBE = MULD2H(ISYMPQ,ISYMAL) 4814 ISYMI = ISYMAL 4815 ISYMJ = ISYMBE 4816C 4817 IF (NRHFFR(ISYMJ) .EQ. 0) GOTO 110 4818C 4819 KSCR1 = 1 4820 KEND1 = KSCR1 + NBAS(ISYMAL)*NRHFFR(ISYMJ) 4821 LWRK1 = LWORK - KEND1 4822C 4823 IF (LWRK1 .LT. 0) THEN 4824 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 4825 CALL QUIT('Insufficient work space in CC_FD2AO') 4826 ENDIF 4827C 4828 CALL DZERO(WORK(KSCR1),KEND1) 4829C 4830 KOFF5 = IGLMRH(ISYMAL,ISYMI) + 1 4831 KOFF6 = ICOFRO(ISYMI,ISYMJ) + 1 4832C 4833 NTOTAL = MAX(NBAS(ISYMAL),1) 4834 NTOTI = MAX(NRHF(ISYMI),1) 4835C 4836 CALL DGEMM('N','N',NBAS(ISYMAL),NRHFFR(ISYMJ),NRHF(ISYMI), 4837 * ONE,XLAMDP(KOFF5),NTOTAL,D2IJ(KOFF6),NTOTI,ZERO, 4838 * WORK(KSCR1),NTOTAL) 4839C 4840 KOFF7 = ILRHSI(ISYMBE) + 1 4841 KOFF8 = IAODIS(ISYMAL,ISYMBE) + 1 4842C 4843 NTOTAL = MAX(NBAS(ISYMAL),1) 4844 NTOTBE = MAX(NBAS(ISYMBE),1) 4845C 4846 CALL DGEMM('N','T',NBAS(ISYMAL),NBAS(ISYMBE),NRHFFR(ISYMJ), 4847 * ONE,WORK(KSCR1),NTOTAL,CMO(KOFF7),NTOTBE,ONE, 4848 * D2AO(KOFF8),NTOTAL) 4849C 4850 110 CONTINUE 4851C 4852C------------------------------------ 4853C Do the frozen-correlated block. 4854C------------------------------------ 4855C 4856 DO 120 ISYMAL = 1,NSYM 4857C 4858 ISYMBE = MULD2H(ISYMPQ,ISYMAL) 4859 ISYMI = ISYMAL 4860 ISYMJ = ISYMBE 4861C 4862 IF (NRHFFR(ISYMI) .EQ. 0) GOTO 120 4863C 4864 KSCR1 = 1 4865 KEND1 = KSCR1 + NBAS(ISYMBE)*NRHFFR(ISYMI) 4866 LWRK1 = LWORK - KEND1 4867C 4868 IF (LWRK1 .LT. 0) THEN 4869 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 4870 CALL QUIT('Insufficient work space in CC_FD2AO') 4871 ENDIF 4872C 4873 CALL DZERO(WORK(KSCR1),KEND1) 4874C 4875 KOFF9 = IGLMRH(ISYMBE,ISYMJ) + 1 4876 KOFF10 = ICOFRO(ISYMJ,ISYMI) + 1 4877C 4878 NTOTBE = MAX(NBAS(ISYMBE),1) 4879 NTOTJ = MAX(NRHF(ISYMJ),1) 4880C 4881 CALL DGEMM('N','N',NBAS(ISYMBE),NRHFFR(ISYMI),NRHF(ISYMJ), 4882 * ONE,XLAMDH(KOFF9),NTOTBE,D2JI(KOFF10),NTOTJ, 4883 * ZERO,WORK(KSCR1),NTOTBE) 4884C 4885 KOFF11 = ILRHSI(ISYMAL) + 1 4886 KOFF12 = IAODIS(ISYMAL,ISYMBE) + 1 4887C 4888 NTOTAL = MAX(NBAS(ISYMAL),1) 4889 NTOTBE = MAX(NBAS(ISYMBE),1) 4890C 4891 CALL DGEMM('N','T',NBAS(ISYMAL),NBAS(ISYMBE),NRHFFR(ISYMI), 4892 * ONE,CMO(KOFF11),NTOTAL,WORK(KSCR1),NTOTBE,ONE, 4893 * D2AO(KOFF12),NTOTAL) 4894C 4895 120 CONTINUE 4896C 4897C--------------------------------- 4898C Do the virtual-frozen block. 4899C--------------------------------- 4900C 4901 DO 130 ISYMAL = 1,NSYM 4902C 4903 ISYMBE = MULD2H(ISYMPQ,ISYMAL) 4904 ISYMA = ISYMAL 4905 ISYMI = ISYMBE 4906C 4907 IF (NRHFFR(ISYMI) .EQ. 0) GOTO 130 4908C 4909 KSCR1 = 1 4910 KEND1 = KSCR1 + NBAS(ISYMAL)*NRHFFR(ISYMI) 4911 LWRK1 = LWORK - KEND1 4912C 4913 IF (LWRK1 .LT. 0) THEN 4914 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 4915 CALL QUIT('Insufficient work space in CC_FD2AO') 4916 ENDIF 4917C 4918 CALL DZERO(WORK(KSCR1),KEND1) 4919C 4920 KOFF13 = IGLMVI(ISYMAL,ISYMA) + 1 4921 KOFF14 = IT1FRO(ISYMA,ISYMI) + 1 4922C 4923 NTOTAL = MAX(NBAS(ISYMAL),1) 4924 NTOTA = MAX(NVIR(ISYMA),1) 4925C 4926 CALL DGEMM('N','N',NBAS(ISYMAL),NRHFFR(ISYMI),NVIR(ISYMA), 4927 * ONE,XLAMDP(KOFF13),NTOTAL,D2AI(KOFF14),NTOTA, 4928 * ZERO,WORK(KSCR1),NTOTAL) 4929C 4930 KOFF15 = ILRHSI(ISYMBE) + 1 4931 KOFF16 = IAODIS(ISYMAL,ISYMBE) + 1 4932C 4933 NTOTAL = MAX(NBAS(ISYMAL),1) 4934 NTOTBE = MAX(NBAS(ISYMBE),1) 4935C 4936 CALL DGEMM('N','T',NBAS(ISYMAL),NBAS(ISYMBE),NRHFFR(ISYMI), 4937 * ONE,WORK(KSCR1),NTOTAL,CMO(KOFF15),NTOTBE,ONE, 4938 * D2AO(KOFF16),NTOTAL) 4939C 4940 130 CONTINUE 4941C 4942C--------------------------------- 4943C Do the frozen-virtual block. 4944C--------------------------------- 4945C 4946 DO 140 ISYMAL = 1,NSYM 4947C 4948 ISYMBE = MULD2H(ISYMPQ,ISYMAL) 4949 ISYMI = ISYMAL 4950 ISYMA = ISYMBE 4951C 4952 IF (NRHFFR(ISYMI) .EQ. 0) GOTO 140 4953C 4954 KSCR1 = 1 4955 KEND1 = KSCR1 + NBAS(ISYMBE)*NRHFFR(ISYMI) 4956 LWRK1 = LWORK - KEND1 4957C 4958 IF (LWRK1 .LT. 0) THEN 4959 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 4960 CALL QUIT('Insufficient work space in CC_FD2AO') 4961 ENDIF 4962C 4963 CALL DZERO(WORK(KSCR1),KEND1) 4964C 4965 KOFF17 = IGLMVI(ISYMBE,ISYMA) + 1 4966 KOFF18 = IT1FRO(ISYMA,ISYMI) + 1 4967C 4968 NTOTBE = MAX(NBAS(ISYMBE),1) 4969 NTOTA = MAX(NVIR(ISYMA),1) 4970C 4971 CALL DGEMM('N','N',NBAS(ISYMBE),NRHFFR(ISYMI),NVIR(ISYMA), 4972 * ONE,XLAMDH(KOFF17),NTOTBE,D2IA(KOFF18),NTOTA, 4973 * ZERO,WORK(KSCR1),NTOTBE) 4974C 4975 KOFF19 = ILRHSI(ISYMAL) + 1 4976 KOFF20 = IAODIS(ISYMAL,ISYMBE) + 1 4977C 4978 NTOTAL = MAX(NBAS(ISYMAL),1) 4979 NTOTBE = MAX(NBAS(ISYMBE),1) 4980C 4981 CALL DGEMM('N','T',NBAS(ISYMAL),NBAS(ISYMBE),NRHFFR(ISYMI), 4982 * ONE,CMO(KOFF19),NTOTAL,WORK(KSCR1),NTOTBE,ONE, 4983 * D2AO(KOFF20),NTOTAL) 4984C 4985 140 CONTINUE 4986C 4987 CALL QEXIT('CC_FD2AO') 4988C 4989 RETURN 4990 END 4991C /* Deck cc_fd2bl */ 4992 SUBROUTINE CC_FD2BL(D2II,D2IJ,D2JI,D2AI,D2IA,DIJ,DAB,DAI,DIA, 4993 * CMO,XLAMDH,XLAMDP,WORK,LWORK,IDEL,ISYMD, 4994 * G,ISYMG) 4995C 4996C Written by Asger Halkier 12/8 - 1998 4997C 4998C Purpose: To calculate the contributions to d(pq,gam;del) 4999C where at least one of the two indices p & q is frozen. 5000C 5001#include "implicit.h" 5002 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 5003 DIMENSION D2II(*), D2IJ(*), D2JI(*), D2AI(*), D2IA(*), DIJ(*) 5004 DIMENSION DAB(*), DAI(*), DIA(*), CMO(*), XLAMDH(*) 5005 DIMENSION XLAMDP(*), WORK(LWORK) 5006#include "priunit.h" 5007#include "ccorb.h" 5008#include "ccsdsym.h" 5009#include "ccfro.h" 5010C 5011 CALL QENTER('CC_FD2BL') 5012C 5013C------------------------------- 5014C Work space allocation one. 5015C------------------------------- 5016C 5017 ISYMB = ISYMD 5018 ISYMK = ISYMD 5019 ISYMA = ISYMB 5020 ISYMI = ISYMG 5021C 5022 KVECA = 1 5023 KVECB = KVECA + NVIR(ISYMA) 5024 KVECK = KVECB + NVIR(ISYMB) 5025 KEND1 = KVECK + NRHF(ISYMK) 5026 LWRK1 = LWORK - KEND1 5027C 5028 IF (LWRK1 .LT. 0) THEN 5029 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 5030 CALL QUIT('Insufficient work space in CC_FD2BL') 5031 ENDIF 5032C 5033 CALL DZERO(WORK(KVECA),KEND1) 5034C 5035 KOFF1 = IGLMVI(ISYMD,ISYMB) + IDEL - IBAS(ISYMD) 5036 KOFF2 = IGLMRH(ISYMD,ISYMK) + IDEL - IBAS(ISYMD) 5037 CALL DCOPY(NVIR(ISYMB),XLAMDH(KOFF1),NBAS(ISYMD),WORK(KVECB),1) 5038 CALL DCOPY(NRHF(ISYMK),XLAMDH(KOFF2),NBAS(ISYMD),WORK(KVECK),1) 5039C 5040C-------------------------------------- 5041C Calculate intermediate vector Va. 5042C-------------------------------------- 5043C 5044 KOFF3 = IMATAB(ISYMA,ISYMB) + 1 5045 KOFF4 = IT1AM(ISYMA,ISYMK) + 1 5046C 5047 NTOTA = MAX(NVIR(ISYMA),1) 5048C 5049 CALL DGEMV('N',NVIR(ISYMA),NVIR(ISYMB),ONE,DAB(KOFF3),NTOTA, 5050 * WORK(KVECB),1,ZERO,WORK(KVECA),1) 5051C 5052 CALL DGEMV('N',NVIR(ISYMA),NRHF(ISYMK),ONE,DAI(KOFF4),NTOTA, 5053 * WORK(KVECK),1,ONE,WORK(KVECA),1) 5054C 5055C------------------------------------ 5056C Calculate virtual-frozen block. 5057C------------------------------------ 5058C 5059 DO 100 I = 1,NRHFFR(ISYMI) 5060C 5061 KOFF5 = ILRHSI(ISYMG) + NBAS(ISYMG)*(I - 1) + G 5062 KOFF6 = IT1FRO(ISYMA,ISYMI) + NVIR(ISYMA)*(I - 1) + 1 5063 CALL DAXPY(NVIR(ISYMA),-CMO(KOFF5),WORK(KVECA),1,D2AI(KOFF6),1) 5064C 5065 100 CONTINUE 5066C 5067C----------------------------------------------------- 5068C Add contribution to frozen-frozen block from Va. 5069C----------------------------------------------------- 5070C 5071 IF (ISYMG .EQ. ISYMD) THEN 5072C 5073 KOFF7 = IGLMVI(ISYMG,ISYMA) + G 5074 FAC = DDOT(NVIR(ISYMA),WORK(KVECA),1,XLAMDP(KOFF7),NBAS(ISYMG)) 5075C 5076 DO 110 ISYMJ = 1,NSYM 5077 DO 120 J = 1,NRHFFR(ISYMJ) 5078 KOFF8 = IFROFR(ISYMJ,ISYMJ) + NRHFFR(ISYMJ)*(J - 1) + J 5079 D2II(KOFF8) = D2II(KOFF8) + TWO*FAC 5080 120 CONTINUE 5081 110 CONTINUE 5082 ENDIF 5083C 5084C------------------------------- 5085C Work space allocation two. 5086C------------------------------- 5087C 5088 ISYMA = ISYMD 5089 ISYML = ISYMD 5090 ISYMI = ISYMA 5091 ISYMJ = ISYMG 5092C 5093 KVECI = 1 5094 KVECA = KVECI + NRHF(ISYMI) 5095 KVECL = KVECA + NVIR(ISYMA) 5096 KEND1 = KVECL + NRHF(ISYML) 5097 LWRK1 = LWORK - KEND1 5098C 5099 IF (LWRK1 .LT. 0) THEN 5100 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 5101 CALL QUIT('Insufficient work space in CC_FD2BL') 5102 ENDIF 5103C 5104 CALL DZERO(WORK(KVECI),KEND1) 5105C 5106 KOFF9 = IGLMVI(ISYMD,ISYMA) + IDEL - IBAS(ISYMD) 5107 KOFF10 = IGLMRH(ISYMD,ISYML) + IDEL - IBAS(ISYMD) 5108 CALL DCOPY(NVIR(ISYMA),XLAMDH(KOFF9),NBAS(ISYMD),WORK(KVECA),1) 5109 CALL DCOPY(NRHF(ISYML),XLAMDH(KOFF10),NBAS(ISYMD),WORK(KVECL),1) 5110C 5111C-------------------------------------- 5112C Calculate intermediate vector Vi. 5113C-------------------------------------- 5114C 5115 KOFF11 = IT1AM(ISYMA,ISYMI) + 1 5116 KOFF12 = IMATIJ(ISYMI,ISYML) + 1 5117C 5118 NTOTA = MAX(NVIR(ISYMA),1) 5119 NTOTI = MAX(NRHF(ISYMI),1) 5120C 5121 CALL DGEMV('T',NVIR(ISYMA),NRHF(ISYMI),ONE,DIA(KOFF11),NTOTA, 5122 * WORK(KVECA),1,ZERO,WORK(KVECI),1) 5123C 5124 CALL DGEMV('N',NRHF(ISYMI),NRHF(ISYML),ONE,DIJ(KOFF12),NTOTI, 5125 * WORK(KVECL),1,ONE,WORK(KVECI),1) 5126C 5127C--------------------------------------- 5128C Calculate correlated-frozen block. 5129C--------------------------------------- 5130C 5131 DO 130 J = 1,NRHFFR(ISYMJ) 5132C 5133 KOFF13 = ILRHSI(ISYMG) + NBAS(ISYMG)*(J - 1) + G 5134 KOFF14 = ICOFRO(ISYMI,ISYMJ) + NRHF(ISYMI)*(J - 1) + 1 5135 CALL DAXPY(NRHF(ISYMI),-CMO(KOFF13),WORK(KVECI),1, 5136 * D2IJ(KOFF14),1) 5137C 5138 130 CONTINUE 5139C 5140C----------------------------------------------------- 5141C Add contribution to frozen-frozen block from Va. 5142C----------------------------------------------------- 5143C 5144 IF (ISYMG .EQ. ISYMD) THEN 5145C 5146 ISYMK = ISYMD 5147 KVECK = KVECI 5148C 5149 KOFF15 = IGLMRH(ISYMG,ISYMK) + G 5150 FAC = DDOT(NRHF(ISYMK),WORK(KVECK),1,XLAMDP(KOFF15), 5151 * NBAS(ISYMG)) 5152C 5153 DO 140 ISYMM = 1,NSYM 5154 DO 150 M = 1,NRHFFR(ISYMM) 5155 KOFF16 = IFROFR(ISYMM,ISYMM) + NRHFFR(ISYMM)*(M - 1) + M 5156 D2II(KOFF16) = D2II(KOFF16) + TWO*FAC 5157 150 CONTINUE 5158 140 CONTINUE 5159 ENDIF 5160C 5161C---------------------------------------------------------------------- 5162C Calculate Hartree-Fock like contributions to frozen-frozen block. 5163C---------------------------------------------------------------------- 5164C 5165 IF (ISYMG .EQ. ISYMD) THEN 5166C 5167 KOFF17 = ILRHSI(ISYMG) + G 5168 KOFF18 = ILRHSI(ISYMD) + IDEL - IBAS(ISYMD) 5169 FAC = TWO*DDOT(NRHFFR(ISYMG),CMO(KOFF17),NBAS(ISYMG), 5170 * CMO(KOFF18),NBAS(ISYMD)) 5171C 5172 DO 160 ISYMI = 1,NSYM 5173 DO 170 I = 1,NRHFFR(ISYMI) 5174 KOFF19 = IFROFR(ISYMI,ISYMI) + NRHFFR(ISYMI)*(I - 1) + I 5175 D2II(KOFF19) = D2II(KOFF19) + TWO*FAC 5176 170 CONTINUE 5177 160 CONTINUE 5178 ENDIF 5179C 5180 ISYMI = ISYMD 5181 ISYMJ = ISYMG 5182 D = IDEL - IBAS(ISYMD) 5183C 5184 DO 180 J = 1,NRHFFR(ISYMJ) 5185 DO 190 I = 1,NRHFFR(ISYMI) 5186 KOFF20 = IFROFR(ISYMI,ISYMJ) + NRHFFR(ISYMI)*(J - 1) + I 5187 KOFF21 = ILRHSI(ISYMD) + NBAS(ISYMD)*(I - 1) + D 5188 KOFF22 = ILRHSI(ISYMG) + NBAS(ISYMG)*(J - 1) + G 5189 D2II(KOFF20) = D2II(KOFF20) - TWO*CMO(KOFF21)*CMO(KOFF22) 5190 190 CONTINUE 5191 180 CONTINUE 5192C 5193C--------------------------------- 5194C Work space allocation three. 5195C--------------------------------- 5196C 5197 ISYMB = ISYMG 5198 ISYMJ = ISYMG 5199 ISYMA = ISYMB 5200 ISYMI = ISYMD 5201C 5202 KVECA = 1 5203 KVECB = KVECA + NVIR(ISYMA) 5204 KVECJ = KVECB + NVIR(ISYMB) 5205 KEND1 = KVECJ + NRHF(ISYMJ) 5206 LWRK1 = LWORK - KEND1 5207C 5208 IF (LWRK1 .LT. 0) THEN 5209 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 5210 CALL QUIT('Insufficient work space in CC_FD2BL') 5211 ENDIF 5212C 5213 CALL DZERO(WORK(KVECI),KEND1) 5214C 5215 KOFF23 = IGLMVI(ISYMG,ISYMB) + G 5216 KOFF24 = IGLMRH(ISYMG,ISYMJ) + G 5217 CALL DCOPY(NVIR(ISYMB),XLAMDP(KOFF23),NBAS(ISYMG),WORK(KVECB),1) 5218 CALL DCOPY(NRHF(ISYMJ),XLAMDP(KOFF24),NBAS(ISYMG),WORK(KVECJ),1) 5219C 5220C-------------------------------------- 5221C Calculate intermediate vector Wa. 5222C-------------------------------------- 5223C 5224 KOFF25 = IMATAB(ISYMB,ISYMA) + 1 5225 KOFF26 = IT1AM(ISYMA,ISYMJ) + 1 5226C 5227 NTOTB = MAX(NVIR(ISYMB),1) 5228 NTOTA = MAX(NVIR(ISYMA),1) 5229C 5230 CALL DGEMV('T',NVIR(ISYMB),NVIR(ISYMA),ONE,DAB(KOFF25),NTOTB, 5231 * WORK(KVECB),1,ZERO,WORK(KVECA),1) 5232C 5233 CALL DGEMV('N',NVIR(ISYMA),NRHF(ISYMJ),ONE,DIA(KOFF26),NTOTA, 5234 * WORK(KVECJ),1,ONE,WORK(KVECA),1) 5235C 5236C------------------------------------ 5237C Calculate frozen-virtual block. 5238C------------------------------------ 5239C 5240 DO 200 I = 1,NRHFFR(ISYMI) 5241C 5242 KOFF27 = ILRHSI(ISYMD) + NBAS(ISYMD)*(I - 1) 5243 * + IDEL - IBAS(ISYMD) 5244 KOFF28 = IT1FRO(ISYMA,ISYMI) + NVIR(ISYMA)*(I - 1) + 1 5245 CALL DAXPY(NVIR(ISYMA),-CMO(KOFF27),WORK(KVECA),1, 5246 * D2IA(KOFF28),1) 5247C 5248 200 CONTINUE 5249C 5250C-------------------------------- 5251C Work space allocation four. 5252C-------------------------------- 5253C 5254 ISYMA = ISYMG 5255 ISYMK = ISYMG 5256 ISYMJ = ISYMA 5257 ISYMI = ISYMD 5258C 5259 KVECJ = 1 5260 KVECA = KVECJ + NRHF(ISYMJ) 5261 KVECK = KVECA + NVIR(ISYMA) 5262 KEND1 = KVECK + NRHF(ISYMK) 5263 LWRK1 = LWORK - KEND1 5264C 5265 IF (LWRK1 .LT. 0) THEN 5266 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 5267 CALL QUIT('Insufficient work space in CC_FD2BL') 5268 ENDIF 5269C 5270 CALL DZERO(WORK(KVECI),KEND1) 5271C 5272 KOFF29 = IGLMVI(ISYMG,ISYMA) + G 5273 KOFF30 = IGLMRH(ISYMG,ISYMk) + G 5274 CALL DCOPY(NVIR(ISYMA),XLAMDP(KOFF29),NBAS(ISYMG),WORK(KVECA),1) 5275 CALL DCOPY(NRHF(ISYMK),XLAMDP(KOFF30),NBAS(ISYMG),WORK(KVECK),1) 5276C 5277C-------------------------------------- 5278C Calculate intermediate vector Wj. 5279C-------------------------------------- 5280C 5281 KOFF31 = IT1AM(ISYMA,ISYMJ) + 1 5282 KOFF32 = IMATIJ(ISYMK,ISYMJ) + 1 5283C 5284 NTOTA = MAX(NVIR(ISYMA),1) 5285 NTOTK = MAX(NRHF(ISYMK),1) 5286C 5287 CALL DGEMV('T',NVIR(ISYMA),NRHF(ISYMJ),ONE,DAI(KOFF31),NTOTA, 5288 * WORK(KVECA),1,ZERO,WORK(KVECJ),1) 5289C 5290 CALL DGEMV('T',NRHF(ISYMK),NRHF(ISYMJ),ONE,DIJ(KOFF32),NTOTK, 5291 * WORK(KVECK),1,ONE,WORK(KVECJ),1) 5292C 5293C--------------------------------------- 5294C Calculate frozen-correlated block. 5295C--------------------------------------- 5296C 5297 DO 210 I = 1,NRHFFR(ISYMI) 5298C 5299 KOFF33 = ILRHSI(ISYMD) + NBAS(ISYMD)*(I - 1) 5300 * + IDEL - IBAS(ISYMD) 5301 KOFF34 = ICOFRO(ISYMJ,ISYMI) + NRHF(ISYMJ)*(I - 1) + 1 5302 CALL DAXPY(NRHF(ISYMJ),-CMO(KOFF33),WORK(KVECJ),1, 5303 * D2JI(KOFF34),1) 5304C 5305 210 CONTINUE 5306C 5307 CALL QEXIT('CC_FD2BL') 5308C 5309 RETURN 5310 END 5311C /* Deck ccs_den2 */ 5312 SUBROUTINE CCS_DEN2(D2IJG,CMO,WORK,LWORK,IDEL,ISYMD) 5313C 5314C Written by Asger Halkier 4/6 - 1998. 5315C 5316C Purpose: Calculate the 2 electron CCS (=SCF) density 5317C d(pq,gam;del) for a given delta (IDEL). 5318C 5319#include "implicit.h" 5320 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 5321 DIMENSION D2IJG(*), CMO(*), WORK(LWORK) 5322#include "priunit.h" 5323#include "ccorb.h" 5324#include "ccsdsym.h" 5325#include "ccfro.h" 5326C 5327 CALL QENTER('CCS_DEN2') 5328C 5329 ISYDEN = ISYMD 5330C 5331C--------------------------------------------------- 5332C Calculate the (occ.occ,occ;del) density block. 5333C--------------------------------------------------- 5334C 5335 KD2IJK = 1 5336 KEND1 = KD2IJK + NFRIJK(ISYDEN) 5337 LWRK1 = LWORK - KEND1 5338C 5339 IF (LWRK1 .LT. 0) THEN 5340 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 5341 CALL QUIT('Insufficient space for allocation in CCS_DEN2') 5342 ENDIF 5343C 5344 CALL DZERO(WORK(KD2IJK),NFRIJK(ISYDEN)) 5345C 5346C--------------------------------- 5347C Calculate the contributions. 5348C--------------------------------- 5349C 5350 CALL FIJK01(WORK(KD2IJK),CMO,IDEL,ISYMD) 5351 CALL FIJK02(WORK(KD2IJK),CMO,IDEL,ISYMD) 5352C 5353 CALL F2_PQAO(D2IJG,WORK(KD2IJK),CMO,ISYDEN) 5354C 5355 CALL QEXIT('CCS_DEN2') 5356C 5357 RETURN 5358 END 5359C /* Deck mp2_den2 */ 5360 SUBROUTINE MP2_DEN2(D2IAG,T2AMM,CMO,WORK,LWORK,IDEL,ISYMD) 5361C 5362C Written by Asger Halkier 5/6 - 1998. 5363C 5364C Purpose: Calculate the non-SCF part of the 2 electron MP2 5365C density d(pq,gam;del) for a given delta (IDEL). 5366C 5367#include "implicit.h" 5368 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 5369 DIMENSION D2IAG(*), T2AMM(*), CMO(*), WORK(LWORK) 5370#include "priunit.h" 5371#include "ccorb.h" 5372#include "ccsdsym.h" 5373C 5374 CALL QENTER('MP2_DEN2') 5375C 5376 ISYDEN = ISYMD 5377C 5378C---------------------------------------------------- 5379C Calculate terms of the (occ.vir,occ;del) block. 5380C Note that we exploit the full permutational 5381C symmetry of the 2-electron integrals later on. 5382C---------------------------------------------------- 5383C 5384 KD2IAJ = 1 5385 KEND1 = KD2IAJ + NT2BCD(ISYDEN) 5386 LWRK1 = LWORK - KEND1 5387C 5388 IF (LWRK1 .LT. 0) THEN 5389 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 5390 CALL QUIT( 5391 * 'Insufficient space for fourth allocation in MP2_DEN2') 5392 ENDIF 5393C 5394 CALL DZERO(WORK(KD2IAJ),NT2BCD(ISYDEN)) 5395C 5396C--------------------------------- 5397C Calculate the contributions. 5398C--------------------------------- 5399C 5400 ISYMT2 = 1 5401 CALL CC_TI(WORK(KD2IAJ),ISYMD,T2AMM,ISYMT2,CMO,1, 5402 * WORK(KEND1),LWRK1,IDEL,ISYMD) 5403C 5404 ICON = 2 5405 CALL CCD2_PQAO(D2IAG,ISYDEN,WORK(KD2IAJ),ISYDEN,CMO,1,ICON) 5406C 5407 CALL QEXIT('MP2_DEN2') 5408C 5409 RETURN 5410 END 5411C /* Deck ccmp_dao */ 5412 SUBROUTINE CCMP_DAO(AODEN,DENIJ,DENIA,CMO,XCMO,WORK,LWORK,ISDEN) 5413C 5414C Written by Asger Halkier 8/6 - 1998 5415C 5416C Version: 1.0 5417C 5418C Purpose: To transform the blocks of the CCS (1 block) and 5419C MP2 (2 blocks) one electron density to AO-basis 5420C and add the results! 5421C 5422#include "implicit.h" 5423 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 5424 DIMENSION AODEN(*), DENIJ(*), DENIA(*) 5425 DIMENSION CMO(*), XCMO(*), WORK(LWORK) 5426#include "priunit.h" 5427#include "ccorb.h" 5428#include "ccsdsym.h" 5429#include "ccsdinp.h" 5430#include "ccfro.h" 5431C 5432 CALL QENTER('CCMP_DAO') 5433C 5434 DO 100 ISYM1 = 1,NSYM 5435C 5436 ISYM2 = MULD2H(ISDEN,ISYM1) 5437C 5438 LNEED = NBAS(ISYM1)*NRHFS(ISYM2) 5439C 5440 IF (LWORK .LT. LNEED) THEN 5441 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', LNEED 5442 CALL QUIT('Insufficient work space in CCMP_DAO') 5443 ENDIF 5444C 5445 CALL DZERO(WORK,LNEED) 5446C 5447C------------------------------------------ 5448C Transform occupied-occupied block. 5449C------------------------------------------ 5450C 5451 KOFF1 = ILRHSI(ISYM1) + 1 5452 KOFF2 = IFROIJ(ISYM1,ISYM2) + 1 5453C 5454 NTOBA1 = MAX(NBAS(ISYM1),1) 5455 NTORH1 = MAX(NRHFS(ISYM1),1) 5456C 5457 CALL DGEMM('N','N',NBAS(ISYM1),NRHFS(ISYM2),NRHFS(ISYM1), 5458 * ONE,CMO(KOFF1),NTOBA1,DENIJ(KOFF2),NTORH1, 5459 * ZERO,WORK,NTOBA1) 5460C 5461 KOFF3 = ILRHSI(ISYM2) + 1 5462 KOFF4 = IAODIS(ISYM1,ISYM2) + 1 5463C 5464 NTOBA1 = MAX(NBAS(ISYM1),1) 5465 NTOBA2 = MAX(NBAS(ISYM2),1) 5466C 5467 CALL DGEMM('N','T',NBAS(ISYM1),NBAS(ISYM2),NRHFS(ISYM2), 5468 * ONE,WORK,NTOBA1,CMO(KOFF3),NTOBA2,ONE, 5469 * AODEN(KOFF4),NTOBA1) 5470C 5471 100 CONTINUE 5472C 5473 IF (MP2) THEN 5474 DO 110 ISYM1 = 1,NSYM 5475C 5476 ISYM2 = MULD2H(ISDEN,ISYM1) 5477C 5478 LNEED = NBAS(ISYM2)*NRHF(ISYM1) 5479C 5480 IF (LWORK .LT. LNEED) THEN 5481 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', LNEED 5482 CALL QUIT('Insufficient work space in CCMP_DAO') 5483 ENDIF 5484C 5485 CALL DZERO(WORK,LNEED) 5486C 5487C------------------------------------------------------------- 5488C Transform occupied-virtual block (stored transposed). 5489C------------------------------------------------------------- 5490C 5491 KOFF1 = IGLMVI(ISYM2,ISYM2) + 1 5492 KOFF2 = IT1AM(ISYM2,ISYM1) + 1 5493C 5494 NTOBA2 = MAX(NBAS(ISYM2),1) 5495 NTOVI2 = MAX(NVIR(ISYM2),1) 5496C 5497 CALL DGEMM('N','N',NBAS(ISYM2),NRHF(ISYM1),NVIR(ISYM2), 5498 * ONE,XCMO(KOFF1),NTOBA2,DENIA(KOFF2),NTOVI2, 5499 * ZERO,WORK,NTOBA2) 5500C 5501 KOFF3 = IGLMRH(ISYM1,ISYM1) + 1 5502 KOFF4 = IAODIS(ISYM1,ISYM2) + 1 5503C 5504 NTOBA1 = MAX(NBAS(ISYM1),1) 5505 NTOBA2 = MAX(NBAS(ISYM2),1) 5506C 5507 CALL DGEMM('N','T',NBAS(ISYM1),NBAS(ISYM2),NRHF(ISYM1), 5508 * ONE,XCMO(KOFF3),NTOBA1,WORK,NTOBA2,ONE, 5509 * AODEN(KOFF4),NTOBA1) 5510C 5511 110 CONTINUE 5512 ENDIF 5513C 5514 CALL QEXIT('CCMP_DAO') 5515C 5516 RETURN 5517 END 5518C /* Deck fijk01 */ 5519 SUBROUTINE FIJK01(D2IJK,CMO,IDEL,ISYMD) 5520C 5521C Written by Asger Halkier 11/6 - 1998 5522C 5523C Purpose: To calculate the first contribution to the 5524C CCS 2-electron density including frozen-core 5525C contributions - based on DIJK01. 5526C 5527#include "implicit.h" 5528 PARAMETER(FOUR = 4.0D0) 5529 DIMENSION D2IJK(*), CMO(*) 5530#include "priunit.h" 5531#include "ccorb.h" 5532#include "ccsdsym.h" 5533#include "ccfro.h" 5534C 5535 CALL QENTER('FIJK01') 5536C 5537 ISYDEN = ISYMD 5538 ISYMK = ISYMD 5539 ISYMIJ = MULD2H(ISYDEN,ISYMK) 5540C 5541C------------------------------------------ 5542C Calculate adresses and add to result. 5543C------------------------------------------ 5544C 5545 DO 100 K = 1,NRHFS(ISYMK) 5546C 5547 DO 110 ISYMJ = 1,NSYM 5548C 5549 DO 120 J = 1,NRHFS(ISYMJ) 5550C 5551 NIJK = IFRIJK(ISYMIJ,ISYMK) + NFROIJ(ISYMIJ)*(K - 1) 5552 * + IFROIJ(ISYMJ,ISYMJ) + NRHFS(ISYMJ)*(J - 1) + J 5553 NDEL = ILRHSI(ISYMK) + NBAS(ISYMD)*(K - 1) 5554 * + IDEL - IBAS(ISYMD) 5555C 5556 D2IJK(NIJK) = D2IJK(NIJK) + FOUR*CMO(NDEL) 5557C 5558 120 CONTINUE 5559 110 CONTINUE 5560 100 CONTINUE 5561C 5562 CALL QEXIT('FIJK01') 5563C 5564 RETURN 5565 END 5566C /* Deck fijk02 */ 5567 SUBROUTINE FIJK02(D2IJK,CMO,IDEL,ISYMD) 5568C 5569C Written by Asger Halkier 11/6 - 1998 5570C 5571C Purpose: To calculate the second contribution to the 5572C CCS 2-electron density including frozen-core 5573C contributions - based on DIJK02. 5574C 5575#include "implicit.h" 5576 PARAMETER(TWO = 2.0D0) 5577 DIMENSION D2IJK(*), CMO(*) 5578#include "priunit.h" 5579#include "ccorb.h" 5580#include "ccsdsym.h" 5581#include "ccfro.h" 5582C 5583 CALL QENTER('FIJK02') 5584C 5585 ISYDEN = ISYMD 5586 ISYMI = ISYMD 5587C 5588C------------------------------------------ 5589C Calculate adresses and add to result. 5590C------------------------------------------ 5591C 5592 DO 100 ISYMK = 1,NSYM 5593C 5594 ISYMIJ = MULD2H(ISYMI,ISYMK) 5595C 5596 DO 110 K = 1,NRHFS(ISYMK) 5597C 5598 DO 120 I = 1,NRHFS(ISYMI) 5599C 5600 NIJK = IFRIJK(ISYMIJ,ISYMK) + NFROIJ(ISYMIJ)*(K - 1) 5601 * + IFROIJ(ISYMI,ISYMK) + NRHFS(ISYMI)*(K - 1) + I 5602 NDEL = ILRHSI(ISYMI) + NBAS(ISYMD)*(I - 1) 5603 * + IDEL - IBAS(ISYMD) 5604C 5605 D2IJK(NIJK) = D2IJK(NIJK) - TWO*CMO(NDEL) 5606C 5607 120 CONTINUE 5608 110 CONTINUE 5609 100 CONTINUE 5610C 5611 CALL QEXIT('FIJK02') 5612C 5613 RETURN 5614 END 5615C /* Deck f2_pqao */ 5616 SUBROUTINE F2_PQAO(D2IJG,D2IJK,CMO,ISYDEN) 5617C 5618C Written by Asger Halkier 11/6 - 1998 5619C 5620C Purpose: To calculate the backtransformation of the CCS 5621C two-electron density d(ij,k;del) to d(ij,gam;del) 5622C including the contributions from frozen core orbitals. 5623C 5624#include "implicit.h" 5625 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 5626 DIMENSION D2IJG(*), D2IJK(*), CMO(*) 5627#include "priunit.h" 5628#include "ccorb.h" 5629#include "ccsdsym.h" 5630#include "ccfro.h" 5631C 5632 CALL QENTER('F2_PQAO') 5633C 5634C---------------------------------- 5635C Calculate the transformation. 5636C---------------------------------- 5637C 5638 DO 100 ISYMIJ = 1,NSYM 5639C 5640 ISYMK = MULD2H(ISYMIJ,ISYDEN) 5641 ISYMG = ISYMK 5642C 5643 KOFF1 = IFRIJK(ISYMIJ,ISYMK) + 1 5644 KOFF2 = ILRHSI(ISYMK) + 1 5645 KOFF3 = IF2IJG(ISYMIJ,ISYMG) + 1 5646C 5647 NTOTIJ = MAX(NFROIJ(ISYMIJ),1) 5648 NTOTG = MAX(NBAS(ISYMG),1) 5649C 5650 CALL DGEMM('N','T',NFROIJ(ISYMIJ),NBAS(ISYMG),NRHFS(ISYMK), 5651 * ONE,D2IJK(KOFF1),NTOTIJ,CMO(KOFF2),NTOTG, 5652 * ONE,D2IJG(KOFF3),NTOTIJ) 5653C 5654 100 CONTINUE 5655C 5656 CALL QEXIT('F2_PQAO') 5657C 5658 RETURN 5659 END 5660C /* Deck cc_d2HFao */ 5661 SUBROUTINE CC_D2HFAO(AODEN,DENIJ,CMO,WORK,LWORK,ISDEN) 5662C 5663C Written by Sonia Coriani 20/11/2008 5664C Used to obtain the correlation gradient in Incremental Scheme 5665C Version: 1.0, based on CCMP_DAO 5666C 5667C Purpose: To transform the blocks of the CCS/HF (1 block) 5668C two electron density (IJ_GAMMA_DELTA) to AO-basis 5669C and subtract the results! 5670C 5671#include "implicit.h" 5672 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 5673 DIMENSION AODEN(*), DENIJ(*) 5674 DIMENSION CMO(*), WORK(LWORK) 5675#include "priunit.h" 5676#include "ccorb.h" 5677#include "ccsdsym.h" 5678#include "ccsdinp.h" 5679#include "ccfro.h" 5680 CALL QENTER('CC_D2HFAO') 5681C 5682 DO 100 ISYM1 = 1,NSYM 5683C 5684 ISYM2 = MULD2H(ISDEN,ISYM1) 5685C 5686 LNEED = NBAS(ISYM1)*NRHFS(ISYM2) 5687C 5688 IF (LWORK .LT. LNEED) THEN 5689 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', LNEED 5690 CALL QUIT('Insufficient work space in CCMP_DAO') 5691 ENDIF 5692C 5693 CALL DZERO(WORK,LNEED) 5694C 5695C------------------------------------------ 5696C Transform occupied-occupied block. 5697C------------------------------------------ 5698C 5699 KOFF1 = ILRHSI(ISYM1) + 1 5700 KOFF2 = IFROIJ(ISYM1,ISYM2) + 1 5701C 5702 NTOBA1 = MAX(NBAS(ISYM1),1) 5703 NTORH1 = MAX(NRHFS(ISYM1),1) 5704C 5705 CALL DGEMM('N','N',NBAS(ISYM1),NRHFS(ISYM2),NRHFS(ISYM1), 5706 * ONE,CMO(KOFF1),NTOBA1,DENIJ(KOFF2),NTORH1, 5707 * ZERO,WORK,NTOBA1) 5708C 5709 KOFF3 = ILRHSI(ISYM2) + 1 5710 KOFF4 = IAODIS(ISYM1,ISYM2) + 1 5711C 5712 NTOBA1 = MAX(NBAS(ISYM1),1) 5713 NTOBA2 = MAX(NBAS(ISYM2),1) 5714C 5715 CALL DGEMM('N','T',NBAS(ISYM1),NBAS(ISYM2),NRHFS(ISYM2), 5716 * -ONE,WORK,NTOBA1,CMO(KOFF3),NTOBA2,ONE, 5717 * AODEN(KOFF4),NTOBA1) 5718C 5719 100 CONTINUE 5720C 5721 CALL QEXIT('CC_D2HFAO') 5722C 5723 RETURN 5724 END 5725C /* Deck cc_protected_division */ 5726 Real*8 Function cc_protected_division(nominator,denominator, 5727 * epsn,epsd) 5728C 5729C Thomas Bondo Pedersen, Jan. 2013. 5730C 5731C Compute 5732C 5733C x=nominator/denominator 5734C 5735C If |nominator|<epsn: 5736C x=0 5737C Else: 5738C If |denominator|<epsd: 5739C quit 5740C Else: 5741C x=nominator/denominator 5742C 5743 Implicit None 5744 Real*8 nominator 5745 Real*8 denominator 5746 Real*8 epsn 5747 Real*8 epsd 5748#include "priunit.h" 5749 5750 Real*8 x 5751 5752 If (abs(nominator).lt.epsn) Then 5753 x=0.0d0 5754 Else 5755 If (abs(denominator).lt.epsd) Then 5756 Write(lupri,'(A)') 'Singularity: nominator/denominator' 5757 Write(lupri,'(A,1P,D25.16,A,D25.16)') 5758 * 'nominator=',nominator,' denominator=',denominator 5759 Write(lupri,'(A,1P,D25.16,A,D25.16)') 5760 * 'epsn= ',epsn,' epsd= ',epsd 5761 Call Quit('Singularity (division by zero)') 5762 x=0.0d0 5763 Else 5764 x=nominator/denominator 5765 End If 5766 End If 5767 5768 cc_protected_division=x 5769 5770 Return 5771 End 5772 5773