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_21i2 */ 21 SUBROUTINE CC_21I2(RHO1, XINT, ISYDIS, XINT2, ISYDIS2, 22 * PINT0,QINT0,ISYVWZ0,PINT2,QINT2,ISYVWZ2, 23 * XLAMP0,XLAMH0,ISYXL0,XLAMP2,ISYXL2, 24 * WORK,LWORK,IOPT,LRHOAO,LDERIV,LX2ISQ) 25C 26C Written by Asger Halkier 31/10 - 1995 27C restructured by Christof Haettig July 1998 28C X2INT already squared in input, Sonia Coriani November 1999 29C 30C Version: 2.0 31C 32C XINT, ISYDIS -- (* *|* delta) batch of integrals, its symmetry 33C XINT2, ISYDIS2 -- the same for derivative integrals (for IOPT=3) 34C PINT0, QINT0 -- P, Q linear combination of ZWV intermediates 35C PINT2, QINT2 -- P, Q calculated from response amplitudes 36C 37C IOPT = 1 : calculate only usual LHTR contributions 38C 2 : include response contrib. from PINT2, QINT2 (F mat.) 39C 3 : include additional orbital relaxation 40C and derivative contributions (for ETA/RHS vectors) 41C 42C LRHOAO = .TRUE. return result with a index in AO 43C LRHOAO = .FALSE. return result with a index in MO 44C LDERIV = .FALSE. omit contribution from derivative integrals 45C LX2ISQ = .TRUE. (a b|* *) of XINT2 is a full matrix (fx LAO) 46C 47C Purpose: To calculate the 21I contribution to RHO1! 48C 49#if defined (IMPLICIT_NONE) 50 IMPLICIT NONE 51#else 52# include "implicit.h" 53#endif 54#include "priunit.h" 55#include "ccorb.h" 56#include "ccsdsym.h" 57C 58 LOGICAL LRHOAO, LDERIV, LX2ISQ 59 INTEGER LWORK,ISYDIS,ISYDIS2,ISYVWZ0,ISYVWZ2,ISYXL0,ISYXL2,IOPT 60C 61#if defined (SYS_CRAY) 62 REAL RHO1(*), PINT0(*), PINT2(*), QINT0(*), QINT2(*), 63 * XINT(*), XINT2(*), WORK(LWORK), 64 * XLAMP0(*), XLAMH0(*),XLAMP2(*) 65 REAL DUMMY, ZERO, ONE, TWO 66#else 67 DOUBLE PRECISION RHO1(*), PINT0(*), PINT2(*), QINT0(*), QINT2(*), 68 * XINT(*), XINT2(*), WORK(LWORK), 69 * XLAMP0(*), XLAMH0(*),XLAMP2(*) 70 DOUBLE PRECISION DUMMY, ZERO, ONE, TWO 71#endif 72C 73 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 74C 75 INTEGER NT3BAO(8), IT3BAO(8,8) 76 INTEGER ISYGIJ, ISYRES, ISY3AO, ISYMG 77 INTEGER ISYALI, ISYDEN 78 INTEGER KRHOAO, KPQAO, KPQMO, KPQDEN, NPQMO 79 INTEGER KEND1,LWRK1 80 INTEGER ICOUNT 81 INTEGER KPQAO0, KPQDEN0, ISYGIJ0, ISYDEN0, IOPTDEN, IOPTCON 82C 83C set some symmetries use globally in this routine: 84C 85 ISYGIJ = MULD2H(ISYVWZ0,ISYXL2) 86 ISYDEN = MULD2H(ISYGIJ,ISYXL0) 87 ISYRES = MULD2H(ISYDEN,ISYDIS) 88 ISYGIJ0 = MULD2H(ISYVWZ0,ISYXL0) 89 ISYDEN0 = MULD2H(ISYGIJ0,ISYXL0) 90C 91 IF ( IOPT.GE.2 .AND. MULD2H(ISYVWZ2,ISYXL0).NE.ISYGIJ ) THEN 92 CALL QUIT('Symmetry mismatch in CC_21I2.') 93 END IF 94C 95C calculate index arrays for intermediates with 3 indeces in AO: 96C 97 DO ISY3AO = 1, NSYM 98 ICOUNT = 0 99 DO ISYMG = 1, NSYM 100 ISYALI = MULD2H(ISY3AO,ISYMG) 101 IT3BAO(ISYALI,ISYMG) = ICOUNT 102 ICOUNT = ICOUNT + NBAS(ISYMG) * NT1AO(ISYALI) 103 END DO 104 NT3BAO(ISY3AO) = ICOUNT 105 END DO 106C 107C---------------------------------------------------------------- 108C allocate work space for 109C -- intermediate with one index backtransformed 110C -- intermediate with one more index backtransformed 111C -- intermediate with two more indeces backtransformed 112C---------------------------------------------------------------- 113C 114 NPQMO = NT2BCD(ISYVWZ0) 115 IF (IOPT.GE.2) NPQMO = MAX(NPQMO,NT2BCD(ISYVWZ2)) 116C 117 KRHOAO = 1 118 KPQMO = KRHOAO + NT1AO(ISYRES) 119 KPQAO = KPQMO + NPQMO 120 KPQDEN = KPQAO + NT2BGD(ISYGIJ) 121 KEND1 = KPQDEN + NT3BAO(ISYDEN) 122 LWRK1 = LWORK - KEND1 123C 124 IF (IOPT.EQ.3) THEN 125 KPQAO0 = KEND1 126 KEND1 = KPQAO0 + NT2BGD(ISYGIJ0) 127 IF (LDERIV) THEN 128 KPQDEN0 = KEND1 129 KEND1 = KPQDEN0 + NT3BAO(ISYDEN0) 130 END IF 131 LWRK1 = LWORK - KEND1 132 END IF 133C 134 IF ( LWRK1 .LT. 0) THEN 135 CALL QUIT('Insufficient work space in CC_21I.') 136 ENDIF 137C 138C-------------------------------------------------------------------- 139C 1. Exchange like contribution: add P and Q intermediates 140C together and transform virtual index back to AO 141C-------------------------------------------------------------------- 142C 143 CALL DCOPY(NT2BCD(ISYVWZ0),PINT0,1,WORK(KPQMO),1) 144 CALL DAXPY(NT2BCD(ISYVWZ0),ONE,QINT0,1,WORK(KPQMO),1) 145C 146 CALL CC_BFAO(WORK(KPQAO),WORK(KPQMO),ISYVWZ0,XLAMP2,ISYXL2, 147 & 1,1,DUMMY,1,DUMMY,1,ZERO,0) 148C 149c XNORM = DDOT(NT2BGD(ISYGIJ),WORK(KPQAO),1,WORK(KPQAO),1) 150c WRITE (LUPRI,*) 'Norm(PQAO)1=',XNORM 151C 152C-------------------------------------------------------------------- 153C if IOPT = 3 calculate PQAO0, backtransformed with XLAMP0 154C for if IOPT = 2/3 include response contribution: 155C-------------------------------------------------------------------- 156C 157 IF (IOPT.EQ.3) THEN 158 CALL CC_BFAO(WORK(KPQAO0),WORK(KPQMO),ISYVWZ0,XLAMP0,ISYXL0, 159 & 1,1,DUMMY,1,DUMMY,1,ZERO,0) 160 END IF 161C 162 IF (IOPT.EQ.2 .OR. IOPT.EQ.3) THEN 163 164 CALL DCOPY(NT2BCD(ISYVWZ2),PINT2,1,WORK(KPQMO),1) 165 CALL DAXPY(NT2BCD(ISYVWZ2),ONE,QINT2,1,WORK(KPQMO),1) 166C 167 CALL CC_BFAO(WORK(KPQAO),WORK(KPQMO),ISYVWZ2,XLAMP0,ISYXL0, 168 & 1,1,DUMMY,1,DUMMY,1,ONE,0) 169 170 END IF 171C 172c XNORM = DDOT(NT2BGD(ISYGIJ),WORK(KPQAO),1,WORK(KPQAO),1) 173c WRITE (LUPRI,*) 'Norm(PQAO)2=',XNORM 174C 175C-------------------------------------------------------------------- 176C transform occupied index j back to AO (alpha) 177C for IOPT=3 add extra relaxation and transform also 178C the j index of the KPQAO0 intermediate 179C-------------------------------------------------------------------- 180C 181 IOPTDEN = 0 182 CALL CC_21IDEN(WORK(KPQAO),ISYGIJ,XLAMP0,ISYXL0,WORK(KPQDEN), 183 * ISYDEN,IT3BAO,WORK(KEND1),LWRK1,'EXC',IOPTDEN) 184 185 IF (IOPT.EQ.3) THEN 186 IOPTDEN = 1 187 CALL CC_21IDEN(WORK(KPQAO0),ISYGIJ0,XLAMP2,ISYXL2, 188 * WORK(KPQDEN),ISYDEN,IT3BAO, 189 * WORK(KEND1),LWRK1,'EXC',IOPTDEN) 190 191 IF (LDERIV) THEN 192 IOPTDEN = 0 193 CALL CC_21IDEN(WORK(KPQAO0),ISYGIJ0,XLAMP0,ISYXL0, 194 * WORK(KPQDEN0),ISYDEN,IT3BAO, 195 * WORK(KEND1),LWRK1,'EXC',IOPTDEN) 196 END IF 197 END IF 198C 199c XNORM = DDOT(NT3BAO(ISYGIJ),WORK(KPQDEN),1,WORK(KPQDEN),1) 200c WRITE (LUPRI,*) 'Norm(PQDEN)1=',XNORM 201C 202C-------------------------------------------------------------------- 203C 2. Coulomb like contribution: scale P intermediate with -2 204C and transform virutal e back to AO 205C-------------------------------------------------------------------- 206C 207 CALL DCOPY(NT2BCD(ISYVWZ0),PINT0,1,WORK(KPQMO),1) 208 CALL DSCAL(NT2BCD(ISYVWZ0),-TWO,WORK(KPQMO),1) 209 210 CALL CC_BFAO(WORK(KPQAO),WORK(KPQMO),ISYVWZ0,XLAMP2,ISYXL2, 211 & 1,1,DUMMY,1,DUMMY,1,ZERO,0) 212C 213c XNORM = DDOT(NT2BGD(ISYGIJ),WORK(KPQAO),1,WORK(KPQAO),1) 214c WRITE (LUPRI,*) 'Norm(PQAO)3=',XNORM 215C 216C-------------------------------------------------------------------- 217C if IOPT = 3 calculate PQAO0, backtransformed with XLAMP0 218C for IOPT = 2/3 include response contribution: 219C-------------------------------------------------------------------- 220C 221 IF (IOPT.EQ.3) THEN 222 CALL CC_BFAO(WORK(KPQAO0),WORK(KPQMO),ISYVWZ0,XLAMP0,ISYXL0, 223 & 1,1,DUMMY,1,DUMMY,1,ZERO,0) 224 END IF 225C 226 IF (IOPT.EQ.2 .OR. IOPT.EQ.3) THEN 227 228 CALL DCOPY(NT2BCD(ISYVWZ2),PINT2,1,WORK(KPQMO),1) 229 CALL DSCAL(NT2BCD(ISYVWZ2),-TWO,WORK(KPQMO),1) 230 231 CALL CC_BFAO(WORK(KPQAO),WORK(KPQMO),ISYVWZ2,XLAMP0,ISYXL0, 232 & 1,1,DUMMY,1,DUMMY,1,ONE,0) 233 234 END IF 235C 236c XNORM = DDOT(NT2BGD(ISYGIJ),WORK(KPQAO),1,WORK(KPQAO),1) 237c WRITE (LUPRI,*) 'Norm(PQAO)4=',XNORM 238C 239C-------------------------------------------------------------------- 240C transform occupied index j back to AO (gamma), add result to 241C the effective density stored in KPQDEN 242C-------------------------------------------------------------------- 243C 244 IOPTDEN = 1 245 CALL CC_21IDEN(WORK(KPQAO),ISYGIJ,XLAMP0,ISYXL0,WORK(KPQDEN), 246 * ISYDEN,IT3BAO,WORK(KEND1),LWRK1,'COU',IOPTDEN) 247 248 IF (IOPT.EQ.3) THEN 249 IOPTDEN = 1 250 CALL CC_21IDEN(WORK(KPQAO0),ISYGIJ0,XLAMP2,ISYXL2, 251 * WORK(KPQDEN),ISYDEN,IT3BAO, 252 * WORK(KEND1),LWRK1,'COU',IOPTDEN) 253 254 IF (LDERIV) THEN 255 IOPTDEN = 1 256 CALL CC_21IDEN(WORK(KPQAO0),ISYGIJ0,XLAMP0,ISYXL0, 257 * WORK(KPQDEN0),ISYDEN,IT3BAO, 258 * WORK(KEND1),LWRK1,'COU',IOPTDEN) 259 END IF 260 END IF 261C 262C-------------------------------------------------------------------- 263C contract the effective density with the 2-el integrals: 264C-------------------------------------------------------------------- 265C 266c XNORM = DDOT(NT3BAO(ISYGIJ),WORK(KPQDEN),1,WORK(KPQDEN),1) 267c WRITE (LUPRI,*) 'Norm(PQDEN)2=',XNORM 268C 269 CALL DZERO(WORK(KRHOAO),NT1AO(ISYRES)) 270 271 CALL CC_21ICON(WORK(KRHOAO),ISYRES,WORK(KPQDEN),ISYDEN, 272 & XINT,ISYDIS,IT3BAO,WORK(KEND1),LWRK1,1) 273 274 IF (IOPT.EQ.3 .AND. LDERIV) THEN 275 276 IF (LX2ISQ) THEN 277 IOPTCON = 2 !X2INT has full (a b| already 278 ELSE 279 IOPTCON = 1 !X2INT is squared inside CC_21ICON 280 END IF 281 282 CALL CC_21ICON(WORK(KRHOAO),ISYRES,WORK(KPQDEN0),ISYDEN0, 283 & XINT2,ISYDIS2,IT3BAO,WORK(KEND1),LWRK1,IOPTCON) 284 END IF 285C 286c XNORM = DDOT(NT1AO(ISYRES),RHO1,1,RHO1,1) 287c WRITE (LUPRI,*) 'XNORM=',XNORM 288c 289C--------------------------------------------- 290C transformation and/or storage in result: 291C--------------------------------------------- 292C 293 IF (LRHOAO) THEN 294 CALL DAXPY(NT1AO(ISYRES),ONE,WORK(KRHOAO),1,RHO1,1) 295 ELSE 296 CALL CC_T1AM(RHO1,ISYRES,WORK(KRHOAO),ISYRES,XLAMH0,ISYXL0,ONE) 297 END IF 298C 299 RETURN 300 END 301C 302*=====================================================================* 303C /* Deck cc_t1am */ 304 SUBROUTINE CC_T1AM(RHO1MO,ISYMO,RHO1AO,ISYAO,XLAMD,ISYLAM,FAC) 305C 306C Purpose: transform AO index of a two-index array with 307C IT1AO structure to MO and store in IT1AM structure 308C 309C RHO1MO(ai) = FAC*RHO1MO(ai) + sum_alp XLAMD(alp a) RHO1AO(alp i) 310C 311C Written by Christof Haettig 16 July 1998 312C 313#include "implicit.h" 314#include "priunit.h" 315#include "ccorb.h" 316#include "ccsdsym.h" 317C 318 DIMENSION RHO1MO(*), RHO1AO(*), XLAMD(*) 319C 320 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 321C 322 IF (MULD2H(ISYAO,ISYLAM).NE.ISYMO ) THEN 323 CALL QUIT('Symmetry mismatch in CC_T1AM.') 324 END IF 325C 326 DO ISYMA = 1, NSYM 327C 328 ISYMI = MULD2H(ISYMO,ISYMA) 329 ISYMBE = MULD2H(ISYLAM,ISYMA) 330C 331 KOFF1 = IGLMVI(ISYMBE,ISYMA) + 1 332 KOFF2 = IT1AO(ISYMBE,ISYMI) + 1 333 KOFF3 = IT1AM(ISYMA,ISYMI) + 1 334C 335 NTOTA = MAX(NVIR(ISYMA),1) 336 NTOTBE = MAX(NBAS(ISYMBE),1) 337C 338 CALL DGEMM('T','N',NVIR(ISYMA),NRHF(ISYMI),NBAS(ISYMBE), 339 * ONE,XLAMD(KOFF1),NTOTBE,RHO1AO(KOFF2),NTOTBE, 340 * FAC,RHO1MO(KOFF3),NTOTA) 341C 342 END DO 343C 344 RETURN 345 END 346*=====================================================================* 347C /* Deck cc_21iden */ 348 SUBROUTINE CC_21IDEN(PQAO,ISYGIJ,XLAMDP,ISYMXL,PQDEN,ISYDEN, 349 * IT3BAO,WORK,LWORK,FLAG,IOPT) 350C 351C Purpose: transform and resort the two-index backtransformed 352C P+Q intermediate used in CC_21I to an density like 353C intermediate with three indeces in AO and one in MO 354C 355C FLAG = 'EXC' do transformation/resort for exchange like contrib. 356C 'COU' do transformation for coulomb like contrib. 357C 358C IOPT = 0 intialize PQDEN with the actual contribution 359C 1 add the actual contribution to what is already in PQDEN 360C 361C Written by Christof Haettig July 1998 362C 363#if defined (IMPLICIT_NONE) 364 IMPLICIT NONE 365#else 366# include "implicit.h" 367#endif 368#include "priunit.h" 369#include "ccorb.h" 370#include "ccsdsym.h" 371C 372 CHARACTER*(*) FLAG 373 INTEGER LWORK, ISYGIJ, ISYDEN, ISYMXL, IOPT, IT3BAO(8,8) 374C 375#if defined (SYS_CRAY) 376 REAL PQAO(*), PQDEN(*), XLAMDP(*), WORK(LWORK) 377 REAL ZERO, ONE, FACT 378#else 379 DOUBLE PRECISION PQAO(*), PQDEN(*), XLAMDP(*), WORK(LWORK) 380 DOUBLE PRECISION ZERO, ONE, FACT 381#endif 382C 383 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 384C 385 INTEGER ISYMJ,ISYMGI,ISYMAL,ISYMI,ISYALI,ISYGAM,KOFF1,KOFF2,KOFF3 386 INTEGER NTOTGA, NTOTAL, NTOTGI, NGI, NAGI, NAIG, KOFF4, ISYMG 387C 388C-------------------------------------------------------------------- 389C transform occupied index j back to AO (alpha) 390C-------------------------------------------------------------------- 391C 392 IF (FLAG(1:3).EQ.'COU') THEN 393 394 FACT = ZERO 395 IF (IOPT.EQ.1) FACT = ONE 396 397 DO ISYMJ = 1, NSYM 398 399 ISYMGI = MULD2H(ISYGIJ,ISYMJ) 400 ISYGAM = MULD2H(ISYMXL,ISYMJ) 401 402 KOFF1 = IGLMRH(ISYGAM,ISYMJ) + 1 403 KOFF2 = IT2BGD(ISYMGI,ISYMJ) + 1 404 KOFF3 = IT3BAO(ISYMGI,ISYGAM) + 1 405 406 NTOTGA = MAX(NBAS(ISYGAM),1) 407 NTOTGI = MAX(NT1AO(ISYMGI),1) 408 409 CALL DGEMM('N','T',NT1AO(ISYMGI),NBAS(ISYGAM),NRHF(ISYMJ), 410 * ONE,PQAO(KOFF2),NTOTGI,XLAMDP(KOFF1),NTOTGA, 411 * FACT,PQDEN(KOFF3),NTOTGI) 412 413 END DO 414C 415C-------------------------------------------------------------------- 416C transform occupied index j back to AO (alpha) 417C-------------------------------------------------------------------- 418C 419 ELSE IF (FLAG(1:3).EQ.'EXC') THEN 420 421 DO ISYMJ = 1, NSYM 422 423 ISYMGI = MULD2H(ISYGIJ,ISYMJ) 424 ISYMAL = MULD2H(ISYMXL,ISYMJ) 425 426 KOFF1 = IGLMRH(ISYMAL,ISYMJ) + 1 427 KOFF2 = IT2BGD(ISYMGI,ISYMJ) + 1 428 429 IF ( LWORK .LT. NBAS(ISYMAL)*NT1AO(ISYMGI) ) THEN 430 CALL QUIT('Insufficient work space in CC_21IDEN.') 431 ENDIF 432 433 NTOTAL = MAX(NBAS(ISYMAL),1) 434 NTOTGI = MAX(NT1AO(ISYMGI),1) 435 436 CALL DGEMM('N','T',NBAS(ISYMAL),NT1AO(ISYMGI),NRHF(ISYMJ), 437 * ONE,XLAMDP(KOFF1),NTOTAL,PQAO(KOFF2),NTOTGI, 438 * ZERO,WORK,NTOTAL) 439C 440C ------------------------------------------------------------ 441C resort from PQDEN(alpha, gamma i) to PQDEN(alpha i, gamma) 442C ------------------------------------------------------------ 443C 444 DO ISYMG = 1, NSYM 445 446 ISYMI = MULD2H(ISYMGI,ISYMG) 447 ISYALI = MULD2H(ISYMAL,ISYMI) 448 449 DO G = 1, NBAS(ISYMG) 450 451 KOFF4 = IT3BAO(ISYALI,ISYMG) 452 * + NT1AO(ISYALI)*(G-1) + IT1AO(ISYMAL,ISYMI) 453 454 IF (IOPT.EQ.1) THEN 455 DO I = 1, NRHF(ISYMI) 456 NGI = IT1AO(ISYMG,ISYMI) + NBAS(ISYMG)*(I-1)+G 457 NAGI = NBAS(ISYMAL)*(NGI-1) + 1 458 NAIG = KOFF4 + NBAS(ISYMAL)*(I-1) + 1 459 460 CALL DAXPY(NBAS(ISYMAL),ONE,WORK(NAGI),1, 461 * PQDEN(NAIG),1) 462 END DO 463 ELSE 464 DO I = 1, NRHF(ISYMI) 465 NGI = IT1AO(ISYMG,ISYMI) + NBAS(ISYMG)*(I-1)+G 466 NAGI = NBAS(ISYMAL)*(NGI-1) + 1 467 NAIG = KOFF4 + NBAS(ISYMAL)*(I-1) + 1 468 469 CALL DCOPY(NBAS(ISYMAL),WORK(NAGI),1, 470 * PQDEN(NAIG),1) 471 END DO 472 END IF 473 474 END DO 475 476 END DO 477 478 END DO 479 480 ELSE 481 CALL QUIT('Illegal FLAG in CC_21IDEN.') 482 END IF 483 484 RETURN 485 486 END 487*=====================================================================* 488C /* Deck cc_21icon */ 489 SUBROUTINE CC_21ICON(RHOAO,ISYRES,PQDEN,ISYDEN,XINT,ISYDIS, 490 * IT3BAO,WORK,LWORK,IOPT) 491*---------------------------------------------------------------------* 492C 493C Purpose: contract the effective density build in CC_21I 494C with the 2-el integrals to RHO1 contribution 495C 496C Written by Christof Haettig 16 July 1998 497C Modified by Sonia Coriani 07 November 1999 to handle 498C full (a b| g d) in input (IOPT = 2) 499C IOPT = 1, usual packed a>= b 500C 501*---------------------------------------------------------------------* 502#if defined (IMPLICIT_NONE) 503 IMPLICIT NONE 504#else 505# include "implicit.h" 506#endif 507#include "priunit.h" 508#include "ccorb.h" 509#include "ccsdsym.h" 510C 511 INTEGER ISYRES, ISYDEN, ISYDIS, IT3BAO(8,8), LWORK 512 INTEGER IOPT 513C 514#if defined (SYS_CRAY) 515 REAL RHOAO(*), PQDEN(*), XINT(*), WORK(LWORK) 516 REAL ZERO, ONE 517#else 518 DOUBLE PRECISION RHOAO(*), PQDEN(*), XINT(*), WORK(LWORK) 519 DOUBLE PRECISION ZERO, ONE 520#endif 521C 522 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 523C 524 INTEGER ISYMG, ISALBE, ISYMAL, ISYMBE, ISYMI, ISYALI 525 INTEGER KOFF1, KOFF3, KOFF4, KOFF5, NTOTAL, NTOTBE 526C 527C Start: check option 528C 529 IF ((IOPT.NE.1) .AND. (IOPT.NE.2)) 530 & CALL QUIT('Unknown IOPT in CC_21ICON (start)') 531C 532 DO ISYMG = 1,NSYM 533C 534 ISALBE = MULD2H(ISYDIS,ISYMG) 535C 536 IF ((IOPT.EQ.1).AND.(LWORK .LT. N2BST(ISALBE))) THEN 537 CALL QUIT('Insufficient work space in CC_21ICON.') 538 ENDIF 539C 540 DO G = 1, NBAS(ISYMG) 541 542C ---------------------------------------------------------- 543C If IOPT = 1 Square up integral distribution (al be| ga de). 544C ---------------------------------------------------------- 545 546 IF (IOPT.EQ.1) THEN 547 548 KOFF1 = IDSAOG(ISYMG,ISYDIS) + 549 & NNBST(ISALBE)*(G - 1) + 1 550 CALL CCSD_SYMSQ(XINT(KOFF1),ISALBE,WORK) 551 552 END IF 553C 554C ------------------------------------------------ 555C Final contraction and storage in result. 556C ------------------------------------------------ 557C 558 DO ISYMAL = 1, NSYM 559 560 ISYMBE = MULD2H(ISALBE,ISYMAL) 561 ISYMI = MULD2H(ISYRES,ISYMBE) 562 ISYALI = MULD2H(ISYMAL,ISYMI) 563 564 IF (MULD2H(ISYALI,ISYMG).NE.ISYDEN) THEN 565 WRITE (LUPRI,*) 'Symmetry mismatch in CC_21ICON:' 566 WRITE (LUPRI,*) ISYALI,ISYMG,ISYDEN 567 CALL QUIT('Symmetry mismatch in CC_21ICON.') 568 END IF 569 570 KOFF4 = IT3BAO(ISYALI,ISYMG) + NT1AO(ISYALI)*(G-1) 571 & + IT1AO(ISYMAL,ISYMI) + 1 572 KOFF5 = IT1AO(ISYMBE,ISYMI) + 1 573 574 NTOTAL = MAX(NBAS(ISYMAL),1) 575 NTOTBE = MAX(NBAS(ISYMBE),1) 576 577 IF (IOPT.EQ.1) THEN 578 579 KOFF3 = IAODIS(ISYMAL,ISYMBE) + 1 580 581 CALL DGEMM('T','N',NBAS(ISYMBE),NRHF(ISYMI), 582 & NBAS(ISYMAL), 583 & -ONE,WORK(KOFF3),NTOTAL,PQDEN(KOFF4),NTOTAL, 584 & ONE,RHOAO(KOFF5),NTOTBE) 585 586 ELSE IF (IOPT.EQ.2) THEN 587 588 KOFF3 = IDSAOGSQ(ISYMG,ISYDIS)+ N2BST(ISALBE)*(G - 1) 589 & + IAODIS(ISYMAL,ISYMBE) + 1 590 591 CALL DGEMM('T','N',NBAS(ISYMBE),NRHF(ISYMI), 592 & NBAS(ISYMAL), 593 & -ONE,XINT(KOFF3),NTOTAL,PQDEN(KOFF4),NTOTAL, 594 & ONE,RHOAO(KOFF5),NTOTBE) 595 ELSE 596 597 CALL QUIT('Unknown option in CC_21ICON') 598 599 END IF 600 601 END DO 602 603 END DO 604 605 END DO 606 607 RETURN 608 END 609*=====================================================================* 610