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 19*=====================================================================* 20C /* Deck cc_bfden */ 21 SUBROUTINE CC_BFDEN( T2AM, ISYTAM, MINT, ISYMIN, 22 * XLAMD1, ISYML1, XLAMD2, ISYML2, 23 * XLAMD3, ISYML3, XLAMD4, ISYML4, 24 * FILNAM, LUBFDEN, IADRBF, IADR, 25 * IVEC, IOPT, LO3BF, WORK, LWORK ) 26*---------------------------------------------------------------------* 27* 28* Purpose: Calculate the effective density used for the BF terms 29* write it to the direct access file FILNAM/LUBFDEN 30* start addresses are stored in IADRBF 31* 32* 33* IOPT .EQ. 0: two-index back transformation of T --> evaluate: 34* 35* M^{g,d}_{i,j} = T^{g,d}_{i,j}(1,2) 36* = sum_ab T^{a,b}_{i,j} XLAMD1(d,b) XLAMD2(g,a) 37* 38* XLAMD1,XLAMD2 : (ordinary) lambda hole matrices 39* XLAMD3,XLAMD4,MINT unused 40* 41* 42* IOPT .EQ. 1: vector function --> evaluate: 43* 44* M^{g,d}_{i,j} = T^{g,d}_{i,j}(1,2) 45* + XLAMD1(g,i)*XLAMD3(d,j) 46* 47* XLAMD1,XLAMD2,XLAMD3 : (ordinary) lambda hole matrices 48* XLAMD4,MINT unused 49* 50* 51* IOPT .EQ. 2: Jacobian left transf. --> evaluate: 52* 53* M^{g,d}_{i,j} = T^{g,d}_{i,j}(1,2) + Mint^{g,d}_{i,j}(1,2) 54* + 2 XLAMD1(g,i)*XLAMD3(d,j) - XLAMD3(d,i)*XLAMD1(g,j) 55* + 2 XLAMD3(g,i)*XLAMD1(d,j) - XLAMD1(d,i)*XLAMD3(g,j) 56* 57* XLAMD1: zeroth-order lambda particle matrix 58* XLAMD2: zeroth-order lambda particle matrix 59* XLAMD3: zeroth-order chi intermediate 60* XLAMD4: unused 61* MINT : M intermediate in left transf. 62* 63* 64* 65* IOPT .EQ. 3: Jacobian right transformation and 66* B matrix transformation --> evaluate: 67* !!! not yet tested for this version !!! 68* 69* 70* M^{g,d}_{i,j} = T^{g,d}_{i,j}(1,2) 71* + XLAMD1(gam,i)*XLAMD3(del,j) 72* + XLAMD3(gam,i)*XLAMD1(del,j) 73* 74* XLAMD1,XLAMD2,XLAMD3 : lambda hole matrices 75* XLAMD4,MINT unused 76* 77* 78* 79* 80* IOPT .EQ. 4: F-matrix transformation --> evaluate 81* !!! not yet tested for this version !!! 82* 83* M^{g,d}_{i,j} = T^{g,d}_{i,j}(1,2) + T^{g,d}_{i,j}(2,1) 84* + 2 XLAMD1(g,i)*XLAMD3(d,j) - XLAMD3(d,i)*XLAMD1(g,j) 85* + 2 XLAMD3(g,i)*XLAMD1(d,j) - XLAMD1(d,i)*XLAMD3(g,j) 86* 87* XLAMD1: zeroth-order lambda particle matrix 88* XLAMD2: first-order lambda particle matrix 89* XLAMD3: first-order one-index backtransf. Zeta1 amplitudes 90* XLAMD4: unused 91* 92* 93* 94* IOPT .EQ. 5: relaxation contrib. to Jacobian right trans. --> 95* 96* M^{g,d}_{i,j} = T^{g,d}_{i,j}(1,2) + T^{g,d}_{i,j}(2,1) 97* + XLAMD1(g,i)*XLAMD3(d,j) 98* + XLAMD2(g,i)*XLAMD4(d,j) 99* 100* XLAMD1: first-order lambda hole matrix 101* XLAMD2: zeroth-order lambda hole matrix 102* XLAMD3: zeroth-order lambda hole matrix 103* XLAMD4: first-order lambda hole matrix 104* 105* 106* IOPT .EQ. 6: relaxation contrib. to Jacobian left trans. --> 107* 108* M^{g,d}_{i,j} = T^{g,d}_{i,j}(1,2) + T^{g,d}_{i,j}(2,1) 109* + 2 XLAMD1(g,i)*XLAMD3(d,j) - XLAMD3(d,i)*XLAMD1(g,j) 110* + 2 XLAMD3(g,i)*XLAMD1(d,j) - XLAMD1(d,i)*XLAMD3(g,j) 111* + 2 XLAMD2(g,i)*XLAMD4(d,j) - XLAMD4(d,i)*XLAMD2(g,j) 112* + 2 XLAMD4(g,i)*XLAMD2(d,j) - XLAMD2(d,i)*XLAMD4(g,j) 113* 114* XLAMD1: first-order lambda hole matrix 115* XLAMD2: zeroth-order lambda hole matrix 116* XLAMD3: zeroth-order chi intermediate 117* XLAMD4: first-order chi intermediate 118* 119* IOPT .EQ. 7: relax. + resp. contribs to eff. dens. for BF(A,QA) 120* 121* M^{g,d}_{i,j} = T^{g,d}_{i,j}(1,2) + T^{g,d}_{i,j}(2,1) 122* + XLAMD1(g,i)*XLAMD3(d,j) + XLAMD3(g,i)*XLAMD1(d,j) 123* + XLAMD2(g,i)*XLAMD4(d,j) + XLAMD4(g,i)*XLAMD2(d,j) 124* 125* XLAMD1: first-order lambda hole matrix 126* XLAMD2: zeroth-order lambda hole matrix 127* XLAMD3: T^A tranformed zero-order lambda hole matrix 128* XLAMD4: T^A tranformed first-order lambda hole matrix 129* 130* 131* IOPT .EQ. 8 used for F matrix transformation 132* 133* IOPT .EQ. 9: symmetrized two-index back transformation of T 134* 135* M^{g,d}_{i,j} = T^{g,d}_{i,j}(1,2) + T^{g,d}_{i,j}(2,1) 136* 137* XLAMD1,XLAMD2 : lambda hole matrices 138* XLAMD3,XLAMD4,MINT unused 139* 140* 141* 142* LO3BF - flag, for special treatment in F matrix 143* (unused at present?) 144* 145* SYMMETRIES: 146* 147* ISYML1,ISYML2,ISYML3,ISYML4 -- XLAMD1, XLAMD2, XLAMD3, XLAMD4 148* ISYTAM -- T2AM 149* ISYMIN -- MINT 150* 151* Christof Haettig, September 1998 152* IOPT = 7 added, Sonia Coriani, June 1999 153* IOPT = 0 added, Christof Haettig, Februar 2000 154* 155*=====================================================================* 156#if defined (IMPLICIT_NONE) 157 IMPLICIT NONE 158#else 159# include "implicit.h" 160#endif 161#include "priunit.h" 162#include "ccinftap.h" 163#include "ccorb.h" 164#include "maxorb.h" 165#include "ccsdsym.h" 166#include "second.h" 167C 168 LOGICAL LO3BF 169 CHARACTER*(*) FILNAM 170 INTEGER LUBFDEN, IVEC, IADR, IADRBF(MXCORB_CC,*) 171 INTEGER ISYML1,ISYML2,ISYML3,ISYML4,ISYTAM,ISYMIN,IOPT,LWORK 172 INTEGER KEND1,ISYMT1,ISYMT2,ISYDEN,ISYMD,ISYMM1,ISYMM2,NSCRM,NMGD 173 INTEGER KSCRM,KMGD,LWRK1,ISYMGD,IOPT2AO,IDEL,IOPBFAO 174C 175#if defined (SYS_CRAY) 176 REAL XLAMD1(*),XLAMD2(*),XLAMD3(*),XLAMD4(*) 177 REAL T2AM(*), MINT(*), WORK(*) 178 REAL TIMIO, DTIME, ZERO, ONE, DDOT, XNORM 179#else 180 DOUBLE PRECISION XLAMD1(*),XLAMD2(*),XLAMD3(*),XLAMD4(*) 181 DOUBLE PRECISION T2AM(*), MINT(*), WORK(*) 182 DOUBLE PRECISION TIMIO, DTIME, ZERO, ONE, DDOT, XNORM 183#endif 184C 185 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 186C 187C----------------------- 188C symmetry analysis: 189C----------------------- 190C 191 ISYMT1 = MULD2H(ISYTAM,ISYML1) 192 ISYMT2 = MULD2H(ISYTAM,ISYML2) 193 ISYDEN = MULD2H(ISYMT1,ISYML2) 194C 195 TIMIO = ZERO 196C 197 IF ( IOPT.NE.0 .AND. IOPT.NE.9) THEN 198 IF ( ISYML3.NE.ISYMT2 ) THEN 199 WRITE (LUPRI,*) 'ISYML3:',ISYML3 200 WRITE (LUPRI,*) 'ISYMT2:',ISYMT2 201 CALL QUIT('CC_BFDEN: Symmetry mismatch') 202 EN DIF 203 END IF 204C 205 IF ( (IOPT.EQ.5) .OR. (IOPT.EQ.6) .OR. (IOPT.EQ.7)) THEN 206 IF ( MULD2H(ISYML4,ISYML2) .NE. ISYDEN) THEN 207 WRITE (LUPRI,*) 'IOPT:',IOPT 208 WRITE (LUPRI,*) 'ISYML4:',ISYML4 209 WRITE (LUPRI,*) 'ISYML2:',ISYML2 210 WRITE (LUPRI,*) 'ISYDEN:',ISYDEN 211 CALL QUIT('CC_BFDEN: Symmetry mismatch') 212 ENDIF 213 END IF 214C 215C--------------------------------- 216C Start loop over Delta index: 217C--------------------------------- 218C 219 DO ISYMD = 1, NSYM 220 221 ISYMM1 = MULD2H(ISYMT1,ISYMD) 222 ISYMM2 = MULD2H(ISYMT2,ISYMD) 223 ISYMGD = MULD2H(ISYDEN,ISYMD) 224 225 NSCRM = MAX(NT2BCD(ISYMM1),NT2BCD(ISYMM2)) 226 NMGD = NT2BGD(ISYMGD) 227 IF (LO3BF .AND. IOPT.EQ.4) NMGD = NMAIJK(ISYMGD) 228 229 KSCRM = 1 230 KMGD = KSCRM + NSCRM 231 KEND1 = KMGD + NMGD 232 LWRK1 = LWORK - KEND1 233 234 IF (LWRK1 .LT. 0) THEN 235 WRITE (LUPRI,*) 'Need : ',KEND1,'Available : ',LWORK 236 WRITE (LUPRI,*) 'Insufficient space in CC_BFDEN.' 237 CALL QUIT('Insufficient space in CC_BFDEN.') 238 ENDIF 239C 240 DO D = 1, NBAS(ISYMD) 241 242 IDEL = D + IBAS(ISYMD) 243C 244C ----------------------------------------------------------- 245C calculate one-index AO back-transformed T2 amplitudes 246C back-transforme Delta index with XLAMD1, 247C ----------------------------------------------------------- 248C 249 IOPT2AO = 1 250 CALL CC_T2AO(T2AM,XLAMD1,ISYML1,WORK(KSCRM),WORK(KEND1), 251 * LWRK1,IDEL,ISYMD,ISYTAM,IOPT2AO) 252C 253C ----------------------------------------------------------- 254C back-transforme Gamma index with XLAMD2, 255C add "F-term" from XLAMD1 x XLAMD3 (MGD is initialized here) 256C ----------------------------------------------------------- 257C 258 IF (IOPT .EQ. 7) THEN 259 IOPBFAO = 3 260 ELSE IF (IOPT.EQ. 9) THEN 261 IOPBFAO = 0 262 ELSE 263 IOPBFAO = IOPT 264 ENDIF 265 266 CALL CC_BFAO(WORK(KMGD),WORK(KSCRM),ISYMM1,XLAMD2,ISYML2, 267 * D,ISYMD,XLAMD3,ISYML3,XLAMD1,ISYML1,ZERO,IOPBFAO) 268C 269C ----------------------------------------------------------- 270C for IOPT=4,5,6,7,9 add extra response (or relax.) contribs. 271C back-transforme Delta index with XLAMD2, 272C back-transforme Gamma index with XLAMD1, 273C add "F-term" from XLAMD2 x XLAMD4 274C ----------------------------------------------------------- 275C 276 IF (IOPT.EQ.4 .OR. IOPT.EQ.5 .OR. IOPT.EQ.6 .OR. 277 * IOPT.EQ.7 .OR. IOPT.EQ.9 ) THEN 278 279 IOPT2AO = 1 280 CALL CC_T2AO(T2AM,XLAMD2,ISYML2,WORK(KSCRM),WORK(KEND1), 281 * LWRK1,IDEL,ISYMD,ISYTAM,IOPT2AO) 282 283 IF (IOPT .EQ. 7) THEN 284 IOPBFAO = 3 285 ELSE IF (IOPT.EQ. 9) THEN 286 IOPBFAO = 0 287 ELSE 288 IOPBFAO = IOPT 289 ENDIF 290 291 CALL CC_BFAO(WORK(KMGD),WORK(KSCRM),ISYMM2,XLAMD1,ISYML1, 292 * D,ISYMD,XLAMD4,ISYML4,XLAMD2,ISYML2,ONE,IOPBFAO) 293 294 ENDIF 295C 296C ----------------------------------------------------------- 297C for IOPT=2,6 include also two-index backtransformed 298C M intermediate (transformed with XLAMD1, XLAMD2) 299C ----------------------------------------------------------- 300C 301 IF ( (IOPT.EQ.2) .OR. (IOPT.EQ.6) .OR. (IOPT.EQ.8) ) THEN 302 303 CALL CC_BFMIAO(WORK(KMGD),MINT,ISYMIN,XLAMD1,ISYML1, 304 & XLAMD2,ISYML2,D,ISYMD,WORK(KEND1),LWRK1) 305 306 IF (IOPT.EQ.6) THEN 307 CALL CC_BFMIAO(WORK(KMGD),MINT,ISYMIN,XLAMD2,ISYML2, 308 & XLAMD1,ISYML1,D,ISYMD,WORK(KEND1),LWRK1) 309 END IF 310 311 END IF 312C 313C ----------------------------------------------------------- 314C save effective density on file: 315C ----------------------------------------------------------- 316C 317 IADRBF(IDEL,IVEC) = IADR 318 319 DTIME = SECOND() 320 CALL PUTWA2(LUBFDEN,FILNAM,WORK(KMGD),IADR,NMGD) 321 TIMIO = TIMIO + SECOND() - DTIME 322 323 IADR = IADR + NMGD 324 325 END DO 326 END DO 327C 328 RETURN 329 END 330*=====================================================================* 331*=====================================================================* 332C /* Deck cc_bfdenf */ 333 SUBROUTINE CC_BFDENF( ZETA2, ISYCTR, MINT, ISYMIN, 334 * XLAMDP, ISYXLP, CHI, ISYCHI, 335 * T1AM, ISYTAM, MGD, WORK, LWORK ) 336*---------------------------------------------------------------------* 337* 338* Purpose: Calculate the effective density used for the BZeta term 339* in the F matrix transformation 340* 341* MGD^{ij}_{alp,k} = Zeta^{k,alp}_{i,j} + M^{k,alp}_{ij} 342* + 2 Chi(k,i) Lambda(alp,j) - Chi(k,j) Lambda(alp,i) 343* 344* Christof Haettig, November 1998 345*=====================================================================* 346 IMPLICIT NONE 347#include "priunit.h" 348#include "ccorb.h" 349#include "ccsdsym.h" 350 351 INTEGER LWORK, ISYCTR, ISYMIN, ISYXLP, ISYCHI, ISYTAM 352 353#if defined (SYS_CRAY) 354 REAL ZETA2(*), MINT(*), CHI(*), T1AM(*), XLAMDP(*) 355 REAL MGD(*), WORK(LWORK) 356 REAL ZERO, ONE, HALF, TWO, DDOT 357#else 358 DOUBLE PRECISION ZETA2(*), MINT(*), CHI(*), T1AM(*), XLAMDP(*) 359 DOUBLE PRECISION MGD(*), WORK(LWORK) 360 DOUBLE PRECISION ZERO, ONE, HALF, TWO, DDOT 361#endif 362 PARAMETER(ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, TWO = 2.0D0) 363 364 INTEGER ISYALP,ISYMM1,ISYMGD,ISYMAK,ISYMI,ISYMK,ISYMJ,ISYMIJ 365 INTEGER KSCRM,KMGD,KOFF1,KOFF2,NIJ,KEND1,LWRK1,IOPT2AO,IOPT 366 INTEGER IALP, ISYMKI 367C 368C----------------------- 369C Symmetry analysis: 370C----------------------- 371C 372 IF ( (ISYMIN.NE.MULD2H(ISYCTR,ISYTAM)) 373 & .OR. (ISYMIN.NE.ISYCHI) ) THEN 374 WRITE (LUPRI,*) 'ISYMIN,ISYCTR,ISYTAM,ISYCHI:', 375 & ISYMIN,ISYCTR,ISYTAM,ISYCHI 376 WRITE (LUPRI,*) 'CC_BFDENF: Symmetry mismatch.' 377 CALL QUIT('CC_BFDENF: Symmetry mismatch.') 378 END IF 379C 380C--------------------------------- 381C Start loop over Alpha index: 382C--------------------------------- 383C 384 DO ISYALP = 1, NSYM 385 386 ISYMM1 = MULD2H(MULD2H(ISYCTR,ISYXLP),ISYALP) 387 ISYMGD = MULD2H(ISYMM1,ISYTAM) 388 389 KSCRM = 1 390 KMGD = KSCRM + NT2BCD(ISYMM1) 391 KEND1 = KMGD + NMAIJK(ISYMGD) 392 LWRK1 = LWORK - KEND1 393 394 IF (LWRK1 .LT. 0) THEN 395 CALL QUIT('Insufficient memory in CC_BFDENF.') 396 END IF 397 398 399 DO A = 1, NBAS(ISYALP) 400C 401C --------------------------------------------------------- 402C calculate one-index AO back-transformed Zeta2 amplitudes: 403C --------------------------------------------------------- 404C 405 IALP = A + IBAS(ISYALP) 406 407 IOPT2AO = 1 408 CALL CC_T2AO(ZETA2,XLAMDP,ISYXLP,WORK(KSCRM),WORK(KEND1), 409 * LWRK1,IALP,ISYALP,ISYCTR,IOPT2AO) 410C 411C --------------------------------------------------------- 412C transform B index with T1AM to occupied and 413C add 2 Chi(ki) Lambda(alp,j) - Chi(kj) Lambda(alp,i) 414C --------------------------------------------------------- 415C 416 IOPT = 4 417 CALL CC_BFMO(WORK(KMGD),WORK(KSCRM),ISYMM1,T1AM,ISYTAM, 418 * A,ISYALP,XLAMDP,ISYXLP,CHI,ISYCHI,ZERO,IOPT) 419C 420C --------------------------------------------------------- 421C add one-index backtransformed M intermediate: 422C --------------------------------------------------------- 423C 424 CALL CC_BFMIMO(WORK(KMGD),MINT,ISYMIN,XLAMDP,ISYXLP,A,ISYALP) 425C 426C ----------------------- 427C sort into result array: 428C ----------------------- 429C 430 DO ISYMJ = 1, NSYM 431 DO ISYMI = 1, NSYM 432 433 ISYMIJ = MULD2H(ISYMI,ISYMJ) 434 ISYMK = MULD2H(ISYMGD,ISYMIJ) 435 ISYMKI = MULD2H(ISYMK,ISYMI) 436 ISYMAK = MULD2H(ISYMK,ISYALP) 437 438 DO J = 1, NRHF(ISYMJ) 439 DO I = 1, NRHF(ISYMI) 440 DO K = 1, NRHF(ISYMK) 441 442 KOFF1 = IMAIJK(ISYMKI,ISYMJ) + NMATIJ(ISYMKI)*(J-1) 443 & + IMATIJ(ISYMK,ISYMI) + NRHF(ISYMK)*(I-1) + K 444 445 NIJ = IMATIJ(ISYMI,ISYMJ) + NRHF(ISYMI)*(J-1) + I 446 447 KOFF2 = IT2AOIJ(ISYMAK,ISYMIJ) + NT1AO(ISYMAK)*(NIJ-1) 448 & + IT1AO(ISYALP,ISYMK) + NBAS(ISYALP)*(K-1) + A 449 450 MGD(KOFF2) = WORK(KMGD-1+KOFF1) 451 452 END DO 453 END DO 454 END DO 455 456 END DO 457 END DO 458 459 END DO 460 END DO 461 462 RETURN 463 END 464*=====================================================================* 465*=====================================================================* 466C /* Deck cc_bfao */ 467 SUBROUTINE CC_BFAO(XMGD,T2AMP,ISYMT,XLAMD0,ISYML0,D,ISYMD, 468 * XLAMD1,ISYML1,XLAMD2,ISYML2,FACT,IOPT) 469*---------------------------------------------------------------------* 470* Purpose: 471* 472* Transform in a three-index batch T2AMP(ci,j) 473* the virtual index 'c' back to AO basis 474* 475* FACT = 0.0 -- initialize XMGD with result 476* FACT = 1.0 -- add result to what is in XMGD 477* 478* IOPT=1,3,5 : add F-term contribution 479* IOPT=2,4,6 : add 2*Coul - Exch combination 480* (else do not add anything) 481* 482*---------------------------------------------------------------------* 483#if defined (IMPLICIT_NONE) 484 IMPLICIT NONE 485#else 486# include "implicit.h" 487#endif 488#include "priunit.h" 489#include "ccorb.h" 490#include "ccsdsym.h" 491C 492#if defined (SYS_CRAY) 493 REAL XMGD(*),T2AMP(*),XLAMD0(*),XLAMD1(*),XLAMD2(*) 494 REAL ONE, TWO, FACT 495#else 496 DOUBLE PRECISION XMGD(*),T2AMP(*),XLAMD0(*),XLAMD1(*),XLAMD2(*) 497 DOUBLE PRECISION ONE, TWO, FACT 498#endif 499 PARAMETER(ONE = 1.0D0, TWO = 2.0D0) 500C 501 INTEGER ISYMT, ISYML0, ISYMD, ISYML1, ISYML2, IOPT 502 INTEGER ISYMJ,ISYMCI,ISYMDJ,ISYMC,ISYMG,ISYMGI,ISYMGJ,ISYMDI 503 INTEGER KOFF1, KOFF2, KOFF3, KOFF4, KOFF5, ISYMI, NVIRC, NBASG 504 INTEGER KOFF6 505C 506 DO ISYMJ = 1,NSYM 507C 508 ISYMCI = MULD2H(ISYMJ,ISYMT) 509 ISYMDJ = MULD2H(ISYMD,ISYMJ) 510C 511 DO ISYMI = 1,NSYM 512C 513 ISYMC = MULD2H(ISYMI,ISYMCI) 514 ISYMG = MULD2H(ISYMC,ISYML0) 515 ISYMGI = MULD2H(ISYMG,ISYMI) 516 ISYMGJ = MULD2H(ISYMG,ISYMJ) 517 ISYMDI = MULD2H(ISYMD,ISYMI) 518C 519c IF (IOPT.EQ.2) THEN 520c KOFF5 = IGLMRH(ISYMD,ISYMJ) + 1 521c WRITE (LUPRI,*) 'XLAMD1: block:',ISYMD,ISYMJ,KOFF5 522c CALL OUTPUT(XLAMD1(KOFF5),1,NBAS(ISYMD),1,NRHF(ISYMJ), 523c & NBAS(ISYMD), NRHF(ISYMJ),1,LUPRI) 524c KOFF5 = IGLMRH(ISYMG,ISYMI) + 1 525c WRITE (LUPRI,*) 'XLAMD2: block:',ISYMG,ISYMI,KOFF5 526c CALL OUTPUT(XLAMD2(KOFF5),1,NBAS(ISYMG),1,NRHF(ISYMI), 527c & NBAS(ISYMG), NRHF(ISYMI),1,LUPRI) 528c END IF 529C 530 NVIRC = MAX(NVIR(ISYMC),1) 531 NBASG = MAX(NBAS(ISYMG),1) 532C 533 KOFF1 = IGLMVI(ISYMG,ISYMC) + 1 534C 535 DO J = 1,NRHF(ISYMJ) 536C 537C -------------------------------------------------- 538C add T2AMP with virtual index backtransformed to AO 539C -------------------------------------------------- 540C 541 KOFF2 = IT2BCD(ISYMCI,ISYMJ) + IT1AM(ISYMC,ISYMI) 542 * + NT1AM(ISYMCI)*(J - 1) + 1 543 KOFF3 = IT2BGD(ISYMGI,ISYMJ) + IT1AO(ISYMG,ISYMI) 544 * + NT1AO(ISYMGI)*(J - 1) + 1 545C 546 CALL DGEMM('N','N',NBAS(ISYMG),NRHF(ISYMI),NVIR(ISYMC), 547 * ONE,XLAMD0(KOFF1),NBASG,T2AMP(KOFF2),NVIRC, 548 * FACT,XMGD(KOFF3),NBASG) 549C 550C ---------------------------------------------------- 551C for IOPT=1,3,5 add F term contribution to density 552C ---------------------------------------------------- 553C 554 IF (IOPT.EQ.1 .OR. IOPT.EQ.3 .OR. IOPT.EQ.5) THEN 555 556 KOFF4 = IGLMRH(ISYMG,ISYMI) + 1 557 KOFF5 = IGLMRH(ISYMD,ISYMJ) + NBAS(ISYMD)*(J-1) + D 558 559 IF (ISYML2.EQ.ISYMGI .AND. ISYML1.EQ.ISYMDJ) THEN 560 CALL DAXPY(NBAS(ISYMG)*NRHF(ISYMI),XLAMD1(KOFF5), 561 * XLAMD2(KOFF4),1,XMGD(KOFF3),1) 562 END IF 563 564 IF (IOPT.EQ.3) THEN 565 IF (ISYML1.EQ.ISYMGI .AND. ISYML2.EQ.ISYMDJ) THEN 566 CALL DAXPY(NBAS(ISYMG)*NRHF(ISYMI),XLAMD2(KOFF5), 567 * XLAMD1(KOFF4),1,XMGD(KOFF3),1) 568 END IF 569 END IF 570 571C 572C ----------------------------------------------- 573C for IOPT=2,4,6,8 add 2*Coul - Exch combination: 574C ----------------------------------------------- 575C 576 ELSE IF (IOPT.EQ.2 .OR. IOPT.EQ.4 .OR. 577 * IOPT.EQ.6 .OR. IOPT.EQ.8 ) THEN 578 579 KOFF4 = IGLMRH(ISYMG,ISYMI) + 1 580 KOFF5 = IGLMRH(ISYMD,ISYMJ) + NBAS(ISYMD)*(J-1) + D 581 582 IF (ISYML2.EQ.ISYMGI .AND. ISYML1.EQ.ISYMDJ) THEN 583 CALL DAXPY(NBAS(ISYMG)*NRHF(ISYMI),TWO*XLAMD1(KOFF5), 584 * XLAMD2(KOFF4),1,XMGD(KOFF3),1) 585 END IF 586 587 IF (ISYML1.EQ.ISYMGI .AND. ISYML2.EQ.ISYMDJ) THEN 588 CALL DAXPY(NBAS(ISYMG)*NRHF(ISYMI),TWO*XLAMD2(KOFF5), 589 * XLAMD1(KOFF4),1,XMGD(KOFF3),1) 590 END IF 591 592 IF (ISYML2.EQ.ISYMGJ .AND. ISYML1.EQ.ISYMDI) THEN 593 DO I = 1, NRHF(ISYMI) 594 595 KOFF4 = IGLMRH(ISYMG,ISYMJ)+ NBAS(ISYMG)*(J-1)+ 1 596 KOFF5 = IGLMRH(ISYMD,ISYMI)+ NBAS(ISYMD)*(I-1)+ D 597 KOFF6 = KOFF3 + NBAS(ISYMG)*(I-1) 598 599 CALL DAXPY(NBAS(ISYMG),-XLAMD1(KOFF5), 600 * XLAMD2(KOFF4),1,XMGD(KOFF6),1) 601 602 END DO 603 END IF 604 605 IF (ISYML1.EQ.ISYMGJ .AND. ISYML2.EQ.ISYMDI) THEN 606 DO I = 1, NRHF(ISYMI) 607 608 KOFF4 = IGLMRH(ISYMG,ISYMJ)+ NBAS(ISYMG)*(J-1)+ 1 609 KOFF5 = IGLMRH(ISYMD,ISYMI)+ NBAS(ISYMD)*(I-1)+ D 610 KOFF6 = KOFF3 + NBAS(ISYMG)*(I-1) 611 612 CALL DAXPY(NBAS(ISYMG),-XLAMD2(KOFF5), 613 * XLAMD1(KOFF4),1,XMGD(KOFF6),1) 614 615 END DO 616 END IF 617C 618 END IF 619 620 END DO 621 END DO 622 END DO 623C 624 RETURN 625 END 626*=====================================================================* 627* END OF SUBROUTINE CC_BFAO * 628*=====================================================================* 629*=====================================================================* 630C /* Deck cc_bfmiao */ 631 SUBROUTINE CC_BFMIAO(XMGD,MINT,ISYMINT,XLAMD1,ISYXL1, 632 & XLAMD2,ISYXL2,D,ISYMD,WORK,LWORK) 633*---------------------------------------------------------------------* 634* Purpose: Transform two indeces of the M intermediate back to AO 635* and add it to the effective density XMGD 636* 637* fixed D index is obtained with XLAMD1 (ISYXL1) 638* free G index is obtained with XLAMD2 (ISYXL2) 639* 640* Christof Haettig, 01.09.1998 641*---------------------------------------------------------------------* 642#if defined (IMPLICIT_NONE) 643 IMPLICIT NONE 644#else 645# include "implicit.h" 646#endif 647#include "priunit.h" 648#include "ccorb.h" 649#include "ccsdsym.h" 650C 651#if defined (SYS_CRAY) 652 REAL XMGD(*), MINT(*), XLAMD1(*), XLAMD2(*), WORK(*) 653 REAL ZERO, ONE, XHALF, XNORM, DDOT 654#else 655 DOUBLE PRECISION XMGD(*), MINT(*), XLAMD1(*), XLAMD2(*), WORK(*) 656 DOUBLE PRECISION ZERO, ONE, XHALF, XNORM, DDOT 657#endif 658 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, XHALF = 0.5D0) 659 660 INTEGER ISYMD,ISYXL1,ISYXL2,ISYMINT, ISYMN, ISYKIJ, ISYMJ, ISYMKI 661 INTEGER ISYKIN, ISYMI, ISYMK, KOFF1, KOFF2, NTOTKI, NTOTD, LWORK 662 INTEGER KOFF4, KOFF5, KOFF6, NTOTG, NTOTK, ISYMG, ISYMGI 663 664 ISYMN = MULD2H(ISYMD,ISYXL1) 665 ISYKIJ = MULD2H(ISYMN,ISYMINT) 666 667 KOFF2 = IGLMRH(ISYMD,ISYMN) + D 668 NTOTD = MAX(NBAS(ISYMD),1) 669 670 IF ( NRHF(ISYMN)*NBAS(ISYMD) .EQ. 0) RETURN 671 672 DO ISYMJ = 1,NSYM 673 674 ISYMKI = MULD2H(ISYKIJ,ISYMJ) 675 ISYKIN = MULD2H(ISYMKI,ISYMN) 676 NTOTKI = MAX(NMATIJ(ISYMKI),1) 677 678 IF (LWORK.LT.NTOTKI) THEN 679 CALL QUIT('Insufficient memory in CC_BFMIAO.') 680 END IF 681 682 DO J = 1, NRHF(ISYMJ) 683 684 KOFF1 = I3ORHF(ISYKIN,ISYMJ) + NMAIJK(ISYKIN)*(J-1) 685 * + IMAIJK(ISYMKI,ISYMN) + 1 686 687 CALL DGEMV('N',NMATIJ(ISYMKI),NRHF(ISYMN),ONE, 688 * MINT(KOFF1),NTOTKI,XLAMD1(KOFF2),NTOTD, 689 * ZERO,WORK,1) 690 691 692 DO ISYMK = 1, NSYM 693 ISYMI = MULD2H(ISYMKI,ISYMK) 694 ISYMG = MULD2H(ISYXL2,ISYMK) 695 ISYMGI = MULD2H(ISYMG,ISYMI) 696 697 KOFF4 = IGLMRH(ISYMG,ISYMK) + 1 698 KOFF5 = IMATIJ(ISYMK,ISYMI) + 1 699 KOFF6 = IT2BGD(ISYMGI,ISYMJ) + NT1AO(ISYMGI)*(J-1) 700 * + IT1AO(ISYMG,ISYMI) + 1 701 702 NTOTG = MAX(NBAS(ISYMG),1) 703 NTOTK = MAX(NRHF(ISYMK),1) 704 705 CALL DGEMM('N','N',NBAS(ISYMG),NRHF(ISYMI),NRHF(ISYMK), 706 * ONE,XLAMD2(KOFF4),NTOTG,WORK(KOFF5),NTOTK, 707 * ONE,XMGD(KOFF6),NTOTG) 708 709 END DO 710 711 END DO 712 713 END DO 714 715 RETURN 716 END 717*=====================================================================* 718* END OF SUBROUTINE CC_BFMIAO * 719*=====================================================================* 720C /* Deck cc_bfmo */ 721 SUBROUTINE CC_BFMO(XMGD,T2AMP,ISYMT,R1AMP,ISYMR,D,ISYMD, 722 * XLAMD,ISYLAM,CHI,ISYCHI,FACT,IOPT) 723*---------------------------------------------------------------------* 724* Purpose: 725* 726* Transform in a three-index batch (ci,j) of T2 amplitudes 727* the virtual index 'c' to occupied using response R1 amplitudes 728* 729* FACT=0.0 : initialize XMGD 730* FACT=1.0 : add contribution to what is already in XMGD 731* 732* IOPT=1,3 : add F-term contribution 733* IOPT=2,4 : add 2*Coul - Exch combination 734* (else do not add anything) 735* 736* Christof Haettig, 04.07.1998 737*---------------------------------------------------------------------* 738#include "implicit.h" 739#include "priunit.h" 740#include "ccorb.h" 741#include "ccsdsym.h" 742C 743 DIMENSION XMGD(*), T2AMP(*), R1AMP(*), XLAMD(*), CHI(*) 744C 745 PARAMETER(ONE = 1.0D0, TWO = 2.0D0, HALF = 0.5D0) 746C 747 DO ISYMJ = 1,NSYM 748C 749 ISYMCI = MULD2H(ISYMJ,ISYMT) 750 ISYMDJ = MULD2H(ISYMJ,ISYMD) 751C 752 DO ISYMI = 1,NSYM 753C 754 ISYMC = MULD2H(ISYMI,ISYMCI) 755 ISYMK = MULD2H(ISYMC,ISYMR) 756 ISYMKI = MULD2H(ISYMK,ISYMI) 757 ISYMKJ = MULD2H(ISYMK,ISYMJ) 758 ISYMDI = MULD2H(ISYMI,ISYMD) 759C 760 NVIRC = MAX(NVIR(ISYMC),1) 761 NRHFK = MAX(NRHF(ISYMK),1) 762C 763 KOFF1 = IT1AM(ISYMC,ISYMK) + 1 764C 765 DO J = 1,NRHF(ISYMJ) 766C 767C ---------------------------------------------------- 768C add T2AMP with virtual index transformed to occupied 769C ---------------------------------------------------- 770C 771 KOFF2 = IT2BCD(ISYMCI,ISYMJ) + NT1AM(ISYMCI)*(J - 1) 772 * + IT1AM(ISYMC,ISYMI) + 1 773 KOFF3 = IMAIJK(ISYMKI,ISYMJ) + NMATIJ(ISYMKI)*(J - 1) 774 * + IMATIJ(ISYMK,ISYMI) + 1 775C 776 CALL DGEMM('T','N',NRHF(ISYMK),NRHF(ISYMI),NVIR(ISYMC), 777 * -ONE,R1AMP(KOFF1),NVIRC,T2AMP(KOFF2),NVIRC, 778 * FACT,XMGD(KOFF3),NRHFK) 779C 780C ---------------------------------------------------- 781C for IOPT=1,3 add F term contribution to density 782C ---------------------------------------------------- 783C 784 IF (IOPT.EQ.1 .OR. IOPT.EQ.3) THEN 785 786 KOFF4 = IMATIJ(ISYMK,ISYMI) + 1 787 KOFF5 = IGLMRH(ISYMD,ISYMJ) + NBAS(ISYMD)*(J-1) + D 788 789 IF (ISYCHI.EQ.ISYMKI .AND. ISYLAM.EQ.ISYMDJ) THEN 790 CALL DAXPY(NRHF(ISYMK)*NRHF(ISYMI),XLAMD(KOFF5), 791 * CHI(KOFF4),1,XMGD(KOFF3),1) 792 END IF 793 794C 795C ------------------------------------------- 796C for IOPT=2,4 add 2*Coul - Exch combination: 797C ------------------------------------------- 798C 799 ELSE IF (IOPT.EQ.2 .OR. IOPT.EQ.4) THEN 800 801 KOFF4 = IMATIJ(ISYMK,ISYMI) + 1 802 KOFF5 = IGLMRH(ISYMD,ISYMJ) + NBAS(ISYMD)*(J-1) + D 803 804 IF (ISYCHI.EQ.ISYMKI .AND. ISYLAM.EQ.ISYMDJ) THEN 805 CALL DAXPY(NRHF(ISYMK)*NRHF(ISYMI),TWO*XLAMD(KOFF5), 806 * CHI(KOFF4),1,XMGD(KOFF3),1) 807 END IF 808 809 IF (ISYCHI.EQ.ISYMKJ .AND. ISYLAM.EQ.ISYMDI) THEN 810 DO I = 1, NRHF(ISYMI) 811 812 KOFF4 = IMATIJ(ISYMK,ISYMJ)+ NRHF(ISYMK)*(J-1)+ 1 813 KOFF5 = IGLMRH(ISYMD,ISYMI)+ NBAS(ISYMD)*(I-1)+ D 814 KOFF6 = KOFF3 + NRHF(ISYMK)*(I-1) 815 816 CALL DAXPY(NRHF(ISYMK),-XLAMD(KOFF5), 817 * CHI(KOFF4),1,XMGD(KOFF6),1) 818 819 END DO 820 END IF 821C 822 END IF 823 824 END DO 825 END DO 826 END DO 827C 828 RETURN 829 END 830*=====================================================================* 831* END OF SUBROUTINE CC_BFMO * 832*=====================================================================* 833C /* Deck cc_bfmimo */ 834 SUBROUTINE CC_BFMIMO(XMGD,MINT,ISYMINT,XLAMD,ISYML,D,ISYMD) 835*---------------------------------------------------------------------* 836* Purpose: Transform one index of the M intermediate back to AO 837* and add it to the effective density XMGD 838* 839* Christof Haettig, 05.07.1998 840*---------------------------------------------------------------------* 841#if defined (IMPLICIT_NONE) 842 IMPLICIT NONE 843#else 844# include "implicit.h" 845#endif 846#include "priunit.h" 847#include "ccorb.h" 848#include "ccsdsym.h" 849C 850#if defined (SYS_CRAY) 851 REAL XMGD(*), MINT(*), XLAMD(*) 852 REAL ONE, XHALF, XNORM, DDOT 853#else 854 DOUBLE PRECISION XMGD(*), MINT(*), XLAMD(*) 855 DOUBLE PRECISION ONE, XHALF, XNORM, DDOT 856#endif 857 PARAMETER (ONE = 1.0D0, XHALF = 0.5D0) 858 859 INTEGER ISYMD, ISYML, ISYMINT, ISYMN, ISYKIJ, ISYMJ, ISYMKI 860 INTEGER ISYKIN, ISYMI, ISYMK, KOFF1, KOFF2, KOFF3, NTOTKI, NTOTD 861 862 ISYMN = MULD2H(ISYMD,ISYML) 863 ISYKIJ = MULD2H(ISYMN,ISYMINT) 864 865 KOFF2 = IGLMRH(ISYMD,ISYMN) + D 866 NTOTD = MAX(NBAS(ISYMD),1) 867 868 DO ISYMJ = 1,NSYM 869 870 ISYMKI = MULD2H(ISYKIJ,ISYMJ) 871 ISYKIN = MULD2H(ISYMKI,ISYMN) 872 873 NTOTKI = MAX(NMATIJ(ISYMKI),1) 874 875 DO J = 1, NRHF(ISYMJ) 876 877 KOFF1 = I3ORHF(ISYKIN,ISYMJ) + NMAIJK(ISYKIN)*(J-1) 878 * + IMAIJK(ISYMKI,ISYMN) + 1 879 880 KOFF3 = IMAIJK(ISYMKI,ISYMJ) + NMATIJ(ISYMKI)*(J-1) + 1 881 882 CALL DGEMV('N',NMATIJ(ISYMKI),NRHF(ISYMN),XHALF, 883 * MINT(KOFF1),NTOTKI,XLAMD(KOFF2),NTOTD, 884 * ONE,XMGD(KOFF3),1) 885 886 END DO 887 888 END DO 889 890 RETURN 891 END 892*=====================================================================* 893* END OF SUBROUTINE CC_BFMIMO * 894*=====================================================================* 895*====================================================================== 896C /* Deck cc_bfi */ 897 SUBROUTINE CC_BFI(BFI,XINT,ISYDIS,XMGD,ISYMGD,D,ISYMD, 898 & SQRINT,FACTOR,WORK,LWORK) 899*---------------------------------------------------------------------- 900* 901* Purpose: contract effective density with AO integrals to 902* the BF intermediate 903* 904* BFI : BF intermediate (symmetry MULD2H(ISYDIS,ISYMGD)) 905* XINT : delta distribution of AO integrals (symmetry ISYDIS) 906* XMGD : effective density for BF term (symmetry ISYMGD) 907* D : delta index (symmetry of the orbital ISYMD) 908* FACTOR : factor with which the contribution is multiplied 909* before it is added to what is already in BFI 910* 911* Christof Haettig, September 1998 912* based on Henrik Kochs CC_BF1 routine 913* Sonia Coriani, October 1999 914* -- option to handle squared XINT distributions in input 915* (SQRINT = true) 916* 917*====================================================================== 918#include "implicit.h" 919#include "priunit.h" 920#include "iratdef.h" 921#include "ccorb.h" 922#include "ccsdsym.h" 923#include "ccsdinp.h" 924 LOGICAL LO3BF, SQRINT 925 INTEGER IOPT, LWORK, ISYDIS, ISYMGD, ISYMD, IOPTSQ 926 DIMENSION XINT(*), BFI(*), XMGD(*), WORK(LWORK) 927 PARAMETER (LO3BF = .FALSE., IOPT = 0) 928 PARAMETER(ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0) 929C 930 ISYRES = MULD2H(ISYDIS,ISYMGD) 931C 932 DO 100 ISYMIJ = 1,NSYM 933C 934 ISYMAB = MULD2H(ISYMIJ,ISYRES) 935 ISYMG = MULD2H(ISYMAB,ISYDIS) 936C 937 KSCRAB = 1 938 KINDV1 = KSCRAB + N2BST(ISYMAB) 939 KINDV2 = KINDV1 + (NNBST(ISYMAB) - 1)/IRAT + 1 940 KEND1 = KINDV2 + (NNBST(ISYMAB) - 1)/IRAT + 1 941 LWRK1 = LWORK - KEND1 942C 943 IF (LWRK1 .LT. 0) THEN 944 CALL QUIT('Insufficient space in CC_BFI') 945 ENDIF 946C 947C ------------------------ 948C Calculate index vectors: 949C ------------------------ 950C 951 CALL CCSD_INDEX(WORK(KINDV1),WORK(KINDV2),ISYMAB) 952C 953C ----------------------------------- 954C Batching and work space allocation. 955C ----------------------------------- 956C 957 IF (ISYMG .EQ. ISYMD) THEN 958 IMAXG = D 959 ELSE IF (ISYMG .LT. ISYMD) THEN 960 IMAXG = NBAS(ISYMG) 961 ELSE 962 IMAXG = 0 963 ENDIF 964C 965 NSIZE = 2*(NNBST(ISYMAB) + NMIJP(ISYMIJ)) 966C 967 IF (IMAXG.EQ.0 .OR. NSIZE.EQ.0) GOTO 100 968C 969 IF (LWRK1.LT.NSIZE) THEN 970 CALL QUIT('Insufficient memory in CC_BFI.') 971 END IF 972C 973 NMAXG = MIN(IMAXG,LWRK1/NSIZE) 974C 975 NBATCH = (IMAXG - 1)/NMAXG + 1 976C 977 DO 110 IBATCH = 1,NBATCH 978C 979 NUMG = NMAXG 980 IF (IBATCH .EQ. NBATCH) THEN 981 NUMG = IMAXG - NMAXG*(NBATCH - 1) 982 ENDIF 983C 984 IG1 = NMAXG*(IBATCH - 1) + 1 985 IG2 = NMAXG*(IBATCH - 1) + NUMG 986C 987 KINTP = KEND1 988 KINTM = KINTP + NNBST(ISYMAB)*NUMG 989 KT2MP = KINTM + NNBST(ISYMAB)*NUMG 990 KT2MM = KT2MP + NUMG*NMIJP(ISYMIJ) 991 KEND2 = KT2MM + NUMG*NMIJP(ISYMIJ) 992C 993 LWRK2 = LWORK - KEND2 994C 995 IF (LWRK2 .LT. 0) THEN 996 CALL QUIT('Insufficient space in CC_BFI') 997 ENDIF 998C 999C ------------------------ 1000C Construct T2MP and T2MM: 1001C ------------------------ 1002C 1003 CALL CCRHS_TPM(WORK(KT2MP),WORK(KT2MM),XMGD,ISYMIJ, 1004 & IG1,NUMG,ISYMG,IOPT,LO3BF) 1005C 1006C ------------------------ 1007C Construct INTP and INTM: 1008C ------------------------ 1009C 1010 IF (SQRINT) THEN 1011 IOPTSQ = 1 1012 ELSE 1013 IOPTSQ = 0 1014 END IF 1015 1016 CALL CCRHS_IPM1(XINT,WORK(KINTP),WORK(KINTM),WORK(KSCRAB), 1017 * WORK(KINDV1),WORK(KINDV2),ISYMAB,ISYMG, 1018 * NUMG,IG1,IG2,IOPTSQ) 1019C 1020C -------------------- 1021C Scale the diagonals: 1022C -------------------- 1023C 1024 IF ((ISYMG .EQ. ISYMD) .AND. (IBATCH .EQ. NBATCH)) THEN 1025 KOFF = KINTP + NNBST(ISYMAB)*(NUMG - 1) 1026 CALL DSCAL(NNBST(ISYMAB),HALF,WORK(KOFF),1) 1027 ENDIF 1028C 1029C ----------------------------- 1030C Add the B-term contributions: 1031C ----------------------------- 1032C 1033 KOFF1 = IT2ORT(ISYMAB,ISYMIJ) + 1 1034 KOFF2 = NT2ORT(ISYRES) + IT2ORT(ISYMAB,ISYMIJ) + 1 1035C 1036 NUMGM = MAX(NUMG,1) 1037 NTOTAB = MAX(NNBST(ISYMAB),1) 1038C 1039 CALL DGEMM('N','N',NNBST(ISYMAB),NMIJP(ISYMIJ),NUMG, 1040 * FACTOR,WORK(KINTP),NTOTAB,WORK(KT2MP),NUMGM, 1041 * ONE,BFI(KOFF1),NTOTAB) 1042C 1043 CALL DGEMM('N','N',NNBST(ISYMAB),NMIJP(ISYMIJ),NUMG, 1044 * FACTOR,WORK(KINTM),NTOTAB,WORK(KT2MM),NUMGM, 1045 * ONE,BFI(KOFF2),NTOTAB) 1046C 1047 110 CONTINUE 1048C 1049 100 CONTINUE 1050C 1051 RETURN 1052C 1053 END 1054*=====================================================================* 1055*=====================================================================* 1056C /* Deck cc_bfib */ 1057 SUBROUTINE CC_BFIB(BFI,BSRHF,ISYRHF,XMGD,ISYMGD,WORK,LWORK) 1058*---------------------------------------------------------------------* 1059* 1060* Purpose: contract effective density with (**|k delta) integrals 1061* to the BF intermediate in the B matrix transformation 1062* (special version of the CC_BFI routine for B matrix) 1063* 1064* BFI : BF(alp i,kj) intermediate (symmetry ISYRHF x ISYMGD) 1065* BSRHF : (**|k delta) integral distribution (symmetry ISYRHF) 1066* XMGD : effective density for BF term (symmetry ISYMGD) 1067* 1068* Written by Christof Haettig July/October 1998 1069*---------------------------------------------------------------------* 1070#if defined (IMPLICIT_NONE) 1071 IMPLICIT NONE 1072#else 1073# include "implicit.h" 1074#endif 1075#include "priunit.h" 1076#include "ccorb.h" 1077#include "maxorb.h" 1078#include "ccsdsym.h" 1079 1080 INTEGER LWORK, ISYMGD, ISYRHF 1081 1082#if defined (SYS_CRAY) 1083 REAL BSRHF(*), BFI(*), XMGD(*), WORK(LWORK) 1084 REAL ONE, HALF 1085#else 1086 DOUBLE PRECISION BSRHF(*), BFI(*), XMGD(*), WORK(LWORK) 1087 DOUBLE PRECISION ONE, HALF 1088#endif 1089 PARAMETER(HALF = 0.5D0, ONE = 1.0D0) 1090C 1091 INTEGER ISYRES, ISYMIJ, ISYMAK, ISYMG 1092 INTEGER ISYMI, ISYMJ, ISYMGI, ISYM, ISYBET, ICOUNT 1093 INTEGER KT2M, KEND1, LWRK1, KOFF1, KOFF3 1094 INTEGER NGIJ, NIJ, NTOTAK, NBASG 1095 INTEGER NBSRHF(8), IBSRHF(8,8) 1096C 1097C -------------------------------------- 1098C precalculate symmetry array for BSRHF: 1099C -------------------------------------- 1100 DO ISYM = 1, NSYM 1101 ICOUNT = 0 1102 DO ISYMAK = 1, NSYM 1103 ISYBET = MULD2H(ISYMAK,ISYM) 1104 IBSRHF(ISYMAK,ISYBET) = ICOUNT 1105 ICOUNT = ICOUNT + NT1AO(ISYMAK)*NBAS(ISYBET) 1106 END DO 1107 NBSRHF(ISYM) = ICOUNT 1108 END DO 1109C 1110 ISYRES = MULD2H(ISYRHF,ISYMGD) 1111C 1112 DO ISYMIJ = 1, NSYM 1113C 1114 ISYMAK = MULD2H(ISYRES,ISYMIJ) 1115 ISYMG = MULD2H(ISYMAK,ISYRHF) 1116C 1117 KT2M = 1 1118 KEND1 = KT2M + NMATIJ(ISYMIJ)*NBAS(ISYMG) 1119 LWRK1 = LWORK - KEND1 1120 1121 IF (LWRK1 .LT. 0) THEN 1122 CALL QUIT('Insufficient work space in CC_BFIB.') 1123 END IF 1124C 1125C --------------------------------------------------- 1126C resort amplitude density: 1127C --------------------------------------------------- 1128C 1129 DO ISYMJ = 1, NSYM 1130C 1131 ISYMI = MULD2H(ISYMJ,ISYMIJ) 1132 ISYMGI = MULD2H(ISYMI,ISYMG) 1133C 1134 DO J = 1, NRHF(ISYMJ) 1135 DO I = 1, NRHF(ISYMI) 1136C 1137 NGIJ = IT2BGD(ISYMGI,ISYMJ) + NT1AO(ISYMGI)*(J-1) 1138 & + IT1AO(ISYMG,ISYMI) + NBAS(ISYMG)*(I-1) + 1 1139 1140 NIJ = IMATIJ(ISYMI,ISYMJ) + NRHF(ISYMI)*(J-1) + I 1141 1142 KOFF1 = KT2M + NBAS(ISYMG)*(NIJ-1) 1143 1144 CALL DCOPY(NBAS(ISYMG),XMGD(NGIJ),1,WORK(KOFF1),1) 1145 1146 END DO 1147 END DO 1148 1149 END DO 1150C 1151C -------------------------------------- 1152C contract with the presorted integrals: 1153C -------------------------------------- 1154C 1155 KOFF1 = IT2AOIJ(ISYMAK,ISYMIJ) + 1 1156 KOFF3 = IBSRHF(ISYMAK,ISYMG) + 1 1157 NTOTAK = MAX(NT1AO(ISYMAK),1) 1158 NBASG = MAX(NBAS(ISYMG),1) 1159C 1160 CALL DGEMM('N','N',NT1AO(ISYMAK),NMATIJ(ISYMIJ),NBAS(ISYMG), 1161 & ONE,BSRHF(KOFF3),NTOTAK,WORK(KT2M),NBASG, 1162 & ONE,BFI(KOFF1), NTOTAK) 1163C 1164 END DO 1165C 1166 RETURN 1167 END 1168*=====================================================================* 1169*=====================================================================* 1170C /* Deck cc_bfbsort */ 1171 SUBROUTINE CC_BFBSORT(DSRHF,BSRHF,ISYRHF) 1172*---------------------------------------------------------------------* 1173* 1174* Purpose: presort DSRHF integral array for the BF intermediate 1175* calculation in the B matrix transformation 1176* 1177* DSRHF : (alp bet|k delta) integrals for a fixed delta 1178* BSRHF : integrals sorted as I(alp k;bet)^del 1179* ISYRHF : symmetry of the integral arrays 1180* 1181* Written by Christof Haettig July/October 1998 1182*---------------------------------------------------------------------* 1183 1184 use dyn_iadrpk 1185 1186#if defined (IMPLICIT_NONE) 1187 IMPLICIT NONE 1188#else 1189# include "implicit.h" 1190#endif 1191#include "priunit.h" 1192#include "ccorb.h" 1193#include "maxorb.h" 1194#include "ccsdsym.h" 1195#include "symsq.h" 1196 1197 INTEGER ISYRHF, ISYM, ISYMAK, ISYBET, ISYMK, ISYMAB, ISYALP 1198 INTEGER ICOUNT, NBSRHF(8), IBSRHF(8,8) 1199 INTEGER NABK, NAKB, NAK, KOFF1, IJSQ 1200 1201#if defined (SYS_CRAY) 1202 REAL DSRHF(*), BSRHF(*) 1203#else 1204 DOUBLE PRECISION DSRHF(*), BSRHF(*) 1205#endif 1206C 1207C -------------------------------------- 1208C precalculate symmetry array for BSRHF: 1209C -------------------------------------- 1210C 1211 DO ISYM = 1, NSYM 1212 ICOUNT = 0 1213 DO ISYMAK = 1, NSYM 1214 ISYBET = MULD2H(ISYMAK,ISYM) 1215 IBSRHF(ISYMAK,ISYBET) = ICOUNT 1216 ICOUNT = ICOUNT + NT1AO(ISYMAK)*NBAS(ISYBET) 1217 END DO 1218 NBSRHF(ISYM) = ICOUNT 1219 END DO 1220C 1221C ------------------- 1222C sort the integrals: 1223C ------------------- 1224C 1225 DO ISYMAK = 1, NSYM 1226 DO ISYMK = 1, NSYM 1227C 1228 ISYBET = MULD2H(ISYMAK,ISYRHF) 1229 ISYALP = MULD2H(ISYMK,ISYMAK) 1230 ISYMAB = MULD2H(ISYALP,ISYBET) 1231C 1232C -------------------------------------------------------- 1233C get (alp k;bet) blocks out of (alp bet|k del) integrals: 1234C -------------------------------------------------------- 1235C 1236 DO K = 1, NRHF(ISYMK) 1237C 1238 KOFF1 = IDSRHF(ISYMAB,ISYMK) + NNBST(ISYMAB)*(K-1) 1239C 1240 DO A = 1, NBAS(ISYALP) 1241 DO B = 1, NBAS(ISYBET) 1242C 1243 IJSQ = IAODIS(ISYALP,ISYBET) + NBAS(ISYALP)*(B-1) + A 1244 NABK = KOFF1 + IADRPK( I2BST(ISYMAB) + IJSQ ) 1245 NAK = IT1AO(ISYALP,ISYMK) + NBAS(ISYALP)*(K-1) + A 1246 NAKB = IBSRHF(ISYMAK,ISYBET) +NT1AO(ISYMAK)*(B-1) + NAK 1247C 1248 BSRHF(NAKB) = DSRHF(NABK) 1249C 1250 END DO 1251 END DO 1252C 1253 END DO 1254C 1255 END DO 1256 END DO 1257C 1258 RETURN 1259 END 1260*=====================================================================* 1261*=====================================================================* 1262C /* Deck ccrhs_tpm */ 1263 SUBROUTINE CCRHS_TPM(XMP,XMM,XMGD,ISYMIJ, 1264 & IG1,NUMG,ISYMG,IOPT,LO3BF) 1265*---------------------------------------------------------------------* 1266* 1267* Purpose: construct M+ and M- 'densities' for a batch of G 1268* 1269*---------------------------------------------------------------------* 1270#include "implicit.h" 1271#include "priunit.h" 1272#include "ccorb.h" 1273#include "ccsdsym.h" 1274#include "ccsdinp.h" 1275C 1276 PARAMETER(ONE = 1.0D0, THREE = 3.0D0) 1277C 1278 LOGICAL LO3BF 1279 DIMENSION XMP(*), XMM(*), XMGD(*) 1280C 1281 INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 1282C 1283 DO ISYMJ = 1,NSYM 1284C 1285 ISYMI = MULD2H(ISYMJ,ISYMIJ) 1286 ISYMGI = MULD2H(ISYMI,ISYMG) 1287 ISYMGJ = MULD2H(ISYMJ,ISYMG) 1288C 1289 IF (ISYMI .LE. ISYMJ) THEN 1290C 1291 NTOTI = NRHF(ISYMI) 1292C 1293 DO J = 1,NRHF(ISYMJ) 1294C 1295 IF (ISYMI .EQ. ISYMJ) NTOTI = J 1296C 1297 DO I = 1,NTOTI 1298C 1299 IF (LO3BF .AND. IOPT.EQ.4) THEN 1300 NGIJ = IMAIJK(ISYMGI,ISYMJ) + NMATIJ(ISYMGI)*(J - 1) 1301 * + IMATIJ(ISYMG,ISYMI) + NRHF(ISYMG)*(I-1) + IG1 1302C 1303 NGJI = IMAIJK(ISYMGJ,ISYMI) + NMATIJ(ISYMGJ)*(I - 1) 1304 * + IMATIJ(ISYMG,ISYMJ) + NRHF(ISYMG)*(J-1) + IG1 1305 ELSE 1306 NGIJ = IT2BGD(ISYMGI,ISYMJ) + NT1AO(ISYMGI)*(J - 1) 1307 * + IT1AO(ISYMG,ISYMI) + NBAS(ISYMG)*(I-1) + IG1 1308C 1309 NGJI = IT2BGD(ISYMGJ,ISYMI) + NT1AO(ISYMGJ)*(I - 1) 1310 * + IT1AO(ISYMG,ISYMJ) + NBAS(ISYMG)*(J-1) + IG1 1311 END IF 1312C 1313 IF (ISYMI .EQ. ISYMJ) THEN 1314 NIJ = IMIJP(ISYMI,ISYMJ) + INDEX(I,J) 1315 ELSE 1316 NIJ = IMIJP(ISYMI,ISYMJ) + NRHF(ISYMI)*(J - 1) + I 1317 ENDIF 1318C 1319 NGIJPM = NUMG*(NIJ - 1) + 1 1320C 1321 IF ( CC2 ) THEN 1322 CALL DZERO(XMP(NGIJPM),NUMG) 1323 CALL DZERO(XMM(NGIJPM),NUMG) 1324 ELSE 1325 CALL DCOPY(NUMG,XMGD(NGIJ),1,XMP(NGIJPM),1) 1326 CALL DCOPY(NUMG,XMGD(NGIJ),1,XMM(NGIJPM),1) 1327C 1328 CALL DAXPY(NUMG, ONE,XMGD(NGJI),1,XMP(NGIJPM),1) 1329 CALL DAXPY(NUMG,-ONE,XMGD(NGJI),1,XMM(NGIJPM),1) 1330 ENDIF 1331 1332 END DO 1333 1334 END DO 1335 1336 END IF 1337 1338 END DO 1339C 1340 RETURN 1341 END 1342*=====================================================================* 1343* END OF SUBROUTINE CCRHS_TPM * 1344*=====================================================================* 1345C /* Deck cc_bfmf */ 1346 SUBROUTINE CC_BFMF(XMP,XMM,NUMG,FACT,IOPT,LO3BF, 1347 & I,ISYMI,J,ISYMJ,IG1,ISYMG,D,ISYMD, 1348 & XLAMD1,ISYML1,XLAMD2,ISYML2) 1349*---------------------------------------------------------------------* 1350* 1351* Purpose: add F-term 'density' to M+ and M- arrays 1352* 1353*---------------------------------------------------------------------* 1354#include "implicit.h" 1355#include "priunit.h" 1356#include "ccorb.h" 1357#include "ccsdsym.h" 1358C 1359 LOGICAL LO3BF 1360 DIMENSION XMP(*), XMM(*), XLAMD1(*), XLAMD2(*) 1361C 1362 IF ((ISYMJ .EQ. MULD2H(ISYML2,ISYMD)) .AND. 1363 & (ISYMI .EQ. MULD2H(ISYML1,ISYMG)) ) THEN 1364C 1365 KOFF1 = IGLMRH(ISYMD,ISYMJ) + NBAS(ISYMD)*(J - 1) + D 1366C 1367 IF (LO3BF .AND. IOPT.EQ.4) THEN 1368 KOFF2 = IMATIJ(ISYMG,ISYMI) + NRHF(ISYMG)*(I - 1) + IG1 1369 ELSE 1370 KOFF2 = IGLMRH(ISYMG,ISYMI) + NBAS(ISYMG)*(I - 1) + IG1 1371 END IF 1372C 1373 CALL DAXPY(NUMG, XLAMD2(KOFF1),XLAMD1(KOFF2),1,XMP,1) 1374 CALL DAXPY(NUMG,FACT*XLAMD2(KOFF1),XLAMD1(KOFF2),1,XMM,1) 1375C 1376 ENDIF 1377C 1378 IF ((ISYMI .EQ. MULD2H(ISYML2,ISYMD)) .AND. 1379 & (ISYMJ .EQ. MULD2H(ISYML1,ISYMG)) ) THEN 1380C 1381 KOFF1 = IGLMRH(ISYMD,ISYMI) + NBAS(ISYMD)*(I - 1) + D 1382C 1383 IF (LO3BF .AND. IOPT.EQ.4) THEN 1384 KOFF2 = IMATIJ(ISYMG,ISYMJ) + NRHF(ISYMG)*(J - 1) + IG1 1385 ELSE 1386 KOFF2 = IGLMRH(ISYMG,ISYMJ) + NBAS(ISYMG)*(J - 1) + IG1 1387 END IF 1388C 1389 CALL DAXPY(NUMG, XLAMD2(KOFF1),XLAMD1(KOFF2),1,XMP,1) 1390 CALL DAXPY(NUMG,-FACT*XLAMD2(KOFF1),XLAMD1(KOFF2),1,XMM,1) 1391C 1392 ENDIF 1393C 1394 RETURN 1395 END 1396*=====================================================================* 1397* END OF SUBROUTINE CC_BFMF * 1398*=====================================================================* 1399