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_ijcb */ 20*=====================================================================* 21 SUBROUTINE CC_IJCB( XINT1, ISY1ALBE, 22 & XINT2, ISY2ALBE, 23 & IDEL, IGAM, 24 & X1IJCB, X1CJIB, 25 & X2IJCB, X2CJIB, 26 & XLAMDP1, XLAMDH1, ISYML1, 27 & XLAMDP2, XLAMDH2, ISYML2, 28 & XLAMDP3, XLAMDH3, ISYML3, 29 & XLAMDP4, XLAMDH4, ISYML4, 30 & WORK, LWORK, 31 & IOPT, LDERIV, LRELAX, LZERO, 32 & LNEWTA, LX2ISQ ) 33*---------------------------------------------------------------------* 34* 35* Purpose: drive the generalized transformation to 36* (ij^|cb) + (ij|c^b) integrals and 37* (cj^|ib) + (c^j|ib) integrals 38* for the two-index (**|gam del) approach 39* assumes three-index arrays XIJCB & XCJIB in core 40* 41* this routine transforms the indices ij and c, the 42* transformation of the delta index to b has to be done 43* from the outside. 44* 45* XINT1 zero order AO integrals 46* XINT2 derivative AO integrals 47* XA1IJCB, XA1CJIB : TA dependent (zero ord) integrals 48* XA2IJCB, XA2CJIB : TA dependent deriv & relax integrals 49* 50* IOPT=0: (ij^|c del) + (ij|c^ del) only 51* 52* IOPT=1: (ij^|c del) + (ij|c^ del) and 53* (cj^|i del) + (c^j|i del) integrals 54* 55* IF LDERIV=.TRUE. transform also the derivative integrals g[1]: 56* with the XLAMD_1 and XLAMD_3 57* IOPT=0: (ij^|c del) + (ij|c^ del) 58* IOPT=1: (ij^|c del) + (ij|c^ del) and (cj^|i del) + (c^j|i del) 59* 60* IF LRELAX=.TRUE. include relaxation contribution to the 61* derivative integrals from the transformation 62* of g[0] with XLAMD_1 * XLAMD_2 * XLAMD_3 * XLAMD_4 63* (or just reorthonormalization if IRELAX = 0) 64* 65* IF LZERO=.FALSE. skip calculation of zero-order integrals 66* IF LX2ISQ = .TRUE. the (al bet| part of X2INT is already full matr. 67* 68* Written by Sonia Coriani, February 1999 69* based on Christof's CC_IAJB 70* Note that no special case is needed for LAO apart LX2ISQ = TRUE 71* 72*=====================================================================* 73#if defined (IMPLICIT_NONE) 74 IMPLICIT NONE 75#else 76# include "implicit.h" 77#endif 78#include "priunit.h" 79#include "ccorb.h" 80#include "ccsdsym.h" 81#include "maxorb.h" 82#include "ccisao.h" 83 84#if defined (SYS_CRAY) 85 REAL ONE, ZERO 86#else 87 DOUBLE PRECISION ONE, ZERO 88#endif 89 PARAMETER (ONE = 1.0d0, ZERO = 0.0d0) 90 91 LOGICAL LDERIV, LRELAX, LZERO, LNEWTA, LX2ISQ 92 INTEGER IDEL, IGAM, ISY1ALBE, ISY2ALBE 93 INTEGER ISYML1, ISYML2, ISYML3, ISYML4, LWORK, IOPT 94 95#if defined (SYS_CRAY) 96 REAL XLAMDP1(*), XLAMDH1(*) 97 REAL XLAMDP2(*), XLAMDH2(*) 98 REAL XLAMDP3(*), XLAMDH3(*) 99 REAL XLAMDP4(*), XLAMDH4(*) 100 REAL XINT1(*), XINT2(*) 101 REAL X1IJCB(*), X1CJIB(*) 102 REAL X2IJCB(*), X2CJIB(*) 103 REAL WORK(LWORK) 104#else 105 DOUBLE PRECISION XLAMDP1(*), XLAMDH1(*) 106 DOUBLE PRECISION XLAMDP2(*), XLAMDH2(*) 107 DOUBLE PRECISION XLAMDP3(*), XLAMDH3(*) 108 DOUBLE PRECISION XLAMDP4(*), XLAMDH4(*) 109 DOUBLE PRECISION XINT1(*), XINT2(*) 110 DOUBLE PRECISION X1IJCB(*), X1CJIB(*) 111 DOUBLE PRECISION X2IJCB(*), X2CJIB(*) 112 DOUBLE PRECISION WORK(LWORK) 113#endif 114 115* local 116 117 INTEGER ISYL11, ISYL12, ISYL13, ISYL14, ISYL23 118 INTEGER ISYDEL, ISYGAM 119 INTEGER KSCR1, KSCR3, KEND1, LWRK1, KDUM 120 PARAMETER (KDUM = + 999 999 999) 121 122*---------------------------------------------------------------------* 123* set some symmetries 124*---------------------------------------------------------------------* 125 126 ISYDEL = ISAO(IDEL) 127 ISYGAM = ISAO(IGAM) 128 129*---------------------------------------------------------------------* 130* work space allocation: 131* 132* KSCR1 -- I^{del,gam}(alp bet) 133* 134* KSCR3 -- I^{del,gam}(alp bet)[1] 135* 136*---------------------------------------------------------------------* 137 KSCR1 = 1 138 KEND1 = KSCR1 + N2BST(ISY1ALBE) 139 140c IF ((LDERIV .OR. LRELAX).AND.(.NOT.LX2ISQ)) THEN 141 142 IF ((LDERIV).AND.(.NOT.LX2ISQ)) THEN 143 KSCR3 = KEND1 144 KEND1 = KSCR3 + N2BST(ISY2ALBE) 145 ELSE 146 KSCR3 = KDUM 147 END IF 148 149 LWRK1 = LWORK - KEND1 150 151 IF ( LWRK1 .LT. 0) THEN 152 CALL QUIT('Insufficient memory in CC_IJCB.') 153 END IF 154 155*---------------------------------------------------------------------* 156* square zero-order integral matrix up (alp and bet) 157* and derivative also if (LDERIV).AND.(.NOT.LX2ISQ) 158*---------------------------------------------------------------------* 159 160 CALL CCSD_SYMSQ(XINT1,ISY1ALBE,WORK(KSCR1)) 161 162 IF ((LDERIV).AND.(.NOT.LX2ISQ)) THEN 163 CALL CCSD_SYMSQ(XINT2,ISY2ALBE,WORK(KSCR3)) 164 END IF 165 166*---------------------------------------------------------------------* 167* call routine for actual calculation of transformed integrals 168*---------------------------------------------------------------------* 169 IF (.NOT.LX2ISQ) THEN 170 CALL CCIJCB(WORK(KSCR1),ISY1ALBE,WORK(KSCR3),ISY2ALBE, 171 & IDEL, IGAM, X1IJCB, X1CJIB, X2IJCB, X2CJIB, 172 & XLAMDP1, XLAMDH1, ISYML1, XLAMDP2, XLAMDH2, ISYML2, 173 & XLAMDP3, XLAMDH3, ISYML3, XLAMDP4, XLAMDH4, ISYML4, 174 & WORK(KEND1), LWRK1, 175 & IOPT, LDERIV, LRELAX, LZERO, LNEWTA ) 176 ELSE 177 CALL CCIJCB(WORK(KSCR1),ISY1ALBE,XINT2,ISY2ALBE, 178 & IDEL, IGAM, X1IJCB, X1CJIB, X2IJCB, X2CJIB, 179 & XLAMDP1, XLAMDH1, ISYML1, XLAMDP2, XLAMDH2, ISYML2, 180 & XLAMDP3, XLAMDH3, ISYML3, XLAMDP4, XLAMDH4, ISYML4, 181 & WORK(KEND1), LWRK1, 182 & IOPT, LDERIV, LRELAX, LZERO, LNEWTA ) 183 END IF 184 185*---------------------------------------------------------------------* 186* return 187*---------------------------------------------------------------------* 188 189 RETURN 190 END 191*=====================================================================* 192* END OF SUBROUTINE CC_IJCB * 193*=====================================================================* 194c /* deck ccijcb */ 195*=====================================================================* 196 SUBROUTINE CCIJCB( X1INT, ISY1ALBE, 197 & X2INT, ISY2ALBE, 198 & IDEL, IGAM, 199 & X1IJCB, X1CJIB, 200 & X2IJCB, X2CJIB, 201 & XLAMDP1, XLAMDH1, ISYML1, 202 & XLAMDP2, XLAMDH2, ISYML2, 203 & XLAMDP3, XLAMDH3, ISYML3, 204 & XLAMDP4, XLAMDH4, ISYML4, 205 & WORK, LWORK, 206 & IOPT, LDERIV, LRELAX, LZERO, 207 & LNEWTA ) 208*---------------------------------------------------------------------* 209* 210* Purpose: realise the generalized transformation to 211* (ij^|cb) + (ij|c^b) integrals and 212* (cj^|ib) + (c^j|ib) integrals 213* for the two-index (**|gam del) approach 214* assumes three-index arrays XIJCB & XCJIB in core 215* 216* See CC_IJCB for details 217* 218* Written by Sonia Coriani, February 1999 219* 220*=====================================================================* 221#if defined (IMPLICIT_NONE) 222 IMPLICIT NONE 223#else 224# include "implicit.h" 225#endif 226#include "priunit.h" 227#include "ccorb.h" 228#include "ccsdsym.h" 229#include "maxorb.h" 230#include "ccisao.h" 231 232#if defined (SYS_CRAY) 233 REAL ONE, ZERO 234#else 235 DOUBLE PRECISION ONE, ZERO 236#endif 237 PARAMETER (ONE = 1.0d0, ZERO = 0.0d0) 238 239 LOGICAL LDERIV, LRELAX, LZERO, LNEWTA 240 INTEGER IDEL, IGAM, ISY1ALBE, ISY2ALBE, ISYML1, ISYML2 241 INTEGER ISYML3, ISYML4, LWORK, IOPT, IDUMMY 242 243#if defined (SYS_CRAY) 244 REAL XLAMDP1(*), XLAMDH1(*) 245 REAL XLAMDP2(*), XLAMDH2(*) 246 REAL XLAMDP3(*), XLAMDH3(*) 247 REAL XLAMDP4(*), XLAMDH4(*) 248 REAL X1INT(*), X2INT(*) 249 REAL X1IJCB(*), X1CJIB(*) 250 REAL X2IJCB(*), X2CJIB(*) 251 REAL WORK(LWORK) 252#else 253 DOUBLE PRECISION XLAMDP1(*), XLAMDH1(*) 254 DOUBLE PRECISION XLAMDP2(*), XLAMDH2(*) 255 DOUBLE PRECISION XLAMDP3(*), XLAMDH3(*) 256 DOUBLE PRECISION XLAMDP4(*), XLAMDH4(*) 257 DOUBLE PRECISION X1INT(*), X2INT(*) 258 DOUBLE PRECISION X1IJCB(*), X1CJIB(*) 259 DOUBLE PRECISION X2IJCB(*), X2CJIB(*) 260 DOUBLE PRECISION WORK(LWORK) 261#endif 262 263* local 264 265 INTEGER ISYL11, ISYL12, ISYL13, ISYL14, ISYL23 266 INTEGER ISYALP, ISYBET, ISYDEL, ISYGAM 267 INTEGER ISYALP1, ISYBET1 268 INTEGER ISYMI, ISYMJ, ISYMIJ, ISYMC 269 INTEGER KSCR2, KSCR4, KSCR41, KSCR5 270 INTEGER KSCR6, KSCR61, KSCR7, KEND1, LWRK1 271 INTEGER KOFF1, KOFF2, KOFF3, KOFF4, KOFF41, KLAMD 272 INTEGER KOFF5, KOFF6, KOFF61, KOFF7 273 INTEGER NBASA, NBASB, NVIRC, NRHFI 274 INTEGER ISYIJ1,ISYIJ2,ISYIJ3,ISYIJ4 275 INTEGER ISYCJ2,ISYCJ4,ISYCJ3 276 277*---------------------------------------------------------------------* 278* set some symmetries 279*---------------------------------------------------------------------* 280 281 ISYDEL = ISAO(IDEL) 282 ISYGAM = ISAO(IGAM) 283 284*---------------------------------------------------------------------* 285* work space allocation: 286* 287* KSCR2 -- I^{del,gam}(alp j); (alp j^); (alp j-bar); (alp j^-bar) 288* KSCR4 -- I^{del,gam}(i j) 289* KSCR41 -- I^{del,gam}(i j^) 290* 291* KSCR6 -- I^{del,gam}([i-bar j + i j-bar]) + I^[1]^{del,gam}(i j) 292* KSCR61 -- I^{del,gam}([i-bar j^ + i j^-bar]) + I^[1]^{del,gam}(i j^) 293* 294* KSCR5 -- I^{del,gam}([c^ j + c j^]) 295* KSCR7 -- I^{del,gam}([c^-bar j + c-bar j^ + c j^-bar + c^ j-bar]) 296* + I^[1]^{del,gam}(c^ j + c j^) 297* 298*---------------------------------------------------------------------* 299 KSCR2 = 1 300 KSCR4 = KSCR2 + NBAST*NRHFT 301 KSCR41 = KSCR4 + NRHFT*NRHFT 302 KEND1 = KSCR41 + NRHFT*NRHFT 303 304 IF (IOPT.EQ.1) THEN 305 KSCR5 = KEND1 306 KEND1 = KSCR5 + NVIRT*NRHFT 307 END IF 308 309 IF (LDERIV .OR. LRELAX) THEN 310 KSCR6 = KEND1 311 KSCR61 = KSCR6 + NRHFT*NRHFT 312 KEND1 = KSCR61 + NRHFT*NRHFT 313 314 IF (IOPT.EQ.1) THEN 315 KSCR7 = KEND1 316 KEND1 = KSCR7 + NVIRT*NRHFT 317 END IF 318 END IF 319 320 LWRK1 = LWORK - KEND1 321 322 IF ( LWRK1 .LT. 0) THEN 323 CALL QUIT('Insufficient memory in CC_IJCB.') 324 END IF 325 326*---------------------------------------------------------------------* 327* transform beta index to J using XLAMDH1 328* -- store (alf J|gam del) in SCR2 329*---------------------------------------------------------------------* 330c** 331c ISYSCR2 = MULD2H(ISY1ALBE,ISYML1) 332c CALL DZERO(NT1AO(ISYSCR2),WORK(KSCR2)) 333c** 334 KOFF2 = KSCR2 335 DO ISYMJ = 1, NSYM 336 337 ISYBET = MULD2H(ISYML1,ISYMJ) 338 ISYALP = MULD2H(ISYBET,ISY1ALBE) 339 340 KOFF1 = 1 + IAODIS(ISYALP,ISYBET) 341 KLAMD = IGLMRH(ISYBET,ISYMJ) + 1 342 343 NBASA = MAX(NBAS(ISYALP),1) 344 NBASB = MAX(NBAS(ISYBET),1) 345 346 CALL DGEMM('N','N',NBAS(ISYALP),NRHF(ISYMJ),NBAS(ISYBET), 347 * ONE,X1INT(KOFF1),NBASA,XLAMDH1(KLAMD), 348 * NBASB,ZERO,WORK(KOFF2),NBASA) 349 350 KOFF2 = KOFF2 + NBAS(ISYALP)*NRHF(ISYMJ) 351 352 END DO 353*---------------------------------------------------------------------* 354* transform alpha index to I using XLAMDP1 355* -- store (i j|gam del) in SCR4 356*---------------------------------------------------------------------* 357 KOFF2 = KSCR2 358 DO ISYMJ = 1, NSYM 359 360 ISYBET = MULD2H(ISYML1,ISYMJ) 361 ISYALP = MULD2H(ISYBET,ISY1ALBE) 362 ISYMI = MULD2H(ISYML1,ISYALP) 363 364 KLAMD = IGLMRH(ISYALP,ISYMI) + 1 365 KOFF4 = KSCR4 + IMATIJ(ISYMI,ISYMJ) 366 367 NBASA = MAX(NBAS(ISYALP),1) 368 NRHFI = MAX(NRHF(ISYMI),1) 369 370 CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),NBAS(ISYALP), 371 * ONE,XLAMDP1(KLAMD),NBASA,WORK(KOFF2), 372 * NBASA,ZERO,WORK(KOFF4),NRHFI) 373 374 KOFF2 = KOFF2 + NBAS(ISYALP)*NRHF(ISYMJ) 375 376 END DO 377 378*---------------------------------------------------------------------* 379* LRELAX transform alpha index to i-bar using XLAMDP2 380* -- store (i-bar j|gam del) in SCR6 381*---------------------------------------------------------------------* 382 IF ( LRELAX ) THEN 383 384 KOFF2 = KSCR2 385 386 DO ISYMJ = 1, NSYM 387 388 ISYBET = MULD2H(ISYML1,ISYMJ) 389 ISYALP = MULD2H(ISYBET,ISY1ALBE) 390 ISYMI = MULD2H(ISYML2,ISYALP) 391 392 KLAMD = IGLMRH(ISYALP,ISYMI) + 1 393 KOFF6 = KSCR6 + IMATIJ(ISYMI,ISYMJ) 394 395 NBASA = MAX(NBAS(ISYALP),1) 396 NRHFI = MAX(NRHF(ISYMI),1) 397 398 CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ), 399 * NBAS(ISYALP),ONE,XLAMDP2(KLAMD),NBASA, 400 * WORK(KOFF2),NBASA,ZERO,WORK(KOFF6),NRHFI) 401 402 KOFF2 = KOFF2 + NBAS(ISYALP)*NRHF(ISYMJ) 403 404 END DO 405 406 END IF 407 408*---------------------------------------------------------------------* 409* for IOPT=1 transform alpha index to C^ using XLAMDP3 410* -- store (c^ j|gam del) in SCR5 411*---------------------------------------------------------------------* 412 IF ( IOPT.EQ.1 ) THEN 413 414 KOFF2 = KSCR2 415 416 DO ISYMJ = 1, NSYM 417 418 ISYBET = MULD2H(ISYML1,ISYMJ) 419 ISYALP = MULD2H(ISYBET,ISY1ALBE) 420 ISYMC = MULD2H(ISYML3,ISYALP) 421 422 KLAMD = IGLMVI(ISYALP,ISYMC) + 1 423 KOFF5 = KSCR5 + IT1AM(ISYMC,ISYMJ) 424 425 NBASA = MAX(NBAS(ISYALP),1) 426 NVIRC = MAX(NVIR(ISYMC),1) 427 428 CALL DGEMM('T','N',NVIR(ISYMC),NRHF(ISYMJ),NBAS(ISYALP), 429 * ONE,XLAMDP3(KLAMD),NBASA,WORK(KOFF2),NBASA, 430 * ZERO,WORK(KOFF5),NVIRC) 431 432 KOFF2 = KOFF2 + NBAS(ISYALP)*NRHF(ISYMJ) 433c 434 END DO 435 436 END IF 437*---------------------------------------------------------------------* 438* for IOPT=1 and LRELAX transform alpha index to C^(bar) 439* using XLAMDP4 -- store (c^-bar j|gam del) in SCR7 440*---------------------------------------------------------------------* 441 IF (( IOPT.EQ.1 ).AND.LRELAX) THEN 442 443 KOFF2 = KSCR2 444 445 DO ISYMJ = 1, NSYM 446 447 ISYBET = MULD2H(ISYML1,ISYMJ) 448 ISYALP = MULD2H(ISYBET,ISY1ALBE) 449 ISYMC = MULD2H(ISYML4,ISYALP) 450 451 KLAMD = IGLMVI(ISYALP,ISYMC) + 1 452 KOFF7 = KSCR7 + IT1AM(ISYMC,ISYMJ) 453 454 NBASA = MAX(NBAS(ISYALP),1) 455 NVIRC = MAX(NVIR(ISYMC),1) 456 457 CALL DGEMM('T','N',NVIR(ISYMC),NRHF(ISYMJ),NBAS(ISYALP), 458 * ONE,XLAMDP4(KLAMD),NBASA,WORK(KOFF2),NBASA, 459 * ZERO,WORK(KOFF7),NVIRC) 460 461 KOFF2 = KOFF2 + NBAS(ISYALP)*NRHF(ISYMJ) 462 463 END DO 464 465 END IF 466 467*---------------------------------------------------------------------* 468* transform beta index to J^ using XLAMDH3 469* -- store (alf J^|gam del) in SCR2 470*---------------------------------------------------------------------* 471 KOFF2 = KSCR2 472 473 DO ISYMJ = 1, NSYM 474 475 ISYBET = MULD2H(ISYML3,ISYMJ) 476 ISYALP = MULD2H(ISYBET,ISY1ALBE) 477 478 KOFF1 = 1 + IAODIS(ISYALP,ISYBET) 479 KLAMD = IGLMRH(ISYBET,ISYMJ) + 1 480 481 NBASA = MAX(NBAS(ISYALP),1) 482 NBASB = MAX(NBAS(ISYBET),1) 483 484 CALL DGEMM('N','N',NBAS(ISYALP),NRHF(ISYMJ),NBAS(ISYBET), 485 * ONE,X1INT(KOFF1),NBASA,XLAMDH3(KLAMD), 486 * NBASB,ZERO,WORK(KOFF2),NBASA) 487 488 KOFF2 = KOFF2 + NBAS(ISYALP)*NRHF(ISYMJ) 489 490 END DO 491 492*---------------------------------------------------------------------* 493* transform alpha index to I using XLAMDP1 494* -- store (i j^|gam del) in SCR41 495*---------------------------------------------------------------------* 496 KOFF2 = KSCR2 497 498 DO ISYMJ = 1, NSYM 499 500 ISYBET = MULD2H(ISYML3,ISYMJ) 501 ISYALP = MULD2H(ISYBET,ISY1ALBE) 502 ISYMI = MULD2H(ISYML1,ISYALP) 503 504 KLAMD = IGLMRH(ISYALP,ISYMI) + 1 505 KOFF41 = KSCR41 + IMATIJ(ISYMI,ISYMJ) 506 507 NBASA = MAX(NBAS(ISYALP),1) 508 NRHFI = MAX(NRHF(ISYMI),1) 509 510 CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),NBAS(ISYALP), 511 * ONE,XLAMDP1(KLAMD),NBASA,WORK(KOFF2), 512 * NBASA,ZERO,WORK(KOFF41),NRHFI) 513 514 KOFF2 = KOFF2 + NBAS(ISYALP)*NRHF(ISYMJ) 515 516 END DO 517 518*---------------------------------------------------------------------* 519* if LRELAX transform alpha index to i-bar using XLAMDP2 520* -- store (i-bar j^| gam del) in SCR61 521*---------------------------------------------------------------------* 522 IF ( LRELAX ) THEN 523 524 KOFF2 = KSCR2 525 526 DO ISYMJ = 1, NSYM 527 528 ISYBET = MULD2H(ISYML3,ISYMJ) 529 ISYALP = MULD2H(ISYBET,ISY1ALBE) 530 ISYMI = MULD2H(ISYML2,ISYALP) 531 532 KLAMD = IGLMRH(ISYALP,ISYMI) + 1 533 KOFF61 = KSCR61 + IMATIJ(ISYMI,ISYMJ) 534 535 NBASA = MAX(NBAS(ISYALP),1) 536 NRHFI = MAX(NRHF(ISYMI),1) 537 538 CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),NBAS(ISYALP), 539 * ONE,XLAMDP2(KLAMD),NBASA,WORK(KOFF2),NBASA, 540 * ZERO,WORK(KOFF61),NRHFI) 541 542 KOFF2 = KOFF2 + NBAS(ISYALP)*NRHF(ISYMJ) 543 544 END DO 545 546 END IF 547 548*---------------------------------------------------------------------* 549* for IOPT=1 transform alpha index to C using XLAMDP1 550* -- add (c j^|gam del) to (c^ j|gam del) in SCR5 551*---------------------------------------------------------------------* 552 IF ( IOPT.EQ.1 ) THEN 553 554 KOFF2 = KSCR2 555 556 DO ISYMJ = 1, NSYM 557 558 ISYBET = MULD2H(ISYML3,ISYMJ) 559 ISYALP = MULD2H(ISYBET,ISY1ALBE) 560 ISYMC = MULD2H(ISYML1,ISYALP) 561 562 KLAMD = IGLMVI(ISYALP,ISYMC) + 1 563 KOFF5 = KSCR5 + IT1AM(ISYMC,ISYMJ) 564 565 NBASA = MAX(NBAS(ISYALP),1) 566 NVIRC = MAX(NVIR(ISYMC),1) 567 568 CALL DGEMM('T','N',NVIR(ISYMC),NRHF(ISYMJ),NBAS(ISYALP), 569 * ONE,XLAMDP1(KLAMD),NBASA,WORK(KOFF2),NBASA, 570 * ONE,WORK(KOFF5),NVIRC) 571 572 KOFF2 = KOFF2 + NBAS(ISYALP)*NRHF(ISYMJ) 573 574 END DO 575 END IF 576*---------------------------------------------------------------------* 577* for IOPT=1.and.LRELAX transform alpha to C-bar using XLAMDP2 578* -- add (c-bar j^|gam del) to (c^-bar j|gam del) in SCR7 579*---------------------------------------------------------------------* 580 IF (( IOPT.EQ.1 ).AND.(LRELAX)) THEN 581 582 KOFF2 = KSCR2 583 584 DO ISYMJ = 1, NSYM 585 586 ISYBET = MULD2H(ISYML3,ISYMJ) 587 ISYALP = MULD2H(ISYBET,ISY1ALBE) 588 ISYMC = MULD2H(ISYML2,ISYALP) 589 590 KLAMD = IGLMVI(ISYALP,ISYMC) + 1 591 KOFF7 = KSCR7 + IT1AM(ISYMC,ISYMJ) 592 593 NBASA = MAX(NBAS(ISYALP),1) 594 NVIRC = MAX(NVIR(ISYMC),1) 595 596 CALL DGEMM('T','N',NVIR(ISYMC),NRHF(ISYMJ),NBAS(ISYALP), 597 * ONE,XLAMDP2(KLAMD),NBASA,WORK(KOFF2),NBASA, 598 * ONE,WORK(KOFF7),NVIRC) 599 600 KOFF2 = KOFF2 + NBAS(ISYALP)*NRHF(ISYMJ) 601 602 END DO 603 604 END IF 605C 606C Finished with SCR2 again 607C 608*---------------------------------------------------------------------* 609* for LRELAX add extra contributions from (alp j-bar|gam del), 610* (alp ^j-bar|gam del) etc 611*---------------------------------------------------------------------* 612 IF ( LRELAX ) THEN 613*---------------------------------------------------------------------* 614* transform beta to J-bar using XLAMDH2 615* -- store (alp j-bar|gam del) in SCR2 616* If (LDERIV) add also (alp j|gam del)[1] 617*---------------------------------------------------------------------* 618 619 KOFF2 = KSCR2 620 621 DO ISYMJ = 1, NSYM 622 623 ISYBET = MULD2H(ISYML2,ISYMJ) 624 ISYALP = MULD2H(ISYBET,ISY1ALBE) 625 626 KOFF1 = 1 + IAODIS(ISYALP,ISYBET) 627 KLAMD = IGLMRH(ISYBET,ISYMJ) + 1 628 629 NBASA = MAX(NBAS(ISYALP),1) 630 NBASB = MAX(NBAS(ISYBET),1) 631 632 CALL DGEMM('N','N',NBAS(ISYALP),NRHF(ISYMJ),NBAS(ISYBET), 633 * ONE,X1INT(KOFF1),NBASA,XLAMDH2(KLAMD), 634 * NBASB,ZERO,WORK(KOFF2),NBASA) 635 636 IF (LDERIV) THEN 637 638 ISYBET1 = MULD2H(ISYML1,ISYMJ) 639 ISYALP1 = MULD2H(ISYBET1,ISY2ALBE) 640 IF (ISYALP1.NE.ISYALP) 641 * CALL QUIT('Symmetry mismatch in CC_IJCB') 642 643 KOFF3 = 1 + IAODIS(ISYALP1,ISYBET1) 644 KLAMD = IGLMRH(ISYBET1,ISYMJ) + 1 645 NBASA = MAX(NBAS(ISYALP1),1) 646 NBASB = MAX(NBAS(ISYBET1),1) 647 CALL DGEMM('N','N',NBAS(ISYALP1),NRHF(ISYMJ), 648 * NBAS(ISYBET1),ONE,X2INT(KOFF3),NBASA, 649 * XLAMDH1(KLAMD),NBASB,ONE,WORK(KOFF2),NBASA) 650 END IF 651 652 653 KOFF2 = KOFF2 + NBAS(ISYALP)*NRHF(ISYMJ) 654 655 END DO 656 657*---------------------------------------------------------------------* 658* transform alpha to I using XLAMDP1 659* -- add (i j-bar|gam del) to (i-bar j|gam del) in SCR6 660*---------------------------------------------------------------------* 661 662 KOFF2 = KSCR2 663 664 DO ISYMJ = 1, NSYM 665 666 ISYBET = MULD2H(ISYML2,ISYMJ) 667 ISYALP = MULD2H(ISYBET,ISY1ALBE) 668 ISYMI = MULD2H(ISYML1,ISYALP) 669 670 KLAMD = IGLMRH(ISYALP,ISYMI) + 1 671 KOFF6 = KSCR6 + IMATIJ(ISYMI,ISYMJ) 672 673 NBASA = MAX(NBAS(ISYALP),1) 674 NRHFI = MAX(NRHF(ISYMI),1) 675 676 CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ), 677 * NBAS(ISYALP),ONE,XLAMDP1(KLAMD),NBASA, 678 * WORK(KOFF2), NBASA,ONE,WORK(KOFF6),NRHFI) 679 680 KOFF2 = KOFF2 + NBAS(ISYALP)*NRHF(ISYMJ) 681 682 END DO 683*---------------------------------------------------------------------* 684* for IOPT=1 transform alpha to C^ using XLAMDP3 685* -- add (c^ j-bar | gam del) in SCR7 686*---------------------------------------------------------------------* 687 IF ( IOPT.EQ.1 ) THEN 688 689 KOFF2 = KSCR2 690 691 DO ISYMJ = 1, NSYM 692 693 ISYBET = MULD2H(ISYML2,ISYMJ) 694 ISYALP = MULD2H(ISYBET,ISY1ALBE) 695 ISYMC = MULD2H(ISYML3,ISYALP) 696 697 KLAMD = IGLMVI(ISYALP,ISYMC) + 1 698 KOFF7 = KSCR7 + IT1AM(ISYMC,ISYMJ) 699 700 NBASA = MAX(NBAS(ISYALP),1) 701 NVIRC = MAX(NVIR(ISYMC),1) 702 703 CALL DGEMM('T','N',NVIR(ISYMC),NRHF(ISYMJ), 704 * NBAS(ISYALP),ONE,XLAMDP3(KLAMD),NBASA, 705 * WORK(KOFF2),NBASA,ONE,WORK(KOFF7),NVIRC) 706 707 KOFF2 = KOFF2 + NBAS(ISYALP)*NRHF(ISYMJ) 708 709 END DO 710 END IF 711 712*---------------------------------------------------------------------* 713* transform beta to J^-bar using XLAMDH4 714* -- store (alp j^-bar|gam del) in SCR2 715* if (LDERIV) add derivative contribution (al j^|gam del)(1) 716*---------------------------------------------------------------------* 717 718 KOFF2 = KSCR2 719 720 DO ISYMJ = 1, NSYM 721 722 ISYBET = MULD2H(ISYML4,ISYMJ) 723 ISYALP = MULD2H(ISYBET,ISY1ALBE) 724 725 KOFF1 = 1 + IAODIS(ISYALP,ISYBET) 726 KLAMD = IGLMRH(ISYBET,ISYMJ) + 1 727 728 NBASA = MAX(NBAS(ISYALP),1) 729 NBASB = MAX(NBAS(ISYBET),1) 730 731 CALL DGEMM('N','N',NBAS(ISYALP),NRHF(ISYMJ),NBAS(ISYBET), 732 * ONE,X1INT(KOFF1),NBASA,XLAMDH4(KLAMD), 733 * NBASB,ZERO,WORK(KOFF2),NBASA) 734 735 IF (LDERIV) THEN 736 737 ISYBET1 = MULD2H(ISYML3,ISYMJ) 738 ISYALP1 = MULD2H(ISYBET1,ISY2ALBE) 739 IF (ISYALP1.NE.ISYALP) 740 * CALL QUIT('Symmetry mismatch in CC_IJCB') 741 742 KOFF3 = 1 + IAODIS(ISYALP1,ISYBET1) 743 KLAMD = IGLMRH(ISYBET1,ISYMJ) + 1 744 NBASA = MAX(NBAS(ISYALP1),1) 745 NBASB = MAX(NBAS(ISYBET1),1) 746 CALL DGEMM('N','N',NBAS(ISYALP1),NRHF(ISYMJ), 747 * NBAS(ISYBET1),ONE,X2INT(KOFF3),NBASA, 748 * XLAMDH3(KLAMD),NBASB,ONE,WORK(KOFF2),NBASA) 749 END IF 750 751 KOFF2 = KOFF2 + NBAS(ISYALP)*NRHF(ISYMJ) 752 753 END DO 754*---------------------------------------------------------------------* 755* transform alpha to I using XLAMDP1 756* -- add (i j^-bar|gam del) to (i-bar j^|gam del) in SCR61 757*---------------------------------------------------------------------* 758 759 KOFF2 = KSCR2 760 761 DO ISYMJ = 1, NSYM 762 763 ISYBET = MULD2H(ISYML4,ISYMJ) 764 ISYALP = MULD2H(ISYBET,ISY1ALBE) 765 ISYMI = MULD2H(ISYML1,ISYALP) 766 767 KLAMD = IGLMRH(ISYALP,ISYMI) + 1 768 KOFF61 = KSCR61 + IMATIJ(ISYMI,ISYMJ) 769 770 NBASA = MAX(NBAS(ISYALP),1) 771 NRHFI = MAX(NRHF(ISYMI),1) 772 773 CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ), 774 * NBAS(ISYALP),ONE,XLAMDP1(KLAMD),NBASA, 775 * WORK(KOFF2),NBASA,ONE,WORK(KOFF61),NRHFI) 776 777 KOFF2 = KOFF2 + NBAS(ISYALP)*NRHF(ISYMJ) 778 779 END DO 780*---------------------------------------------------------------------* 781* for IOPT=1 transform alpha to C using XLAMDP1 782* -- add (c j^-bar | gam del) in SCR7 783*---------------------------------------------------------------------* 784 IF ( IOPT.EQ.1 ) THEN 785 786 KOFF2 = KSCR2 787 788 DO ISYMJ = 1, NSYM 789 790 ISYBET = MULD2H(ISYML4,ISYMJ) 791 ISYALP = MULD2H(ISYBET,ISY1ALBE) 792 ISYMC = MULD2H(ISYML1,ISYALP) 793 794 KLAMD = IGLMVI(ISYALP,ISYMC) + 1 795 KOFF7 = KSCR7 + IT1AM(ISYMC,ISYMJ) 796 797 NBASA = MAX(NBAS(ISYALP),1) 798 NVIRC = MAX(NVIR(ISYMC),1) 799 800 CALL DGEMM('T','N',NVIR(ISYMC),NRHF(ISYMJ), 801 * NBAS(ISYALP),ONE,XLAMDP1(KLAMD),NBASA, 802 * WORK(KOFF2),NBASA,ONE,WORK(KOFF7),NVIRC) 803 804 KOFF2 = KOFF2 + NBAS(ISYALP)*NRHF(ISYMJ) 805 806 END DO 807 END IF 808 END IF 809*---------------------------------------------------------------------* 810* Add the contribution to the result X1IJCB and X2IJCB vector: 811* (transform the gamma index to VIRTUAL) 812* Note that IJCB integrals are sorted as I_cj,i,del 813*---------------------------------------------------------------------* 814 ISYL11 = MULD2H(ISYML1,ISYML1) 815 ISYL12 = MULD2H(ISYML1,ISYML2) 816 ISYL13 = MULD2H(ISYML1,ISYML3) 817 ISYL14 = MULD2H(ISYML1,ISYML4) 818 ISYL23 = MULD2H(ISYML2,ISYML3) 819 820 IF ( LNEWTA ) THEN 821C ------------------------- 822C add (ij|c^ del) to X1IJCB: 823C ------------------------- 824 ISYIJ1 = MULD2H(ISY1ALBE,ISYL11) 825 826 CALL CC_IJCB2(IGAM, WORK(KSCR4), ISYIJ1, ISYGAM, 827 & XLAMDP3, ISYML3, X1IJCB) 828C ------------------------- 829C add (ij^|c del) to X1IJCB: 830C ------------------------- 831 ISYIJ3 = MULD2H(ISY1ALBE,ISYL13) 832 CALL CC_IJCB2(IGAM, WORK(KSCR41), ISYIJ3, ISYGAM, 833 & XLAMDP1, ISYML1, X1IJCB) 834 END IF 835 836 837 IF ( LRELAX ) THEN 838C ------------------------------ 839C add (ij|c^-bar del) to X2IJCB: 840C ------------------------------ 841 ISYIJ1 = MULD2H(ISY1ALBE,ISYL11) 842 843 CALL CC_IJCB2(IGAM, WORK(KSCR4), ISYIJ1, ISYGAM, 844 & XLAMDP4, ISYML4, X2IJCB) 845C ------------------------------ 846C add (ij^|c-bar del) to X2IJCB: 847C ------------------------------ 848 ISYIJ3 = MULD2H(ISY1ALBE,ISYL13) 849 850 CALL CC_IJCB2(IGAM, WORK(KSCR41), ISYIJ3, ISYGAM, 851 & XLAMDP2, ISYML2, X2IJCB) 852 END IF 853 IF ( LDERIV .OR. LRELAX ) THEN 854C ------------------------------------------------ 855C add ([i-bar j^ + i j^-bar]|c del) to X2IJCB: 856C ------------------------------------------------ 857 ISYIJ4 = MULD2H(ISY1ALBE,ISYL14) 858 859 CALL CC_IJCB2(IGAM, WORK(KSCR61), ISYIJ4, ISYGAM, 860 & XLAMDP1, ISYML1, X2IJCB) 861C ------------------------------------------------ 862C add ([i-bar j + i j-bar]|c^ del) to X2IJCB: 863C ------------------------------------------------ 864 ISYIJ2 = MULD2H(ISY1ALBE,ISYL12) 865 866 CALL CC_IJCB2(IGAM, WORK(KSCR6), ISYIJ2, ISYGAM, 867 & XLAMDP3, ISYML3, X2IJCB) 868 END IF 869 870*---------------------------------------------------------------------* 871* Add the contribution to the result X1CJIB and X2CJIB vector: 872* transform gamma to OCCUPIED 873* Integrals sorted as I_cj,i,del 874*---------------------------------------------------------------------* 875 IF ( IOPT.EQ.1 ) THEN 876 877 IF ( LNEWTA ) THEN 878C ------------------------- 879C add ([c^ j + c j^]|i del) to X1CJIB: 880C ------------------------- 881c*** 882 ISYCJ3 = MULD2H(ISY1ALBE,ISYL13) 883 884 CALL CC_IAJB1(IGAM, WORK(KSCR5), ISYCJ3, ISYGAM, 885 & XLAMDP1, ISYML1, X1CJIB, .FALSE., IDUMMY) 886 END IF 887 888 889 IF ( LRELAX ) THEN 890C ------------------------------ 891C add ([c^ j + c j^]|i-bar del) to X2CJIB: 892C ------------------------------ 893 ISYCJ3 = MULD2H(ISY1ALBE,ISYL13) 894 895 CALL CC_IAJB1(IGAM, WORK(KSCR5), ISYCJ3, ISYGAM, 896 & XLAMDP2, ISYML2, X2CJIB, .FALSE., IDUMMY) 897 END IF 898 899 IF ( LDERIV .OR. LRELAX ) THEN 900C ------------------------------------------------ 901C add ([c-bar j^ + c j^-bar + c^-bar j + c^ j-bar] | i del) to X2IABJ: 902C ------------------------------------------------ 903 ISYCJ4 = MULD2H(ISY1ALBE,ISYL14) 904 905 CALL CC_IAJB1(IGAM, WORK(KSCR7), ISYCJ4, ISYGAM, 906 & XLAMDP1, ISYML1, X2CJIB, .FALSE., IDUMMY) 907 END IF 908 909 END IF 910 911 RETURN 912 END 913*=====================================================================* 914* END OF SUBROUTINE CCIJCB * 915*=====================================================================* 916c /* deck cc_ijcb2 */ 917*=====================================================================* 918 SUBROUTINE CC_IJCB2(IGAM, XIJG, ISYMIJ, ISYGAM, 919 & XLAMDA, ISYLAM, XIJCB ) 920*---------------------------------------------------------------------* 921* 922* Purpose: transform (ij|gam del) to (ij|c del), with sorting 923* I_cj,i,del 924* Sonia, March 1999 925*---------------------------------------------------------------------* 926#if defined (IMPLICIT_NONE) 927 IMPLICIT NONE 928#else 929# include "implicit.h" 930#endif 931#include "ccsdsym.h" 932#include "ccorb.h" 933#include "maxorb.h" 934#include "ccisao.h" 935 936#if defined (SYS_CRAY) 937 REAL XIJG(*), XIJCB(*), XLAMDA(*) 938#else 939 DOUBLE PRECISION XIJG(*), XIJCB(*), XLAMDA(*) 940#endif 941 942 INTEGER IGAM, ISYMIJ, ISYGAM, ISYLAM, ISYMC, KLAMD, KRES 943 INTEGER ISYMCJ, ISYMI, ISYMJ, KOFF5, NBASG 944 INTEGER KXIJ 945 946* transform integral batch: 947 ISYMC = MULD2H(ISYGAM,ISYLAM) 948 G = IGAM - IBAS(ISYGAM) 949 NBASG = MAX(NBAS(ISYGAM),1) 950 951 DO ISYMJ = 1, NSYM 952 ISYMI = MULD2H(ISYMIJ,ISYMJ) 953 ISYMCJ = MULD2H(ISYMC,ISYMJ) 954 955 DO J = 1, NRHF(ISYMJ) 956 DO I = 1, NRHF(ISYMI) 957 958 KLAMD = IGLMVI(ISYGAM,ISYMC) + G !offs Lambda_gam,c 959 !offset I_ij,gam 960 KXIJ = IMATIJ(ISYMI,ISYMJ) + NRHF(ISYMI)*(J-1) + I 961 !offset I_cj,i 962 KRES = IT2BCD(ISYMCJ,ISYMI) + NT1AM(ISYMCJ)*(I-1) 963 & + IT1AM(ISYMC,ISYMJ) + NVIR(ISYMC)*(J-1) + 1 964 965 CALL DAXPY(NVIR(ISYMC),XIJG(KXIJ),XLAMDA(KLAMD), 966 & NBASG,XIJCB(KRES),1) 967 END DO 968 END DO 969 970 END DO 971 972 RETURN 973 END 974*=====================================================================* 975* END OF SUBROUTINE CC_IJCB2 * 976*=====================================================================* 977