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_21i3 */ 20*-------------------------------------------------------------- 21 SUBROUTINE CC_21I3(RHO1, XINT0, ISYDIS0, XINT1, ISYDIS1, 22 * PINT0,QINT0,ISYPQ0,PINT1,QINT1,ISYPQ1, 23 * PAINT0,QAINT0,ISYPQA0, 24 * PAINT1,QAINT1,ISYPQA1, 25 * XLAMP0,XLAMH0, ISYML0,XLAMPQ,ISYMLQ, 26 * XLAMPA,ISYMLA,XLAMPQA,ISYMLQA, 27 * WORK,LWORK,IOPT,LRHOAO,LDERIV,LXI1SQ) 28*-------------------------------------------------------------- 29* 30* Generalization of the CC_21I and CC_21I2 routines 31* for the F contribution (21I contr) to the F^bT^A result vector 32* (Calculation of the FQA intermediate) 33* 34* Sonia Coriani February 1999 35* 36* Version: 2.0 37* 38* XINT0, ISYDIS0 -- (* *|* delta) batch of integrals, its symmetry 39* XINT1, ISYDIS1 -- the same for derivative integrals 40* PINT0, QINT0 -- P_ci,j (del), Q_ci,j (del) (ISYPQ0) 41* PINT1, QINT1 -- Pbar_ci,j (del), Qbar_ci,j(del) (ISYPQ1) 42* PAINT0,QAINT0 -- PA_ci,j (del) QA_ci,j(del) (ISYPQA0) 43* PAINT1,QAINT1 -- PAbar_ci,j QAbar_ci,j (del) (ISYPQA1) 44* 45* LRHOAO = .TRUE. return result with a index in AO 46* LRHOAO = .FALSE. return result with a index in MO 47* LDERIV = .FALSE. omit contribution from derivative integrals 48* LXI1SQ = .TRUE. (a b|* *) of XINT1 is a full matrix (fx LAO) 49* 50* IOPT = 1 The FA intermediate (test case) 51* IOPT = 2 The FQA intermediate 52*------------------------------------------------------------------ 53* 54#if defined (IMPLICIT_NONE) 55 IMPLICIT NONE 56#else 57# include "implicit.h" 58#endif 59#include "priunit.h" 60#include "ccorb.h" 61#include "ccsdsym.h" 62C 63 LOGICAL LRHOAO, LDERIV, LXI1SQ 64 INTEGER LWORK,ISYDIS0,ISYDIS1,IOPT 65 INTEGER ISYPQ0,ISYPQ1,ISYPQA0,ISYPQA1, ip 66 INTEGER ISYML0, ISYMLQ, ISYMLA, ISYMLQA, ISYM 67C 68#if defined (SYS_CRAY) 69 REAL RHO1(*), PINT0(*), PINT1(*), QINT0(*), QINT1(*) 70 REAL PAINT0(*), PAINT1(*), QAINT0(*), QAINT1(*) 71 REAL XINT0(*), XINT1(*), WORK(LWORK) 72 REAL XLAMP0(*), XLAMPA(*), XLAMPQ(*), XLAMPQA(*) 73 REAL XLAMH0(*) 74 REAL DUMMY, ZERO, ONE, TWO, DDOT, XNORM 75#else 76 DOUBLE PRECISION RHO1(*), PINT0(*), PINT1(*), QINT0(*), QINT1(*) 77 DOUBLE PRECISION PAINT0(*), PAINT1(*), QAINT0(*), QAINT1(*) 78 DOUBLE PRECISION XINT0(*), XINT1(*), WORK(LWORK) 79 DOUBLE PRECISION XLAMP0(*), XLAMPA(*), XLAMPQ(*), XLAMPQA(*) 80 DOUBLE PRECISION XLAMH0(*) 81 DOUBLE PRECISION DUMMY, ZERO, ONE, TWO, DDOT, XNORM 82#endif 83C 84 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 85C 86 INTEGER NT3BAO(8), IT3BAO(8,8) 87 INTEGER ISYGIJ1, ISYGIJ2, ISYDENBA, ISYDENA 88 INTEGER ISYRES, ISYMJ, ISYMGI, ISYMAL, ISY3AO, ISYMG 89 INTEGER ISYALI, ISYMI, ISALBE, ISYMBE, ISYMA, ISYGAM 90 INTEGER KRHOAO, KPQMO, KPQONE, KPQTWO, KDENBA, KDENA 91 INTEGER KAOINT,NTOTA,NPQMO, NPQMO01, NPQMOA1, NTOTGA 92 INTEGER KEND1,LWRK1,KOFF1,KOFF2,KOFF3,KOFF4,KOFF5 93 INTEGER ICOUNT, NTOTAL, NTOTGI, NGI, NAGI, NAIG, NTOTBE, NTOTAJ 94 INTEGER IOPTDEN, IOPTCON 95C 96C set some symmetries used globally in this routine: 97C 98 ISYGIJ1 = MULD2H(ISYPQ0,ISYMLA) 99* IF (IOPT.EQ.2) THEN 100* ISYGIJ2 = MULD2H(ISYPQ0,ISYMLQA) 101* END IF 102 ISYGIJ2 = MULD2H(ISYPQA0,ISYMLQ) 103 ISYDENBA = MULD2H(ISYGIJ1,ISYMLQ) 104 ISYDENA = MULD2H(ISYGIJ1,ISYML0) 105 ISYRES = MULD2H(ISYDENBA,ISYDIS0) 106C 107* IF ( IOPT.GE.2 .AND. MULD2H(ISYVWZ2,ISYXL0).NE.ISYGIJ ) THEN 108* CALL QUIT('Symmetry mismatch in CC_21I3NEW.') 109* END IF 110C 111C calculate index arrays for intermediates with 3 indices in AO: 112C 113 DO ISY3AO = 1, NSYM 114 ICOUNT = 0 115 DO ISYMG = 1, NSYM 116 ISYALI = MULD2H(ISY3AO,ISYMG) 117 IT3BAO(ISYALI,ISYMG) = ICOUNT 118 ICOUNT = ICOUNT + NBAS(ISYMG) * NT1AO(ISYALI) 119 END DO 120 NT3BAO(ISY3AO) = ICOUNT 121 END DO 122C 123C---------------------------------------------------------------- 124C allocate work space for 125C -- intermediate with one index backtransformed 126C -- intermediate with one more index backtransformed 127C -- intermediate with two more indices backtransformed 128C---------------------------------------------------------------- 129C 130 NPQMO = 1 131 DO ISYM = 1, NSYM 132 NPQMO = MAX(NPQMO,NT2BCD(ISYM)) 133 END DO 134 135C 136C IF (IOPT.GE.2) NPQMO = MAX(NPQMO,NT2BCD(ISYVWZ2)) 137C 138 KRHOAO = 1 139 KPQMO = KRHOAO + NT1AO(ISYRES) 140 KPQONE = KPQMO + NPQMO 141 IF (IOPT.EQ.2) THEN 142 KPQTWO = KPQONE + NT2BGD(ISYGIJ1) 143 ELSE 144 KPQTWO = KPQONE 145* ISYGIJ2 = ISYGIJ1 146 END IF 147 148 KDENBA = KPQTWO + NT2BGD(ISYGIJ2) 149 KEND1 = KDENBA + NT3BAO(ISYDENBA) 150 151 LWRK1 = LWORK - KEND1 152 153 IF ((IOPT.EQ.2).AND.LDERIV) THEN 154 KDENA = KEND1 155 KEND1 = KDENA + NT3BAO(ISYDENA) 156 END IF 157 LWRK1 = LWORK - KEND1 158C 159 IF ( LWRK1 .LT. 0) THEN 160 CALL QUIT('Insufficient work space in CC_21I3.') 161 ENDIF 162C 163C-------------------------------------------------------------------- 164C 1. Exchange like contribution: add P and Q intermediates 165C together and transform virtual index "c" back to AO 166C with the LambdaA^p matrix 167C add PA and QA intermediates together and backtransform 168C "c" with the Lambda^p matrix. Add to previous 169C Repeat for Pbar (P1) and Qbar (Q1) 170C Repeat for PAbar (PA1) and QAbar (QA1) 171C-------------------------------------------------------------------- 172C 173 CALL DCOPY(NT2BCD(ISYPQ0),PINT0,1,WORK(KPQMO),1) 174 CALL DAXPY(NT2BCD(ISYPQ0),ONE,QINT0,1,WORK(KPQMO),1) 175C 176 !LambdaA*PQ -> KPQONE 177 CALL CC_BFAO(WORK(KPQONE),WORK(KPQMO),ISYPQ0,XLAMPA,ISYMLA, 178 & 1,1,DUMMY,1,DUMMY,1,ZERO,0) 179c 180c XNORM = DDOT(NT2BGD(ISYGIJ1),WORK(KPQONE),1,WORK(KPQONE),1) 181c WRITE (LUPRI,*) 'CC_21I3> Norm(PQ-AO)1=',XNORM 182c 183 IF (IOPT.EQ.2) THEN 184 !LambdaQA*PQ -> KPQTWO 185 CALL CC_BFAO(WORK(KPQTWO),WORK(KPQMO),ISYPQ0,XLAMPQA,ISYMLQA, 186 & 1,1,DUMMY,1,DUMMY,1,ZERO,0) 187 END IF 188 189 CALL DCOPY(NT2BCD(ISYPQA0),PAINT0,1,WORK(KPQMO),1) 190c WRITE(LUPRI,*) 'CC_21I3> test PQ' 191c do ip = 1, 5 192c WRITE(LUPRI,*) PAINT0(ip) 193c end do 194 CALL DAXPY(NT2BCD(ISYPQA0),ONE,QAINT0,1,WORK(KPQMO),1) 195 196 !Lambda0*calPQ -> + KPQONE 197 CALL CC_BFAO(WORK(KPQONE),WORK(KPQMO),ISYPQA0,XLAMP0,ISYML0, 198 & 1,1,DUMMY,1,DUMMY,1,ONE,0) 199c 200c XNORM = DDOT(NT2BGD(ISYGIJ1),WORK(KPQONE),1,WORK(KPQONE),1) 201c WRITE (LUPRI,*) 'CC_21I3> Norm(PQA-AO)1=',XNORM 202c 203 204 IF (IOPT.EQ.2) THEN 205 !LambdaQ*calPQ -> + KPQTWO 206 CALL CC_BFAO(WORK(KPQTWO),WORK(KPQMO),ISYPQA0,XLAMPQ,ISYMLQ, 207 & 1,1,DUMMY,1,DUMMY,1,ONE,0) 208C-------------------------------------------------------------------- 209C Backtransform also barred PQ and calPQ and add to KPQTWO 210C-------------------------------------------------------------------- 211 212 CALL DCOPY(NT2BCD(ISYPQ1),PINT1,1,WORK(KPQMO),1) 213 CALL DAXPY(NT2BCD(ISYPQ1),ONE,QINT1,1,WORK(KPQMO),1) 214C 215 !LambdaA*barPQ -> + KPQTWO 216 CALL CC_BFAO(WORK(KPQTWO),WORK(KPQMO),ISYPQ1,XLAMPA,ISYMLA, 217 & 1,1,DUMMY,1,DUMMY,1,ONE,0) 218 219 CALL DCOPY(NT2BCD(ISYPQA1),PAINT1,1,WORK(KPQMO),1) 220 !LambdaA*bar-calPQ -> + KPQTWO 221 CALL DAXPY(NT2BCD(ISYPQA1),ONE,QAINT1,1,WORK(KPQMO),1) 222 223 CALL CC_BFAO(WORK(KPQTWO),WORK(KPQMO),ISYPQA1,XLAMP0,ISYML0, 224 & 1,1,DUMMY,1,DUMMY,1,ONE,0) 225 END IF 226C 227C-------------------------------------------------------------------- 228C transform occupied index j back to AO (alpha) 229C-------------------------------------------------------------------- 230C 231 IOPTDEN = 0 !initialize density with actual contribution 232 !backt j with XLAMPQ 233 CALL CC_21IDEN(WORK(KPQONE),ISYGIJ1,XLAMPQ,ISYMLQ,WORK(KDENBA), 234 * ISYDENBA,IT3BAO,WORK(KEND1),LWRK1,'EXC',IOPTDEN) 235 236 IF (IOPT.EQ.2) THEN 237 IOPTDEN = 1 !add to density 238 !backt j with XLAMP0 239 CALL CC_21IDEN(WORK(KPQTWO),ISYGIJ2,XLAMP0,ISYML0, 240 * WORK(KDENBA),ISYDENBA,IT3BAO, 241 * WORK(KEND1),LWRK1,'EXC',IOPTDEN) 242 243 IF (LDERIV) THEN 244 IOPTDEN = 0 245 CALL CC_21IDEN(WORK(KPQONE),ISYGIJ1,XLAMP0,ISYML0, 246 * WORK(KDENA),ISYDENA,IT3BAO, 247 * WORK(KEND1),LWRK1,'EXC',IOPTDEN) 248 END IF 249 END IF 250C 251c XNORM = DDOT(NT3BAO(ISYGIJ1),WORK(KDENBA),1,WORK(KPQDEN),1) 252c WRITE (LUPRI,*) 'Norm(PQDEN)1=',XNORM 253C 254C-------------------------------------------------------------------- 255C 2. Coulomb like contribution: scale P intermediate with -2 256C and transform virtual "c" back to AO 257C-------------------------------------------------------------------- 258C 259 CALL DCOPY(NT2BCD(ISYPQ0),PINT0,1,WORK(KPQMO),1) 260 CALL DSCAL(NT2BCD(ISYPQ0),-TWO,WORK(KPQMO),1) 261 262 CALL CC_BFAO(WORK(KPQONE),WORK(KPQMO),ISYPQ0,XLAMPA,ISYMLA, 263 & 1,1,DUMMY,1,DUMMY,1,ZERO,0) 264C 265c XNORM = DDOT(NT2BGD(ISYGIJ),WORK(KPQAO),1,WORK(KPQAO),1) 266c WRITE (LUPRI,*) 'Norm(PQAO)3=',XNORM 267C 268C 269 IF (IOPT.EQ.2) THEN 270 CALL CC_BFAO(WORK(KPQTWO),WORK(KPQMO),ISYPQ0,XLAMPQA, 271 & ISYMLQA,1,1,DUMMY,1,DUMMY,1,ZERO,0) 272 END IF 273C 274 275 CALL DCOPY(NT2BCD(ISYPQA0),PAINT0,1,WORK(KPQMO),1) 276 CALL DSCAL(NT2BCD(ISYPQA0),-TWO,WORK(KPQMO),1) 277 278 CALL CC_BFAO(WORK(KPQONE),WORK(KPQMO),ISYPQA0,XLAMP0,ISYML0, 279 & 1,1,DUMMY,1,DUMMY,1,ONE,0) 280 281 IF (IOPT.EQ.2) THEN 282 CALL CC_BFAO(WORK(KPQTWO),WORK(KPQMO),ISYPQA0,XLAMPQ,ISYMLQ, 283 & 1,1,DUMMY,1,DUMMY,1,ONE,0) 284 285C 286c XNORM = DDOT(NT2BGD(ISYGIJ),WORK(KPQAO),1,WORK(KPQAO),1) 287c WRITE (LUPRI,*) 'Norm(PQAO)4=',XNORM 288C 289 290 CALL DCOPY(NT2BCD(ISYPQ1),PINT1,1,WORK(KPQMO),1) 291 CALL DSCAL(NT2BCD(ISYPQ1),-TWO,WORK(KPQMO),1) 292 293 CALL CC_BFAO(WORK(KPQTWO),WORK(KPQMO),ISYPQ1,XLAMPA,ISYMLA, 294 & 1,1,DUMMY,1,DUMMY,1,ONE,0) 295 296 CALL DCOPY(NT2BCD(ISYPQA1),PAINT1,1,WORK(KPQMO),1) 297 CALL DSCAL(NT2BCD(ISYPQA1),-TWO,WORK(KPQMO),1) 298 299 CALL CC_BFAO(WORK(KPQTWO),WORK(KPQMO),ISYPQA1,XLAMP0,ISYML0, 300 & 1,1,DUMMY,1,DUMMY,1,ONE,0) 301 302 END IF 303C-------------------------------------------------------------------- 304C transform occupied index j back to AO (gamma), add result to 305C the effective density stored in KPQDEN 306C-------------------------------------------------------------------- 307C 308 309 IOPTDEN = 1 310 CALL CC_21IDEN(WORK(KPQONE),ISYGIJ1,XLAMPQ,ISYMLQ,WORK(KDENBA), 311 * ISYDENBA,IT3BAO,WORK(KEND1),LWRK1,'COU',IOPTDEN) 312 313 IF (IOPT.EQ.2) THEN 314 IOPTDEN = 1 315 CALL CC_21IDEN(WORK(KPQTWO),ISYGIJ2,XLAMP0,ISYML0, 316 * WORK(KDENBA),ISYDENBA,IT3BAO, 317 * WORK(KEND1),LWRK1,'COU',IOPTDEN) 318 319 IF (LDERIV) THEN 320 IOPTDEN = 1 321 CALL CC_21IDEN(WORK(KPQONE),ISYGIJ1,XLAMP0,ISYML0, 322 * WORK(KDENA),ISYDENA,IT3BAO, 323 * WORK(KEND1),LWRK1,'COU',IOPTDEN) 324 END IF 325 END IF 326C 327C-------------------------------------------------------------------- 328C contract the effective density with the 2-el integrals: 329C-------------------------------------------------------------------- 330C 331c XNORM = DDOT(NT3BAO(ISYGIJ),WORK(KPQDEN),1,WORK(KPQDEN),1) 332c WRITE (LUPRI,*) 'Norm(PQDEN)2=',XNORM 333C 334 CALL DZERO(WORK(KRHOAO),NT1AO(ISYRES)) 335 336 CALL CC_21ICON(WORK(KRHOAO),ISYRES,WORK(KDENBA),ISYDENBA, 337 * XINT0,ISYDIS0,IT3BAO,WORK(KEND1),LWRK1,1) 338 339 IF ((IOPT.EQ.2).AND.LDERIV) THEN 340 341 IF (LXI1SQ) THEN 342 IOPTCON = 2 !X2INT has full (a b| already 343 ELSE 344 IOPTCON = 1 !X2INT is squared inside CC_21ICON1 345 END IF 346 347 CALL CC_21ICON(WORK(KRHOAO),ISYRES,WORK(KDENA),ISYDENA, 348 * XINT1,ISYDIS1,IT3BAO,WORK(KEND1),LWRK1,IOPTCON) 349 END IF 350C 351c XNORM = DDOT(NT1AO(ISYRES),RHO1,1,RHO1,1) 352c WRITE (LUPRI,*) 'XNORM=',XNORM 353c 354C--------------------------------------------- 355C transformation and/or storage in result: 356C--------------------------------------------- 357C 358 IF (LRHOAO) THEN 359 CALL DAXPY(NT1AO(ISYRES),ONE,WORK(KRHOAO),1,RHO1,1) 360 ELSE 361 CALL CC_T1AM(RHO1,ISYRES,WORK(KRHOAO),ISYRES,XLAMH0,ISYML0,ONE) 362 END IF 363C 364 RETURN 365 END 366C 367*=====================================================================* 368