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*=====================================================================* 20 SUBROUTINE CCSDT_FMAT_NODDY(LISTL,IDLSTL,LISTB,IDLSTB, 21 & IOPTRES,NEW_NODDY_F_CONT, 22 & OMEGA1,OMEGA2, 23 & OMEGA1EFF,OMEGA2EFF, 24 & IDOTS,DOTPROD,LISTDP,ITRAN, 25 & NFTRAN,MXVEC,WORK,LWORK) 26*---------------------------------------------------------------------* 27* 28* Purpose: compute triples contribution to F matrix transformation 29* 30* (F T^B)^eff_1,2 = (F T^B)_1,2(CCSD) 31* + (F T^B)_1,2(L3,T^B3) 32* - (F T^B)_3 A_3;1,2 (w_3 - w)^1 33* 34* 35* Written by Christof Haettig, April 2002 36* based on different other noddy codes 37* 38*=====================================================================* 39 IMPLICIT NONE 40#include "priunit.h" 41#include "ccsdinp.h" 42#include "maxorb.h" 43#include "ccsdsym.h" 44#include "ccfield.h" 45#include "ccorb.h" 46#include "inftap.h" 47#include "dummy.h" 48#include "ccnoddy.h" 49 50 LOGICAL LOCDBG, FD_TEST, XI_ONLY 51 PARAMETER (LOCDBG=.FALSE., FD_TEST=.FALSE., XI_ONLY=.FALSE.) 52 53 INTEGER ISYM0 54 PARAMETER (ISYM0 = 1) 55 56 LOGICAL NEW_NODDY_F_CONT 57 CHARACTER*3 LISTL, LISTB, LISTDP 58 INTEGER LWORK, IDLSTL, IDLSTB, ITRAN, NFTRAN, MXVEC 59 INTEGER IDOTS(MXVEC,NFTRAN), IOPTRES 60#if defined (SYS_CRAY) 61 REAL DOTPROD(MXVEC,NFTRAN), DDOT 62 REAL WORK(LWORK), FREQB, FREQL, FREQF, FREQC 63 REAL OMEGA1(NT1AMX), OMEGA2(NT2AMX) 64 REAL OMEGA1EFF(NT1AMX), OMEGA2EFF(NT2AMX) 65 REAL SIXTH, ONE, TWO, TCON, SCON, DCON, FF, SIGN 66#else 67 DOUBLE PRECISION DOTPROD(MXVEC,NFTRAN), DDOT 68 DOUBLE PRECISION WORK(LWORK), FREQB, FREQL, FREQF, FREQC 69 DOUBLE PRECISION OMEGA1(NT1AMX), OMEGA2(NT2AMX) 70 DOUBLE PRECISION OMEGA1EFF(NT1AMX), OMEGA2EFF(NT2AMX) 71 DOUBLE PRECISION SIXTH, ONE, TWO, TCON, SCON, DCON, FF, SIGN 72#endif 73 PARAMETER(SIXTH=1.0D0/6.0D0, ONE=1.0D0, TWO=2.0D0) 74 75 CHARACTER*10 MODEL 76 LOGICAL SKIP_T3_CONT, L2INCL, SKIP_F3 77 INTEGER KXINT, KXIAJB, KYIAJB, KT1AMP0, KLAMP0, KLAMH0, KEND1, 78 & LWRK1, KINT1SB, KINT2SB, KEND2, LWRK2, KTB3AM, KL1AM, 79 & KL2AM, KINT1T0, KINT2T0, KINT1TB, KINT2TB, KEND3, LWRK3, 80 & K0IOVVO, K0IOOVV, K0IOOOO, K0IVVVV, 81 & KBIOVVO, KBIOOVV, KBIOOOO, KBIVVVV, 82 & KL3AM, KT02AM, KF3AM, IRECNR, KINT1S0, KINT2S0, KSCR1, 83 & KFOCK0, KFOCKB, KLAMPB, KLAMHB, KOME1, KOME2, KDUM, 84 & LUFOCK, ISYML, KFOCKD, ISYMB, KEND4, LWRK4, 85 & KTC3AM, KINT1SC, KINT2SC, KLAMPC, KLAMHC, KFOCKC, IDLSTC, 86 & IVEC, ISYMC, KFFD1AM, KFFD2AM, KFFD3AM, KTB1AM, KTB2AM, 87 & ISYMD, ILLL, IDEL, ISYDIS, IJ, NIJ, INDEX, LUTEMP, IDX, 88 & IOPT, ILSTSYM, KTC1AM, KTC2AM, KFIELD, KFLDB1, KFIELDAO, 89 & KT03AM, KT3SCR 90 91 INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 92 93 CALL QENTER('CCSDT_FMAT_NODDY') 94 95 IF (DIRECT) CALL QUIT('CCSDT_FMAT_NODDY: DIRECT NOT IMPLEMENTED') 96 97 IF (LOCDBG) THEN 98 WRITE(LUPRI,*) 'CCSDT_FMAT_NODDY> F matrix result on entry:' 99 CALL CC_PRP(OMEGA1,OMEGA2,1,1,1) 100 END IF 101 102 SKIP_T3_CONT = .FALSE. 103 IF (NEW_NODDY_F_CONT) THEN 104 SKIP_T3_CONT = (IOPTRES.EQ.5) 105 IF (LOCDBG .OR. DEBUG) THEN 106 WRITE(LUPRI,*) 'Use new noddy F matrix implementation...' 107 IF (SKIP_T3_CONT) 108 & WRITE(LUPRI,*) 'Contributions from T3 will be skipped here...' 109 END IF 110 ENDIF 111*---------------------------------------------------------------------* 112* Memory allocation: 113*---------------------------------------------------------------------* 114 KSCR1 = 1 115 KFOCKD = KSCR1 + NT1AMX 116 KEND1 = KFOCKD + NORBT 117 118 KFOCK0 = KEND1 119 KEND1 = KFOCK0 + NORBT*NORBT 120 121 IF (NONHF) THEN 122 KFIELD = KEND1 123 KFLDB1 = KFIELD + NORBT*NORBT 124 KFIELDAO = KFLDB1 + NORBT*NORBT 125 KEND1 = KFIELDAO + NORBT*NORBT 126 END IF 127 128 KT1AMP0 = KEND1 129 KLAMP0 = KT1AMP0 + NT1AMX 130 KLAMH0 = KLAMP0 + NLAMDT 131 KEND1 = KLAMH0 + NLAMDT 132 133 KTB1AM = KEND1 134 KFOCKB = KTB1AM + NT1AMX 135 KLAMPB = KFOCKB + NORBT*NORBT 136 KLAMHB = KLAMPB + NLAMDT 137 KEND1 = KLAMHB + NLAMDT 138 139 KOME1 = KEND1 140 KOME2 = KOME1 + NT1AMX 141 KEND1 = KOME2 + NT1AMX*NT1AMX 142 143 KL1AM = KEND1 144 KL2AM = KL1AM + NT1AMX 145 KEND1 = KL2AM + NT1AMX*NT1AMX 146 147 KINT1T0 = KEND1 148 KINT2T0 = KINT1T0 + NT1AMX*NVIRT*NVIRT 149 KEND1 = KINT2T0 + NRHFT*NRHFT*NT1AMX 150 151 KINT1S0 = KEND1 152 KINT2S0 = KINT1S0 + NT1AMX*NVIRT*NVIRT 153 KEND1 = KINT2S0 + NRHFT*NRHFT*NT1AMX 154 155 KXIAJB = KEND1 156 KYIAJB = KXIAJB + NT1AMX*NT1AMX 157 KEND1 = KYIAJB + NT1AMX*NT1AMX 158 159 KINT1SB = KEND1 160 KINT2SB = KINT1SB + NT1AMX*NVIRT*NVIRT 161 KEND1 = KINT2SB + NRHFT*NRHFT*NT1AMX 162 163 LWRK1 = LWORK - KEND1 164 IF (LWRK1 .LT. 0) THEN 165 write(lupri,*) 'Need:',kend1,' have:',lwork 166 CALL QUIT('Insufficient space in CCSDT_FMAT_NODDY (1)') 167 ENDIF 168 169 CALL DZERO(WORK(KOME1),NT1AMX) 170 CALL DZERO(WORK(KOME2),NT1AMX*NT1AMX) 171 172*---------------------------------------------------------------------* 173* Get zeroth-order Lambda matrices: 174*---------------------------------------------------------------------* 175 IOPT = 1 176 KDUM = KEND1 177 Call CC_RDRSP('R0',0,ISYM0,IOPT,MODEL,WORK(KT1AMP0),WORK(KDUM)) 178 179 Call LAMMAT(WORK(KLAMP0),WORK(KLAMH0),WORK(KT1AMP0), 180 & WORK(KEND1),LWRK1) 181 182*---------------------------------------------------------------------* 183* Read precalculated integrals from file: 184* XINT1S0 = (CK|BD) 185* XINT2S0 = (CK|LJ) 186* XINT1T0 = (KC|BD) 187* XINT2T0 = (KC|LJ) 188* XIAJB = 2(IA|JB) - (IB|JA) 189* YIAJB = (IA|JB) 190*---------------------------------------------------------------------* 191 CALL CCSDT_READ_NODDY(.TRUE.,WORK(KFOCKD),WORK(KFOCK0), 192 & WORK(KFIELD),WORK(KFIELDAO), 193 & .TRUE.,WORK(KXIAJB),WORK(KYIAJB), 194 & .TRUE.,WORK(KINT1S0),WORK(KINT2S0), 195 & .TRUE.,WORK(KINT1T0),WORK(KINT2T0), 196 & .FALSE.,DUMMY,DUMMY,DUMMY,DUMMY, 197 & NORBT,NLAMDT,NRHFT,NVIRT,NT1AMX) 198 199*---------------------------------------------------------------------* 200* If needed, compute zeroth-order triples amplitudes: 201*---------------------------------------------------------------------* 202 IF (NONHF) THEN 203 KT03AM = KEND1 204 KEND2 = KT03AM + NT1AMX*NT1AMX*NT1AMX 205 206 LWRK2 = LWORK - KEND2 207 IF (LWRK2 .LT. 0) THEN 208 CALL QUIT('Insufficient space in CCSDT_FMAT_NODDY (T03)') 209 END IF 210 211 LUTEMP = -1 212 CALL GPOPEN(LUTEMP,FILNODT30,'UNKNOWN',' ','UNFORMATTED', 213 & IDUMMY,.FALSE.) 214 READ(LUTEMP) (WORK(KT03AM+I-1), I=1,NT1AMX*NT1AMX*NT1AMX) 215 CALL GPCLOSE(LUTEMP,'KEEP') 216 217 218 ! store triples amplitudes on file 219 LUTEMP = -1 220 CALL GPOPEN(LUTEMP,'T3AMPL','UNKNOWN',' ','UNFORMATTED', 221 & IDUMMY,.FALSE.) 222 REWIND LUTEMP 223 WRITE (LUTEMP) (WORK(KT03AM-1+IDX),IDX=1,NT1AMX*NT1AMX*NT1AMX) 224 CALL GPCLOSE(LUTEMP,'KEEP') 225 END IF 226 227*---------------------------------------------------------------------* 228* Compute T^B_3 amplitudes and some more integrals: 229* 230* LISTB = "R1" --> first-order response amplitudes 231* = "RE" --> right eigenvectors 232* 233* XINT1SB = (CK|BD)-bar 234* XINT2SB = (CK|LJ)-bar 235*---------------------------------------------------------------------* 236 KTB3AM = KEND1 237 KEND2 = KTB3AM + NT1AMX*NT1AMX*NT1AMX 238 239 ! for finite difference we need T3 until the end 240 IF (FD_TEST) KEND1 = KEND2 241 242 LWRK2 = LWORK - KEND2 243 IF (LWRK2 .LT. 0) THEN 244 CALL QUIT('Insufficient space in CCSDT_FMAT_NODDY (2)') 245 ENDIF 246 247 IF (LISTB(1:3).EQ.'R1 ' .OR. LISTB(1:3).EQ.'RE ' .OR. 248 & LISTB(1:3).EQ.'RC ' ) THEN 249 KDUM = KEND2 250 CALL CCSDT_T31_NODDY(WORK(KTB3AM),LISTB,IDLSTB,FREQB,XI_ONLY, 251 & .FALSE.,WORK(KINT1S0),WORK(KINT2S0), 252 & .FALSE.,WORK(KINT1T0),WORK(KINT2T0), 253 & .FALSE.,WORK(KXIAJB),WORK(KYIAJB), 254 & WORK(KINT1SB),WORK(KINT2SB), 255 & WORK(KLAMPB),WORK(KLAMHB),WORK(KFOCKB), 256 & WORK(KLAMP0),WORK(KLAMH0),WORK(KFOCK0), 257 & WORK(KDUM),WORK(KFOCKD), 258 & WORK(KEND2),LWRK2) 259 ELSE IF (LISTB(1:3).EQ.'R2 ' .OR. LISTB(1:3).EQ.'ER1') THEN 260 CALL CCSDT_T32_NODDY(WORK(KTB3AM),LISTB,IDLSTB,FREQB, 261 & WORK(KINT1S0),WORK(KINT2S0), 262 & WORK(KLAMP0),WORK(KLAMH0),WORK(KFOCK0), 263 & WORK(KFOCKD),WORK(KFIELDAO),WORK(KFIELD), 264 & WORK(KSCR1),WORK(KEND2),LWRK2) 265 IOPT = 1 266 ISYMB = ILSTSYM(LISTB,IDLSTB) 267 CALL CC_RDRSP(LISTB,IDLSTB,ISYMB,IOPT,MODEL,WORK(KTB1AM),DUMMY) 268 CALL CCLR_LAMTRA(WORK(KLAMP0),WORK(KLAMPB),WORK(KLAMH0), 269 & WORK(KLAMHB),WORK(KTB1AM),ISYMB) 270 CALL CCSDT_INTS1_NODDY(.TRUE.,WORK(KINT1SB),WORK(KINT2SB), 271 & .FALSE.,DUMMY,DUMMY, 272 & WORK(KLAMP0),WORK(KLAMH0), 273 & WORK(KLAMPB),WORK(KLAMHB), 274 & WORK(KEND2),LWRK2) 275 ELSE 276 CALL QUIT('Unknown or illegal list in CCSDT_FMAT_NODDY.') 277 END IF 278 279 IF (NONHF .OR. FD_TEST) THEN 280 LUTEMP = -1 281 CALL GPOPEN(LUTEMP,'T3BARAM','UNKNOWN',' ','UNFORMATTED', 282 & IDUMMY,.FALSE.) 283 REWIND LUTEMP 284 WRITE (LUTEMP) (WORK(KTB3AM-1+IDX),IDX=1,NT1AMX*NT1AMX*NT1AMX) 285 CALL GPCLOSE(LUTEMP,'KEEP') 286 IF (LOCDBG) WRITE(LUPRI,*) 'Dumped R1_3... norm^2=', 287 & DDOT(NT1AMX*NT1AMX*NT1AMX,WORK(KTB3AM),1,WORK(KTB3AM),1) 288 END IF 289 290*---------------------------------------------------------------------* 291* Compute contribution from <L_2|[[H,T^B_3],\tau_nu_1|HF>: 292*---------------------------------------------------------------------* 293 IF (LWRK2 .LT. NT2AMX) THEN 294 CALL QUIT('Insufficient space in CCSDT_FMAT_NODDY (3)') 295 ENDIF 296 297 ISYML = ILSTSYM(LISTL,IDLSTL) 298 299 ! read left amplitudes from file and square up the doubles part 300 IOPT = 3 301 Call CC_RDRSP(LISTL,IDLSTL,ISYML,IOPT,MODEL, 302 & WORK(KL1AM),WORK(KEND2)) 303 CALL CC_T2SQ(WORK(KEND2),WORK(KL2AM),ISYML) 304 305 IF (.NOT.SKIP_T3_CONT) THEN ! might be done via CCSDT_FBC_NODDY 306 CALL T3_LEFT1(WORK(KOME1),WORK(KYIAJB),WORK(KXIAJB), 307 & WORK(KL2AM),WORK(KTB3AM)) 308 309 CALL T3_LEFT2(WORK(KOME1),WORK(KYIAJB),WORK(KXIAJB), 310 & WORK(KL2AM),WORK(KTB3AM)) 311 312 CALL T3_LEFT3(WORK(KOME1),WORK(KXIAJB),WORK(KL2AM),WORK(KTB3AM)) 313 END IF 314 315*---------------------------------------------------------------------* 316* Compute integrals needed for the following contributions: 317*---------------------------------------------------------------------* 318 KEND2 = KEND1 319 320 K0IOVVO = KEND2 321 K0IOOVV = K0IOVVO + NRHFT*NVIRT*NVIRT*NRHFT 322 K0IOOOO = K0IOOVV + NRHFT*NVIRT*NVIRT*NRHFT 323 K0IVVVV = K0IOOOO + NRHFT*NRHFT*NRHFT*NRHFT 324 KEND2 = K0IVVVV + NVIRT*NVIRT*NVIRT*NVIRT 325 326 KBIOVVO = KEND2 327 KBIOOVV = KBIOVVO + NRHFT*NVIRT*NVIRT*NRHFT 328 KBIOOOO = KBIOOVV + NRHFT*NVIRT*NVIRT*NRHFT 329 KBIVVVV = KBIOOOO + NRHFT*NRHFT*NRHFT*NRHFT 330 KEND2 = KBIVVVV + NVIRT*NVIRT*NVIRT*NVIRT 331 332 KINT1TB = KEND2 333 KINT2TB = KINT1TB + NT1AMX*NVIRT*NVIRT 334 KEND2 = KINT2TB + NRHFT*NRHFT*NT1AMX 335 336 LWRK2 = LWORK - KEND2 337 IF (LWRK2 .LT. 0) THEN 338 WRITE(LUPRI,*) 'Need : ',KEND2,'Available : ',LWORK 339 CALL QUIT('Insufficient space in CCSDT_FMAT_NODDY (4)') 340 ENDIF 341 342 CALL CCSDT_READ_NODDY(.FALSE.,WORK(KFOCKD),WORK(KFOCK0), 343 & WORK(KFIELD),WORK(KFIELDAO), 344 & .FALSE.,WORK(KXIAJB),WORK(KYIAJB), 345 & .FALSE.,WORK(KINT1S0),WORK(KINT2S0), 346 & .FALSE.,WORK(KINT1T0),WORK(KINT2T0), 347 & .TRUE.,WORK(K0IOVVO),WORK(K0IOOVV), 348 & WORK(K0IOOOO),WORK(K0IVVVV), 349 & NORBT,NLAMDT,NRHFT,NVIRT,NT1AMX) 350 351 CALL DZERO(WORK(KBIOVVO),NRHFT*NVIRT*NVIRT*NRHFT) 352 CALL DZERO(WORK(KBIOOVV),NRHFT*NVIRT*NVIRT*NRHFT) 353 CALL DZERO(WORK(KBIOOOO),NRHFT*NRHFT*NRHFT*NRHFT) 354 CALL DZERO(WORK(KBIVVVV),NVIRT*NVIRT*NVIRT*NVIRT) 355 356 CALL DZERO(WORK(KINT1TB),NT1AMX*NVIRT*NVIRT) 357 CALL DZERO(WORK(KINT2TB),NRHFT*NRHFT*NT1AMX) 358 359 DO ISYMD = 1, NSYM 360 DO ILLL = 1,NBAS(ISYMD) 361 IDEL = IBAS(ISYMD) + ILLL 362 ISYDIS = MULD2H(ISYMD,ISYMOP) 363 364C ---------------------------- 365C Work space allocation no. 2. 366C ---------------------------- 367 KXINT = KEND2 368 KEND3 = KXINT + NDISAO(ISYDIS) 369 LWRK3 = LWORK - KEND3 370 IF (LWRK3 .LT. 0) THEN 371 WRITE(LUPRI,*) 'Need : ',KEND3,'Available : ',LWORK 372 CALL QUIT('Insufficient space in CCSDT_FMAT_NODDY (5)') 373 ENDIF 374 375C --------------------------- 376C Read in batch of integrals. 377C --------------------------- 378 CALL CCRDAO(WORK(KXINT),IDEL,1,WORK(KEND3),LWRK3, 379 * IRECNR,DIRECT) 380 381C ---------------------------------- 382C Calculate integrals needed in CC3: 383C ---------------------------------- 384 385 CALL CCSDT_TRAN1_R(WORK(KINT1TB),WORK(KINT2TB), 386 & WORK(KLAMP0),WORK(KLAMH0), 387 & WORK(KLAMPB),WORK(KLAMHB), 388 & WORK(KXINT),IDEL) 389 390 CALL CCFOP_TRAN1_R(WORK(KBIOVVO),WORK(KBIOOVV), 391 & WORK(KBIOOOO),WORK(KBIVVVV), 392 & WORK(KLAMP0),WORK(KLAMH0), 393 & WORK(KLAMPB),WORK(KLAMHB), 394 & WORK(KLAMP0),WORK(KLAMH0), 395 & WORK(KXINT),IDEL) 396 397 CALL CCFOP_TRAN1_R(WORK(KBIOVVO),WORK(KBIOOVV), 398 & WORK(KBIOOOO),WORK(KBIVVVV), 399 & WORK(KLAMP0),WORK(KLAMH0), 400 & WORK(KLAMP0),WORK(KLAMH0), 401 & WORK(KLAMPB),WORK(KLAMHB), 402 & WORK(KXINT),IDEL) 403 404 END DO 405 END DO 406 407*---------------------------------------------------------------------* 408* Compute L_3 multipliers: 409*---------------------------------------------------------------------* 410 KEND3 = KEND2 411 412 KL3AM = KEND3 413 KEND3 = KL3AM + NT1AMX*NT1AMX*NT1AMX 414 415 LWRK3 = LWORK - KEND3 416 IF (LWRK3 .LT. 0) THEN 417 CALL QUIT('Insufficient space in CCSDT_FMAT_NODDY (6)') 418 ENDIF 419 420 421 IF (LISTL(1:3).EQ.'L0 ') THEN 422 423 FREQL = 0.0D0 424 425 IF (NONHF .AND. LWRK3.LT.NT1AMX*NT1AMX*NT1AMX) 426 * CALL QUIT('Out of memory in CCSDT_F_NODDY.') 427 428 CALL CCSDT_L03AM(WORK(KL3AM),WORK(KINT1T0),WORK(KINT2T0), 429 * WORK(KXIAJB),WORK(KFOCK0),WORK(KL1AM), 430 * WORK(KL2AM),WORK(KSCR1),WORK(KFOCKD), 431 * WORK(KFIELD),WORK(KEND3)) 432 433 CALL DSCAL(NT1AMX*NT1AMX*NT1AMX,-1.0D0,WORK(KL3AM),1) 434 435 ELSE IF (LISTL(1:3).EQ.'L1 ' .OR. LISTL(1:3).EQ.'LE ' .OR. 436 & LISTL(1:3).EQ.'M1 ' .OR. LISTL(1:3).EQ.'N2 ' .OR. 437 & LISTL(1:3).EQ.'E0 ' ) THEN 438 439 CALL CCSDT_TBAR31_NODDY(WORK(KL3AM),FREQL,LISTL,IDLSTL, 440 & WORK(KLAMP0),WORK(KLAMH0), 441 & WORK(KFOCK0),WORK(KFOCKD),WORK(KSCR1), 442 & WORK(KXIAJB),WORK(KINT1T0),WORK(KINT2T0), 443 & WORK(KEND3),LWRK3) 444 445 CALL DSCAL(NT1AMX*NT1AMX*NT1AMX,-1.0D0,WORK(KL3AM),1) 446 447 ELSE 448 449 ! FREQL = ?? 450 451 CALL QUIT('CCSDT_FMAT_NODDY> LISTL NOT AVAILABLE:'//LISTL) 452 453 END IF 454 455*---------------------------------------------------------------------* 456* Compute contribution from <L_3|[[H^B,T^0_2],\tau_nu_1]|HF> 457* and <L_3|[[H^0,T^B_2],\tau_nu_1]|HF> 458* for finite difference include "L_3|[V,T^B_3],\tau_nu_1]|HF> 459*---------------------------------------------------------------------* 460 KT02AM = KEND3 461 KTB2AM = KT02AM + NT1AMX*NT1AMX 462 KEND3 = KTB2AM + NT1AMX*NT1AMX 463 LWRK3 = LWORK - KEND3 464 IF (LWRK3 .LT. NT2AMX) THEN 465 CALL QUIT('Insufficient space in CCSDT_FMAT_NODDY (7)') 466 ENDIF 467 468 ! read T^0 doubles amplitudes from file and square up 469 IOPT = 2 470 KDUM = KEND3 471 Call CC_RDRSP('R0',0,ISYM0,IOPT,MODEL,WORK(KDUM),WORK(KEND3)) 472 CALL CC_T2SQ(WORK(KEND3),WORK(KT02AM),ISYM0) 473 474 CALL CC3_L3_OMEGA1_NODDY(WORK(KOME1),WORK(KL3AM), 475 & WORK(KBIOOOO),WORK(KBIOVVO), 476 & WORK(KBIOOVV),WORK(KBIVVVV), 477 & WORK(KT02AM)) 478 479 480 481 ! read T^B amplitudes from file and square doubles up 482 ISYMB = ILSTSYM(LISTB,IDLSTB) 483 IOPT = 3 484 Call CC_RDRSP(LISTB,IDLSTB,ISYMB,IOPT,MODEL, 485 & WORK(KTB1AM),WORK(KEND3)) 486 Call CCLR_DIASCL(WORK(KEND3),TWO,ISYMB) 487 CALL CC_T2SQ(WORK(KEND3),WORK(KTB2AM),ISYMB) 488 489 CALL CC3_L3_OMEGA1_NODDY(WORK(KOME1),WORK(KL3AM), 490 & WORK(K0IOOOO),WORK(K0IOVVO), 491 & WORK(K0IOOVV),WORK(K0IVVVV), 492 & WORK(KTB2AM)) 493 494 495 IF (NONHF) THEN 496 KTB3AM = KEND3 497 KEND4 = KTB3AM + NT1AMX*NT1AMX*NT1AMX 498 LWRK4 = LWORK - KEND4 499 IF (LWRK4 .LT. 0) THEN 500 CALL QUIT('Insufficient space in CCSDT_FMAT_NODDY (nhf1)') 501 ENDIF 502 503 LUTEMP = -1 504 CALL GPOPEN(LUTEMP,'T3BARAM','UNKNOWN',' ','UNFORMATTED', 505 & IDUMMY,.FALSE.) 506 REWIND LUTEMP 507 READ (LUTEMP) (WORK(KTB3AM-1+IDX),IDX=1,NT1AMX*NT1AMX*NT1AMX) 508 CALL GPCLOSE(LUTEMP,'KEEP') 509 510 CALL CCSDT_E1AM(WORK(KOME1),WORK(KL3AM), 511 & WORK(KTB3AM),WORK(KFIELD)) 512 END IF 513 514 515 DO I = 1,NT1AMX 516 OMEGA1(I) = OMEGA1(I) + WORK(KOME1+I-1) 517 END DO 518 519*---------------------------------------------------------------------* 520* Compute contribution from <L_3|[H^B,\tau_nu_2]|HF> 521* for finite difference include "L_3|[[V,T^B_2],\tau_nu_2]|HF> 522*---------------------------------------------------------------------* 523 524 CALL CC3_L3_OMEGA2_NODDY(WORK(KOME2),WORK(KL3AM), 525 * WORK(KINT1SB),WORK(KINT2SB)) 526 527 IF (NONHF) THEN 528 CALL CCSDT_E2AM(WORK(KOME2),WORK(KL3AM), 529 * WORK(KTB2AM),WORK(KFIELD)) 530 END IF 531 532 DO I = 1,NT1AMX 533 DO J = 1,I 534 IJ = NT1AMX*(I-1) + J 535 NIJ = INDEX(I,J) 536 OMEGA2(NIJ) = OMEGA2(NIJ) + WORK(KOME2+IJ-1) 537 END DO 538 END DO 539 540 IF (LOCDBG) THEN 541 WRITE(LUPRI,*) 'CCSDT_FMAT_NODDY> F matrix result on output:' 542 CALL CC_PRP(OMEGA1,OMEGA2,1,1,1) 543 END IF 544 545*---------------------------------------------------------------------* 546* Compute response Fock matrix F^B(kc) = sum_ia L_kcia T^B_ia 547*---------------------------------------------------------------------* 548 CALL DZERO(WORK(KFOCKB),NORBT*NORBT) 549 550 CALL CCSDT_FCK_R(WORK(KFOCKB),WORK(KXIAJB),WORK(KTB1AM)) 551 552*---------------------------------------------------------------------* 553* Compute triples result vector <L_2|[H^B,\tau_nu_3]|HF>: 554*---------------------------------------------------------------------* 555 KEND3 = KEND2 556 557 IF (NONHF) THEN 558 ! in case of non-HF fields we need to keep L3 559 KEND3 = KL3AM + NT1AMX*NT1AMX*NT1AMX 560 END IF 561 562 KF3AM = KEND3 563 KEND3 = KF3AM + NT1AMX*NT1AMX*NT1AMX 564 565 LWRK3 = LWORK - KEND3 566 IF (LWRK3 .LT. 0) THEN 567 CALL QUIT('Insufficient space in CCSDT_FMAT_NODDY (8)') 568 ENDIF 569 570 CALL DZERO(WORK(KF3AM),NT1AMX*NT1AMX*NT1AMX) 571 572 IF (.NOT.SKIP_T3_CONT) THEN ! might be done via CCSDT_FBC_NODDY 573 CALL CCSDT_F3AM(WORK(KF3AM),WORK(KFOCKB),WORK(KINT1TB), 574 & WORK(KINT2TB),WORK(KL2AM)) 575 END IF 576 577 IF (NONHF) THEN 578c ------------------------------------------- 579c compute one-index transformed field [V,R1]: 580c ------------------------------------------- 581 CALL DCOPY(NORBT*NORBT,WORK(KFIELDAO),1,WORK(KFLDB1),1) 582 CALL CC_FCKMO(WORK(KFLDB1),WORK(KLAMP0),WORK(KLAMHB), 583 * WORK(KEND3),LWRK3,1,1,1) 584 CALL CC_FCKMO(WORK(KFIELDAO),WORK(KLAMPB),WORK(KLAMH0), 585 * WORK(KEND3),LWRK3,1,1,1) 586 CALL DAXPY(NORBT*NORBT,ONE,WORK(KFIELDAO),1,WORK(KFLDB1),1) 587 IF (LOCDBG) WRITE(LUPRI,*) 'Norm^2(FLDB1):', 588 * DDOT(NORBT*NORBT,WORK(KFLDB1),1,WORK(KFLDB1),1) 589 590 L2INCL = .FALSE. 591 KDUM = KEND3 592 CALL CCSDT_E3AM(WORK(KF3AM),WORK(KDUM),WORK(KL3AM), 593 & WORK(KFLDB1),L2INCL) 594 END IF 595 596*---------------------------------------------------------------------* 597* Compute finite difference result and compare: 598*---------------------------------------------------------------------* 599 IF (FD_TEST) THEN 600 KTB2AM = KEND3 601 KTB3AM = KTB2AM + NT1AMX*NT1AMX 602 KFFD1AM = KTB3AM + NT1AMX*NT1AMX*NT1AMX 603 KFFD2AM = KFFD1AM + NT1AMX 604 KFFD3AM = KFFD2AM + NT1AMX*NT1AMX 605 KEND4 = KFFD3AM + NT1AMX*NT1AMX*NT1AMX 606 LWRK4 = LWORK - KEND4 607 IF (LWRK4 .LT. NT2AMX) THEN 608 CALL QUIT('Insufficient space in CCSDT_FMAT_NODDY (9)') 609 ENDIF 610 611 LUTEMP = -1 612 CALL GPOPEN(LUTEMP,'T3BARAM','UNKNOWN',' ','UNFORMATTED', 613 & IDUMMY,.FALSE.) 614 REWIND LUTEMP 615 READ (LUTEMP) (WORK(KTB3AM-1+IDX),IDX=1,NT1AMX*NT1AMX*NT1AMX) 616 CALL GPCLOSE(LUTEMP,'DELETE') 617 618 ! read T^B doubles amplitudes from file and square up 619 ISYMB = ILSTSYM(LISTB,IDLSTB) 620 IOPT = 3 621 Call CC_RDRSP(LISTB,IDLSTB,ISYMB,IOPT,MODEL, 622 & WORK(KTB1AM),WORK(KEND4)) 623 Call CCLR_DIASCL(WORK(KEND4),TWO,ISYMB) 624 CALL CC_T2SQ(WORK(KEND4),WORK(KTB2AM),ISYMB) 625 626 CALL CCSDT_FMAT_FD(WORK(KTB1AM),WORK(KTB2AM),WORK(KTB3AM), 627 & WORK(KFFD1AM),WORK(KFFD2AM),WORK(KFFD3AM), 628 & WORK(KFOCKB),WORK(KFOCK0),WORK(KEND4),LWRK4) 629 630 631 WRITE(LUPRI,*) 'CCSDT_F_NODDY> my F1AM:' 632 CALL OUTPUT(WORK(KOME1),1,NVIRT,1,NRHFT, 633 & NVIRT,NRHFT,1,LUPRI) 634 WRITE(LUPRI,*) 'CCSDT_F_NODDY> finite difference F1AM:' 635 CALL OUTPUT(WORK(KFFD1AM),1,NVIRT,1,NRHFT, 636 & NVIRT,NRHFT,1,LUPRI) 637 CALL DAXPY(NT1AMX,-1.0D0,WORK(KOME1),1,WORK(KFFD1AM),1) 638 WRITE(LUPRI,*) 'CCSDT_F_NODDY> norm of difference:', 639 & DDOT(NT1AMX,WORK(KFFD1AM),1,WORK(KFFD1AM),1) 640 641 WRITE(LUPRI,*) 'CCSDT_F_NODDY> my F2AM:' 642 CALL OUTPUT(WORK(KOME2),1,NT1AMX,1,NT1AMX, 643 & NT1AMX,NT1AMX,1,LUPRI) 644 WRITE(LUPRI,*) 'CCSDT_F_NODDY> finite difference F2AM:' 645 CALL OUTPUT(WORK(KFFD2AM),1,NT1AMX,1,NT1AMX, 646 & NT1AMX,NT1AMX,1,LUPRI) 647 CALL DAXPY(NT1AMX*NT1AMX,-1.0D0,WORK(KOME2),1,WORK(KFFD2AM),1) 648 WRITE(LUPRI,*) 'CCSDT_F_NODDY> norm of difference:', 649 & DDOT(NT1AMX*NT1AMX,WORK(KFFD2AM),1,WORK(KFFD2AM),1) 650 651 WRITE(LUPRI,*) 'CCSDT_F_NODDY> my F3AM:' 652 CALL OUTPUT(WORK(KF3AM),1,NT1AMX*NT1AMX,1,NT1AMX, 653 & NT1AMX*NT1AMX,NT1AMX,1,LUPRI) 654 WRITE(LUPRI,*) 'CCSDT_F_NODDY> finite difference F3AM:' 655 CALL OUTPUT(WORK(KFFD3AM),1,NT1AMX*NT1AMX,1,NT1AMX, 656 & NT1AMX*NT1AMX,NT1AMX,1,LUPRI) 657 CALL DAXPY(NT1AMX*NT1AMX*NT1AMX,-1.0D0,WORK(KF3AM),1, 658 & WORK(KFFD3AM),1) 659 WRITE(LUPRI,*) 'CCSDT_F_NODDY> norm of difference:', 660 & DDOT(NT1AMX*NT1AMX*NT1AMX,WORK(KFFD3AM),1,WORK(KFFD3AM),1) 661 662 WRITE(LUPRI,*) 'CCSDT_F_NODDY> difference vector:' 663 CALL OUTPUT(WORK(KFFD3AM),1,NT1AMX*NT1AMX,1,NT1AMX, 664 & NT1AMX*NT1AMX,NT1AMX,1,LUPRI) 665 666 END IF 667 668*---------------------------------------------------------------------* 669* Now we split: 670* for IOPTRES < 5 we compute the effective F transformed vector 671* for IOPTRES = 5 we compute the contractions F T^B T^C 672*---------------------------------------------------------------------* 673 IF (IOPTRES.GE.1 .AND. IOPTRES.LE.4) THEN 674 675 KT02AM = KEND3 676 KEND3 = KT02AM + NT1AMX*NT1AMX 677 LWRK3 = LWORK - KEND3 678 IF (LWRK3 .LT. NT2AMX) THEN 679 CALL QUIT('Insufficient space in CCSDT_FMAT_NODDY (10)') 680 ENDIF 681 682 ! read T^0 doubles amplitudes from file and square up 683 IOPT = 2 684 KDUM = KEND3 685 Call CC_RDRSP('R0',0,ISYM0,IOPT,MODEL,WORK(KDUM),WORK(KEND3)) 686 CALL CC_T2SQ(WORK(KEND3),WORK(KT02AM),ISYM0) 687 688 CALL DCOPY(NT1AMX,OMEGA1,1,OMEGA1EFF,1) 689 CALL DCOPY(NT2AMX,OMEGA2,1,OMEGA2EFF,1) 690 691 FREQF = FREQL + FREQB 692 693 CALL CC_LHPART_NODDY(OMEGA1EFF,OMEGA2EFF,WORK(KF3AM),-FREQF, 694 & WORK(KFOCKD),WORK(KFIELD), 695 & WORK(K0IOOOO),WORK(K0IOVVO), 696 & WORK(K0IOOVV),WORK(K0IVVVV), 697 & WORK(KT02AM),WORK(KINT1S0),WORK(KINT2S0), 698 & WORK(KEND3),LWRK3) 699 700 701 ELSE IF (IOPTRES.EQ.5) THEN 702 703 SIGN = +1.0D0 704 705 SKIP_F3 = SKIP_T3_CONT .AND. (.NOT.NONHF) 706 707 CALL CCDOTRSP_NODDY(WORK(KOME1),WORK(KOME2),WORK(KF3AM),SIGN, 708 & ITRAN,LISTDP,IDOTS,DOTPROD,MXVEC, 709 & WORK(KLAMP0),WORK(KLAMH0), 710 & WORK(KFOCK0),WORK(KFOCKD), 711 & WORK(KXIAJB), WORK(KYIAJB), 712 & WORK(KINT1T0),WORK(KINT2T0), 713 & WORK(KINT1S0),WORK(KINT2S0), 714 & 'CCSDT_F_NODDY',LOCDBG,LOCDBG,SKIP_F3, 715 & WORK(KEND3),LWRK3) 716 717 ELSE 718 CALL QUIT('Illegal value for IOPTRES IN CCSDT_FMAT_NODDY') 719 END IF 720 721 CALL QEXIT('CCSDT_FMAT_NODDY') 722 723 RETURN 724 END 725*---------------------------------------------------------------------* 726* END OF SUBROUTINE CCSDT_FMAT_NODDY * 727*---------------------------------------------------------------------* 728*=====================================================================* 729 SUBROUTINE CCSDT_F3AM(F3AM,FOCK,XINT1,XINT2,L2AM) 730*---------------------------------------------------------------------* 731* 732* Purpose: compute triples exictation amplitudes of 733* F matrix transformed vector 734* 735* (F T^B)_nu_3 = <L_2|[H^B,tau_nu_3]|HF> 736* 737* Written by Christof Haettig, April 2002 738*=====================================================================* 739#if defined (IMPLICIT_NONE) 740 IMPLICIT NONE 741#else 742# include "implicit.h" 743#endif 744#include "priunit.h" 745#include "ccsdinp.h" 746#include "maxorb.h" 747#include "ccsdsym.h" 748#include "ccfield.h" 749#include "ccorb.h" 750 751#if defined (SYS_CRAY) 752 REAL AIBJCK 753 REAL F3AM(NT1AMX,NT1AMX,NT1AMX) 754 REAL L2AM(NT1AMX,NT1AMX) 755 REAL XINT1(NT1AMX,NVIRT,NVIRT) 756 REAL XINT2(NT1AMX,NRHFT,NRHFT) 757 REAL FOCK(NORBT,NORBT) 758 REAL TWO 759#else 760 DOUBLE PRECISION AIBJCK 761 DOUBLE PRECISION F3AM(NT1AMX,NT1AMX,NT1AMX) 762 DOUBLE PRECISION L2AM(NT1AMX,NT1AMX) 763 DOUBLE PRECISION XINT1(NT1AMX,NVIRT,NVIRT) 764 DOUBLE PRECISION XINT2(NT1AMX,NRHFT,NRHFT) 765 DOUBLE PRECISION FOCK(NORBT,NORBT) 766 DOUBLE PRECISION TWO 767#endif 768 PARAMETER ( TWO = 2.0D0 ) 769 770 INTEGER NI, NA, NAI, NK, NC, NCK, NJ, NCJ, NB, NBJ, NBK, NBI, 771 & NAJ, ND, NDJ, NDK, NL, NAL, NBL, NCI, NAK 772C 773 DO NK = 1,NRHFT 774 DO NC = 1,NVIRT 775 NCK = NVIRT*(NK-1) + NC 776C 777 DO NJ = 1,NRHFT 778 NCJ = NVIRT*(NJ-1) + NC 779 780 DO NB = 1,NVIRT 781 NBJ = NVIRT*(NJ-1) + NB 782 NBK = NVIRT*(NK-1) + NB 783C 784 DO NI = 1,NRHFT 785 NBI = NVIRT*(NI-1) + NB 786 787 DO NA = 1,NVIRT 788 NAI = NVIRT*(NI-1) + NA 789 NAJ = NVIRT*(NJ-1) + NA 790C 791 AIBJCK = 0.0D0 792 793 AIBJCK = AIBJCK + FOCK(NK,NRHFT+NC) * L2AM(NAI,NBJ) 794C 795 AIBJCK = AIBJCK - FOCK(NJ,NRHFT+NC) * L2AM(NAI,NBK) 796C 797 DO ND = 1,NVIRT 798 NDJ = NVIRT*(NJ-1) + ND 799 NDK = NVIRT*(NK-1) + ND 800C 801 AIBJCK = AIBJCK + 802 & (TWO*XINT1(NCK,ND,NB)-XINT1(NBK,ND,NC))* 803 & L2AM(NDJ,NAI) 804C 805 AIBJCK = AIBJCK - XINT1(NBI,ND,NC) * L2AM(NDK,NAJ) 806 END DO 807C 808 DO NL = 1,NRHFT 809 NAL = NVIRT*(NL-1) + NA 810 NBL = NVIRT*(NL-1) + NB 811C 812 AIBJCK = AIBJCK - 813 & (TWO*XINT2(NCK,NJ,NL)-XINT2(NCJ,NK,NL))* 814 & L2AM(NBL,NAI) 815C 816 AIBJCK = AIBJCK + XINT2(NCJ,NI,NL) * L2AM(NBK,NAL) 817 END DO 818C 819 F3AM(NAI,NBJ,NCK) = F3AM(NAI,NBJ,NCK) + AIBJCK 820 F3AM(NAI,NCK,NBJ) = F3AM(NAI,NCK,NBJ) + AIBJCK 821 F3AM(NBJ,NAI,NCK) = F3AM(NBJ,NAI,NCK) + AIBJCK 822 F3AM(NCK,NAI,NBJ) = F3AM(NCK,NAI,NBJ) + AIBJCK 823 F3AM(NBJ,NCK,NAI) = F3AM(NBJ,NCK,NAI) + AIBJCK 824 F3AM(NCK,NBJ,NAI) = F3AM(NCK,NBJ,NAI) + AIBJCK 825C 826 END DO 827 END DO 828 END DO 829 END DO 830 END DO 831 END DO 832C 833C------------------------------------------------------ 834C Get rid of amplitudes that are not allowed. 835C------------------------------------------------------ 836 837 CALL CCSDT_CLEAN_T3(F3AM,NT1AMX,NVIRT,NRHFT) 838 839 RETURN 840 END 841*---------------------------------------------------------------------* 842* END OF SUBROUTINE CCSDT_F3AM * 843*---------------------------------------------------------------------* 844*=====================================================================* 845 SUBROUTINE CCSDT_FMAT_FD(TB1AM,TB2AM,TB3AM,F1AM,F2AM,F3AM, 846 & FOCKA,FOCK0,WORK,LWORK) 847*---------------------------------------------------------------------* 848* 849* Purpose: Construct F matrix transformed vector 850* by finite difference on left hand transformation 851* 852* Written by Christof Haettig, April 2002 853*=====================================================================* 854#if defined (IMPLICIT_NONE) 855 IMPLICIT NONE 856#else 857# include "implicit.h" 858#endif 859#include "priunit.h" 860#include "ccsdinp.h" 861#include "maxorb.h" 862#include "ccsdsym.h" 863#include "ccorb.h" 864 865#if defined (SYS_CRAY) 866 REAL DELTA, ZERO 867#else 868 DOUBLE PRECISION DELTA, ZERO, ONE 869#endif 870 PARAMETER( DELTA = 1.0D-6, ONE = 1.0D0, ZERO = 0.0D0) 871 872 INTEGER LWORK 873 874#if defined (SYS_CRAY) 875 REAL WORK(LWORK), 876 & FOCK0(NORBT,NORBT), FOCKA(NORBT,NORBT), DDOT, 877 & TB1AM(NT1AMX),TB2AM(NT1AMX,NT1AMX),TB3AM(NT1AMX,NT1AMX,NT1AMX), 878 & F1AM(NT1AMX), F2AM(NT1AMX,NT1AMX), F3AM(NT1AMX,NT1AMX,NT1AMX) 879#else 880 DOUBLE PRECISION WORK(LWORK), 881 & FOCK0(NORBT,NORBT), FOCKA(NORBT,NORBT), DDOT, 882 & TB1AM(NT1AMX),TB2AM(NT1AMX,NT1AMX),TB3AM(NT1AMX,NT1AMX,NT1AMX), 883 & F1AM(NT1AMX), F2AM(NT1AMX,NT1AMX), F3AM(NT1AMX,NT1AMX,NT1AMX) 884#endif 885 886 CHARACTER*10 MODEL 887 INTEGER KT1AM, KT2AM, KT3AM, KL1AM, KL2AM, KL3AM, KFOCK, 888 & KLAMDP, KLAMDH, KEND1, LWRK1, LUSIFC, IOPT, IDUMMY 889 890 CALL QUIT('CCSDT_FMAT_FD not adapted for intermediates!') 891*---------------------------------------------------------------------* 892* memory allocations: 893*---------------------------------------------------------------------* 894 KT1AM = 1 895 KT2AM = KT1AM + NT1AMX 896 KT3AM = KT2AM + NT1AMX*NT1AMX 897 KL1AM = KT3AM + NT1AMX*NT1AMX*NT1AMX 898 KL2AM = KL1AM + NT1AMX 899 KL3AM = KL2AM + NT1AMX*NT1AMX 900 KFOCK = KL3AM + NT1AMX*NT1AMX*NT1AMX 901 KLAMDP= KFOCK + NORBT*NORBT 902 KLAMDH= KLAMDP+ NLAMDT 903 KEND1 = KLAMDH+ NLAMDT 904 LWRK1 = LWORK - KEND1 905 IF (LWRK1 .LT. 0) THEN 906 CALL QUIT('Insufficient space in CCSDT_FD') 907 ENDIF 908 909 WRITE(LUPRI,*) 'CCSDT_FMAT_FD> norm^2(TB3AM):', 910 & DDOT(NT1AMX**3,TB3AM,1,TB3AM,1) 911 912*---------------------------------------------------------------------* 913* Get T0 singles & doubles and L0 singles & doubles and 914* zeroth-order lambda matrices: 915*---------------------------------------------------------------------* 916 IOPT = 3 917 Call CC_RDRSP('R0',0,1,IOPT,MODEL,WORK(KT1AM),WORK(KT3AM)) 918 CALL CC_T2SQ(WORK(KT3AM),WORK(KT2AM),1) 919 920 IOPT = 3 921 Call CC_RDRSP('L0',0,1,IOPT,MODEL,WORK(KL1AM),WORK(KL3AM)) 922 CALL CC_T2SQ(WORK(KL3AM),WORK(KL2AM),1) 923 924 CALL LAMMAT(WORK(KLAMDP),WORK(KLAMDH),WORK(KT1AM), 925 & WORK(KEND1),LWRK1) 926 927 IOPT = 1 928 CALL CC_LHTR_NODDY(0.0D0,F1AM,F2AM,WORK(KT1AM),WORK(KT2AM), 929 & WORK(KL1AM),WORK(KL2AM), 930 & WORK(KLAMDP),WORK(KLAMDH), 931 & WORK(KEND1),LWRK1,IOPT) 932 933*---------------------------------------------------------------------* 934* Get T0 singles, doubles & triples 935*---------------------------------------------------------------------* 936 IOPT = 3 937 Call CC_RDRSP('R0',0,1,IOPT,MODEL,WORK(KT1AM),WORK(KT3AM)) 938 CALL CC_T2SQ(WORK(KT3AM),WORK(KT2AM),1) 939 940 LUSIFC = -1 941 CALL GPOPEN(LUSIFC,'T3AMPL','UNKNOWN',' ','UNFORMATTED',IDUMMY, 942 & .FALSE.) 943 REWIND LUSIFC 944 READ (LUSIFC) (WORK(KT3AM-1+I), I=1,NT1AMX*NT1AMX*NT1AMX) 945 CALL GPCLOSE(LUSIFC,'KEEP') 946 947 WRITE(LUPRI,*) 'CCSDT_FMAT_FD> norm^2(T3AM):', 948 & DDOT(NT1AMX**3,WORK(KT3AM),1,WORK(KT3AM),1) 949 950*---------------------------------------------------------------------* 951* Get L0 singles, doubles & triples 952*---------------------------------------------------------------------* 953 IOPT = 3 954 Call CC_RDRSP('L0',0,1,IOPT,MODEL,WORK(KL1AM),WORK(KL3AM)) 955 CALL CC_T2SQ(WORK(KL3AM),WORK(KL2AM),1) 956 957 LUSIFC = -1 958 CALL GPOPEN(LUSIFC,'L3AMPL','UNKNOWN',' ','UNFORMATTED',IDUMMY, 959 & .FALSE.) 960 REWIND LUSIFC 961 READ (LUSIFC) (WORK(KL3AM-1+I), I=1,NT1AMX*NT1AMX*NT1AMX) 962 CALL GPCLOSE(LUSIFC,'KEEP') 963 964 WRITE(LUPRI,*) 'CCSDT_FMAT_FD> norm^2(L3AM):', 965 & DDOT(NT1AMX**3,WORK(KL3AM),1,WORK(KL3AM),1) 966 967*---------------------------------------------------------------------* 968* Compute finite difference vector: 969*---------------------------------------------------------------------* 970 CALL DZERO(F1AM,NT1AMX ) 971 CALL DZERO(F2AM,NT1AMX*NT1AMX ) 972 CALL DZERO(F3AM,NT1AMX*NT1AMX*NT1AMX) 973 974 CALL DCOPY(NORBT*NORBT,FOCK0,1,WORK(KFOCK),1) 975 976 CALL DAXPY(NT1AMX, -DELTA,TB1AM,1,WORK(KT1AM),1) 977 CALL DAXPY(NT1AMX*NT1AMX, -DELTA,TB2AM,1,WORK(KT2AM),1) 978 CALL DAXPY(NT1AMX*NT1AMX*NT1AMX,-DELTA,TB3AM,1,WORK(KT3AM),1) 979 980 CALL DAXPY(NORBT*NORBT,-DELTA,FOCKA,1,WORK(KFOCK),1) 981 982 CALL LAMMAT(WORK(KLAMDP),WORK(KLAMDH),WORK(KT1AM), 983 & WORK(KEND1),LWRK1) 984 985 986 CALL CC_LHTR_NODDY2(F1AM,F2AM,F3AM, 987 & WORK(KT1AM),WORK(KT2AM),WORK(KT3AM), 988 & WORK(KL1AM),WORK(KL2AM),WORK(KL3AM), 989 & WORK(KLAMDP),WORK(KLAMDH), 990 & WORK(KFOCK),WORK(KEND1),LWRK1) 991 992 CALL DSCAL(NT1AMX, -ONE,F1AM,1) 993 CALL DSCAL(NT1AMX*NT1AMX, -ONE,F2AM,1) 994 CALL DSCAL(NT1AMX*NT1AMX*NT1AMX,-ONE,F3AM,1) 995 996 CALL DAXPY(NT1AMX, 2.D0*DELTA,TB1AM,1,WORK(KT1AM),1) 997 CALL DAXPY(NT1AMX*NT1AMX, 2.D0*DELTA,TB2AM,1,WORK(KT2AM),1) 998 CALL DAXPY(NT1AMX*NT1AMX*NT1AMX,2.D0*DELTA,TB3AM,1,WORK(KT3AM),1) 999 1000 CALL DAXPY(NORBT*NORBT,2.D0*DELTA,FOCKA,1,WORK(KFOCK),1) 1001 1002 CALL LAMMAT(WORK(KLAMDP),WORK(KLAMDH),WORK(KT1AM), 1003 & WORK(KEND1),LWRK1) 1004 1005 CALL CC_LHTR_NODDY2(F1AM,F2AM,F3AM, 1006 & WORK(KT1AM),WORK(KT2AM),WORK(KT3AM), 1007 & WORK(KL1AM),WORK(KL2AM),WORK(KL3AM), 1008 & WORK(KLAMDP),WORK(KLAMDH), 1009 & WORK(KFOCK),WORK(KEND1),LWRK1) 1010 1011 1012 CALL DSCAL(NT1AMX ,0.5D0/DELTA,F1AM,1) 1013 CALL DSCAL(NT1AMX*NT1AMX ,0.5D0/DELTA,F2AM,1) 1014 CALL DSCAL(NT1AMX*NT1AMX*NT1AMX,0.5D0/DELTA,F3AM,1) 1015 1016 RETURN 1017 END 1018*---------------------------------------------------------------------* 1019* END OF SUBROUTINE CCSDT_F_FD * 1020*---------------------------------------------------------------------* 1021*=====================================================================* 1022 SUBROUTINE CCSDT_TBAR31_NODDY(TBAR3AM,FREQB,LISTB,IDLSTB, 1023 & XLAMDP,XLAMDH,FOCK0,FOCKD,SCR1, 1024 & XIAJB,XINT1T0,XINT2T0, 1025 & WORK,LWORK) 1026*---------------------------------------------------------------------* 1027* 1028* Purpose: compute triples part of first-order response multipliers 1029* 1030* Implemented for L1, X1B, E0, N2, M1 and LE vectors. 1031* 1032* Input: LISTB, IDLSTB 1033* XLAMDP, XLAMDH 1034* FOCK0, FOCKD 1035* XIAJB 1036* XINT1T0, XINT2T0 1037* 1038* Output: TBAR3AM 1039* FREQB 1040* 1041* Scratch: SCR1 1042* 1043* Written by Christof Haettig, April 2002 1044* 1045*=====================================================================* 1046 IMPLICIT NONE 1047#include "priunit.h" 1048#include "ccsdinp.h" 1049#include "maxorb.h" 1050#include "ccsdsym.h" 1051#include "ccorb.h" 1052#include "ccl1rsp.h" 1053#include "ccn2rsp.h" 1054#include "cclrmrsp.h" 1055#include "ccexgr.h" 1056#include "ccexci.h" 1057#include "ccfield.h" 1058 1059 LOGICAL LOCDBG, PRINT_T3 1060 PARAMETER (LOCDBG=.false.,PRINT_T3=.false.) 1061 1062 1063 CHARACTER LISTB*3 1064 LOGICAL RHS_ONLY 1065 INTEGER LWORK, IDLSTB 1066#if defined (SYS_CRAY) 1067 REAL WORK(LWORK), FOCKD(*), FOCK0(*), SCR1(*) 1068 REAL FREQB, DUMMY, DDOT, ZERO, FF, ONE 1069 REAL XLAMDP(*), XLAMDH(*) 1070 REAL TBAR3AM(NT1AMX,NT1AMX,NT1AMX) 1071 REAL XINT1T0(*), XINT2T0(*), XIAJB(*), YIAJB(*) 1072#else 1073 DOUBLE PRECISION WORK(LWORK), FOCKD(*), FOCK0(*), SCR1(*) 1074 DOUBLE PRECISION FREQB, DUMMY, DDOT, ZERO, FF, ONE 1075 DOUBLE PRECISION XLAMDP(*), XLAMDH(*) 1076 DOUBLE PRECISION TBAR3AM(NT1AMX,NT1AMX,NT1AMX) 1077 DOUBLE PRECISION XINT1T0(*), XINT2T0(*), XIAJB(*) 1078#endif 1079 LOGICAL L2INCL 1080 INTEGER ISYM0 1081 PARAMETER( L2INCL=.TRUE. , ISYM0=1, ZERO=0.0D0, ONE=1.0D0) 1082 1083 CHARACTER MODEL*10, LABELB*8, LISTR*3, LISTL*3 1084 LOGICAL LORXB, DO_INTR 1085 INTEGER ISYMB, KFOCKB, KL1AM, KL2AM, KL3AM, KEND1, LWRK1, 1086 & IOPT, IRREP, IERR, NAI, NBJ, NCK, NAIBJCK, ILSTNR, 1087 & KTB1AM, KLAMPB, KLAMHB, KINT1TB, KINT2TB, ISYMD, ILLL, 1088 & IDEL, ISYDIS, KXINT, KEND2, LWRK2, IRECNR, IR1TAMP, 1089 & ILSTSYM, ILSTNL, ISYML, ISYMR, IMAT, KFIELD, KFIELDAO, 1090 & KFLDB1, KT3SCR 1091 1092 1093 CALL QENTER('CCTBAR31') 1094*---------------------------------------------------------------------* 1095* Do some initial tests and memory allocations: 1096*---------------------------------------------------------------------* 1097 if (locdbg) then 1098 write(lupri,*) 'entered ccsdt_tbar31_noddy for "',LISTB(1:3), 1099 & '" ',idlstb 1100 end if 1101 1102 if (NSYM.NE.1) CALL QUIT('NO SYMMETRY FOR CCSDT_TBAR31_NODDY') 1103 1104 ISYMB = ILSTSYM(LISTB,IDLSTB) 1105 1106c ------------------------------------------------------------ 1107c first-order response of ground state lagrangian multipliers: 1108c ------------------------------------------------------------ 1109 IF (LISTB(1:3).EQ.'L1 '.OR.LISTB(1:3).EQ.'X1B') THEN 1110 RHS_ONLY = LISTB.EQ.'X1B' 1111 1112 ! get symmetry, frequency and integral label from common blocks 1113 ! defined in ccl1rsp.h 1114 FREQB = FRQLRZ(IDLSTB) 1115 LABELB = LRZLBL(IDLSTB) 1116 LORXB = LORXLRZ(IDLSTB) 1117 1118 IF (LORXB) CALL QUIT('NO ORBITAL RELAX. IN CCSDT_TBAR31_NODDY') 1119 1120 ! find corresponding first-order right amplitudes 1121 LISTR = 'R1 ' 1122 ILSTNR = IR1TAMP(LABELB,LORXB,FREQB,ISYMB) 1123 1124 ! set zero-order left vector 1125 LISTL = 'L0 ' 1126 ILSTNL = 0 1127 1128c ------------------------------------------------------- 1129c zero-order lagrangian multipliers for ground to excited 1130c state transitions: 1131c ------------------------------------------------------- 1132 ELSE IF (LISTB(1:3).EQ.'M1 ') THEN 1133 RHS_ONLY = .FALSE. 1134 1135 FREQB = FRQLRM(IDLSTB) 1136 LABELB = '- none -' 1137 1138 ! find corresponding right eigenvector 1139 LISTR = 'RE ' 1140 ILSTNR = ILRM(IDLSTB) 1141 1142 ! set zero-order left vector 1143 LISTL = 'L0 ' 1144 ILSTNL = 0 1145 1146c -------------------------------------------------------- 1147c zero-order lagrangian multipliers for excited states: 1148c (Diagonal case for N2, but with Tbar^(0) added. The 1149c latter does here only contribute indirectly via the 1150c singles and doubles part of the solution vector "E0") 1151c -------------------------------------------------------- 1152 ELSE IF (LISTB(1:3).EQ.'E0 ') THEN 1153 RHS_ONLY = .FALSE. 1154 1155 FREQB = ZERO 1156 LABELB = '- none -' 1157 1158 ! find corresponding right eigenvector 1159 LISTR = 'RE ' 1160 ILSTNR = IXGRST(IDLSTB) 1161 1162 ! set zero-order left vector 1163 LISTL = 'LE ' 1164 ILSTNL = IXGRST(IDLSTB) 1165 1166c -------------------------------------------------------- 1167c zero-order lagrangian multipliers for excited to excited 1168c state transitions: 1169c -------------------------------------------------------- 1170 ELSE IF (LISTB(1:3).EQ.'N2 ') THEN 1171 RHS_ONLY = .FALSE. 1172 1173 FREQB = FRQIN2(IDLSTB) + FRQFN2(IDLSTB) 1174 LABELB = '- none -' 1175 1176 ! find corresponding right eigenvector 1177 LISTR = 'RE ' 1178 ILSTNR = IFN2(IDLSTB) 1179 1180 ! set zero-order left vector 1181 LISTL = 'LE ' 1182 ILSTNL = IIN2(IDLSTB) 1183 1184c ----------------------------- 1185c zero-order left eigenvectors: 1186c ----------------------------- 1187 ELSE IF (LISTB(1:3).EQ.'LE ') THEN 1188 RHS_ONLY = .FALSE. 1189 1190 FREQB = -EIGVAL(IDLSTB) 1191 LABELB = '- none -' 1192 1193 ! we don't have any "right" vector entering a right hand side 1194 LISTR = '---' 1195 ILSTNR = -99 1196 1197 ! we don't have any "zero-order left" vector entering a rh side 1198 LISTL = '---' 1199 ILSTNL = -99 1200 1201 ELSE 1202 WRITE(LUPRI,*) 'LISTB,IDLSTB:',LISTB(1:3),IDLSTB 1203 CALL QUIT('Unknown/Illegal list in CCSDT_TBAR31_NODDY.') 1204 END IF 1205 1206 IF (LISTB(1:3).NE.'LE ') THEN 1207 DO_INTR = .TRUE. 1208 ISYML = ILSTSYM(LISTL,ILSTNL) 1209 ISYMR = ILSTSYM(LISTR,ILSTNR) 1210 IF (MULD2H(ISYML,ISYMR).NE.ISYMB) CALL QUIT( 1211 & 'Symmetry mismatch in CCSDT_TBAR31_NODDY.') 1212 ELSE 1213 DO_INTR = .FALSE. 1214 END IF 1215 1216 1217 KL1AM = 1 1218 KL2AM = KL1AM + NT1AMX 1219 KL3AM = KL2AM + NT1AMX*NT1AMX 1220 KEND1 = KL3AM + NT1AMX*NT1AMX*NT1AMX 1221 1222 KTB1AM = KEND1 1223 KFOCKB = KTB1AM + NT1AMX 1224 KLAMPB = KFOCKB + NORBT*NORBT 1225 KLAMHB = KLAMPB + NLAMDT 1226 KEND1 = KLAMHB + NLAMDT 1227 1228 IF (NONHF) THEN 1229 KFIELD = KEND1 1230 KFLDB1 = KFIELD + NORBT*NORBT 1231 KFIELDAO = KFLDB1 + NORBT*NORBT 1232 KEND1 = KFIELDAO + NORBT*NORBT 1233 END IF 1234 1235 KINT1TB = KEND1 1236 KINT2TB = KINT1TB + NT1AMX*NVIRT*NVIRT 1237 KEND1 = KINT2TB + NRHFT*NRHFT*NT1AMX 1238 1239 LWRK1 = LWORK - KEND1 1240 IF (LWRK1 .LT. NT2AMX) THEN 1241 CALL QUIT('Insufficient space in CCSDT_TBAR31_NODDY') 1242 ENDIF 1243 1244*---------------------------------------------------------------------* 1245* Calculate response Lambda matrices from the right vectors: 1246*---------------------------------------------------------------------* 1247 IF (DO_INTR) THEN 1248 IOPT = 1 1249 CALL CC_RDRSP(LISTR,ILSTNR,ISYMB,IOPT,MODEL,WORK(KTB1AM),WORK) 1250 1251 CALL CCLR_LAMTRA(XLAMDP,WORK(KLAMPB), 1252 & XLAMDH,WORK(KLAMHB),WORK(KTB1AM),ISYMB) 1253 END IF 1254 1255*---------------------------------------------------------------------* 1256* Get the one electron integrals from the fields 1257*---------------------------------------------------------------------* 1258 IF (NONHF) THEN 1259 CALL DZERO(WORK(KFIELDAO),NORBT*NORBT) 1260 DO I = 1, NFIELD 1261 FF = EFIELD(I) 1262 CALL CC_ONEP(WORK(KFIELDAO),WORK(KEND1),LWRK1, 1263 * FF,1,LFIELD(I)) 1264 ENDDO 1265 CALL DCOPY(NORBT*NORBT,WORK(KFIELDAO),1,WORK(KFIELD),1) 1266 1267 ! calculate external field in zero-order lambda basis 1268 CALL CC_FCKMO(WORK(KFIELD),XLAMDP,XLAMDH, 1269 * WORK(KEND1),LWRK1,1,1,1) 1270 1271 ! one-index transformed field integrals for [V,R1] 1272 CALL DCOPY(NORBT*NORBT,WORK(KFIELDAO),1,WORK(KFLDB1),1) 1273 CALL CC_FCKMO(WORK(KFLDB1),XLAMDP,WORK(KLAMHB), 1274 * WORK(KEND1),LWRK1,1,1,1) 1275 1276 CALL CC_FCKMO(WORK(KFIELDAO),WORK(KLAMPB),XLAMDH, 1277 * WORK(KEND1),LWRK1,1,1,1) 1278 CALL DAXPY(NORBT*NORBT,ONE,WORK(KFIELDAO),1,WORK(KFLDB1),1) 1279 ENDIF 1280 1281*---------------------------------------------------------------------* 1282* Compute integrals needed for the following contributions: 1283*---------------------------------------------------------------------* 1284 IF (DO_INTR) THEN 1285 1286 CALL DZERO(WORK(KINT1TB),NT1AMX*NVIRT*NVIRT) 1287 CALL DZERO(WORK(KINT2TB),NRHFT*NRHFT*NT1AMX) 1288 1289 DO ISYMD = 1, NSYM 1290 DO ILLL = 1,NBAS(ISYMD) 1291 IDEL = IBAS(ISYMD) + ILLL 1292 ISYDIS = MULD2H(ISYMD,ISYMOP) 1293 1294C ---------------------------- 1295C Work space allocation no. 2. 1296C ---------------------------- 1297 KXINT = KEND1 1298 KEND2 = KXINT + NDISAO(ISYDIS) 1299 LWRK2 = LWORK - KEND2 1300 IF (LWRK2 .LT. 0) THEN 1301 WRITE(LUPRI,*) 'Need : ',KEND2,'Available : ',LWORK 1302 CALL QUIT('Insufficient space in CCSDT_TBAR31_NODDY') 1303 ENDIF 1304 1305C --------------------------- 1306C Read in batch of integrals. 1307C --------------------------- 1308 CALL CCRDAO(WORK(KXINT),IDEL,1,WORK(KEND2),LWRK2, 1309 & IRECNR,DIRECT) 1310 1311C ---------------------------------- 1312C Calculate integrals needed in CC3: 1313C ---------------------------------- 1314 CALL CCSDT_TRAN1_R(WORK(KINT1TB),WORK(KINT2TB), 1315 & XLAMDP,XLAMDH, 1316 & WORK(KLAMPB),WORK(KLAMHB), 1317 & WORK(KXINT),IDEL) 1318 END DO 1319 END DO 1320 END IF 1321 1322*---------------------------------------------------------------------* 1323* triples part of the zero-order lagrangian multipliers: 1324* (note: in case of non-HF external fields TBAR3AM is used 1325* here as scratch array) 1326*---------------------------------------------------------------------* 1327 IF ( LISTB(1:3).EQ.'L1 ' .OR. LISTB(1:3).EQ.'X1B' .OR. 1328 & (LISTB(1:3).NE.'LE '.AND.NONHF) ) THEN 1329 1330C IF (.NOT.( LISTB(1:3).EQ.'L1 ' .OR. LISTB(1:3).EQ.'X1B' )) 1331C & CALL QUIT('CCSDT_TBAR31_NODDY not tested for this case!') 1332 1333 IF (LWRK1 .LT. NT2AMX) 1334 & CALL QUIT('Insufficient space in CCSDT_TBAR31_NODDY') 1335 1336 ! read L^0 multipliers from file and square up doubles part 1337 IOPT = 3 1338 Call CC_RDRSP('L0 ',0,ISYM0,IOPT,MODEL, 1339 & WORK(KL1AM),WORK(KEND1)) 1340 CALL CC_T2SQ(WORK(KEND1),WORK(KL2AM),ISYM0) 1341 1342 ! remember that CCSDT_L03AM returns -L3 !! 1343 CALL CCSDT_L03AM(WORK(KL3AM),XINT1T0,XINT2T0,XIAJB,FOCK0, 1344 & WORK(KL1AM),WORK(KL2AM),SCR1,FOCKD, 1345 & WORK(KFIELD),TBAR3AM) 1346 1347 CALL DSCAL(NT1AMX*NT1AMX*NT1AMX,-1.0D0,WORK(KL3AM),1) 1348 END IF 1349 1350*---------------------------------------------------------------------* 1351* Initialize result vector TBAR3AM: 1352*---------------------------------------------------------------------* 1353 CALL DZERO(TBAR3AM,NT1AMX*NT1AMX*NT1AMX) 1354 1355*=====================================================================* 1356* Compute Eta_3 contribution to Tbar^B_3 (only for "L1" vectors): 1357* for this contribution we need to compute first the triples part 1358* of the ground state zero-order lagrangian multipliers L^0_3 1359* and one-electron property integral matrix: 1360*=====================================================================* 1361 IF (LISTB(1:3).EQ.'L1 '.OR.LISTB(1:3).EQ.'X1B') THEN 1362 1363c ------------------------------------------- 1364c Matrix with property integrals in MO basis: 1365c ------------------------------------------- 1366 ! read property integrals from file: 1367 CALL CCPRPAO(LABELB,.TRUE.,WORK(KFOCKB),IRREP,IMAT,IERR, 1368 & WORK(KEND1),LWRK1) 1369 IF ((IERR.GT.0) .OR. (IERR.EQ.0 .AND. IRREP.NE.ISYMB)) THEN 1370 WRITE(LUPRI,*) 'IERR,IRREP,ISYMB:',IERR,IRREP,ISYMB 1371 CALL QUIT('CCSDT_TBAR31_NODDY: error reading oper. '//LABELB) 1372 ELSE IF (IERR.LT.0) THEN 1373 CALL DZERO(WORK(KFOCKB),N2BST(ISYMB)) 1374 END IF 1375 1376 ! transform property integrals to Lambda-MO basis 1377 CALL CC_FCKMO(WORK(KFOCKB),XLAMDP,XLAMDH, 1378 & WORK(KEND1),LWRK1,ISYMB,1,1) 1379 1380c ------------------------------- 1381c Eta_3 contribution to Tbar^B_3: 1382c ------------------------------- 1383 CALL CCSDT_E3AM(TBAR3AM,WORK(KL2AM),WORK(KL3AM), 1384 & WORK(KFOCKB),L2INCL) 1385 1386 IF (LOCDBG) THEN 1387 WRITE(LUPRI,*) 'CCSDT_TBAR31_NODDY> norm^2(E3AM):', 1388 & DDOT(NT1AMX*NT1AMX*NT1AMX,TBAR3AM,1,TBAR3AM,1) 1389 END IF 1390 1391 END IF 1392 1393*=====================================================================* 1394* Compute F matrix (FT^B)_3 contribution to Tbar^B_3 (for all 1395* vectors, apart from the eigenvectors "LE"): 1396* (Note, that here we overwrite FOCKB with a new matrix!!!) 1397*=====================================================================* 1398 IF (LISTB(1:3).NE.'LE ') THEN 1399 ! read zero-order left vector from file and square up 1400 IOPT = 3 1401 Call CC_RDRSP(LISTL,ILSTNL,ISYML,IOPT,MODEL, 1402 & WORK(KL1AM),WORK(KEND1)) 1403 CALL CC_T2SQ(WORK(KEND1),WORK(KL2AM),ISYM0) 1404 1405 CALL DZERO(WORK(KFOCKB),NORBT*NORBT) 1406 1407 CALL CCSDT_FCK_R(WORK(KFOCKB),XIAJB,WORK(KTB1AM)) 1408 1409 CALL CCSDT_F3AM(TBAR3AM,WORK(KFOCKB),WORK(KINT1TB), 1410 & WORK(KINT2TB),WORK(KL2AM)) 1411 1412 IF (NONHF) THEN 1413 CALL CCSDT_E3AM(TBAR3AM,DUMMY,WORK(KL3AM), 1414 & WORK(KFLDB1),.FALSE.) 1415 END IF 1416 1417 IF (LOCDBG) THEN 1418 WRITE(LUPRI,*) 'CCSDT_TBAR31_NODDY> norm^2(F3AM+E3AM):', 1419 & DDOT(NT1AMX*NT1AMX*NT1AMX,TBAR3AM,1,TBAR3AM,1) 1420 END IF 1421 1422 END IF 1423 1424*=====================================================================* 1425* Compute contribution from jacobian transformation: 1426* <Tbar^B_1|[H,tau_nu_3]|HF>/eps_3 1427* + <Tbar^B_2|[H,tau_nu_3]|HF>/eps_3 1428* Note: here we overwrite L1AM and L2AM with response multipliers 1429*=====================================================================* 1430 IF (.NOT. RHS_ONLY) THEN 1431 IOPT = 3 1432 Call CC_RDRSP(LISTB,IDLSTB,ISYMB,IOPT,MODEL, 1433 & WORK(KL1AM),WORK(KEND1)) 1434 if (locdbg) write(lupri,*)'ccsdt_tbar31_noddy> norm^2 of L2:', 1435 & ddot(nt2amx,work(kend1),1,work(kend1),1) 1436 CALL CC_T2SQ(WORK(KEND1),WORK(KL2AM),ISYMB) 1437 1438 IF (LOCDBG) THEN 1439 WRITE(LUPRI,*) 'CCSDT_TBAR31_NODDY> norm^2(F3AM+E3AM)-2:', 1440 & DDOT(NT1AMX*NT1AMX*NT1AMX,TBAR3AM,1,TBAR3AM,1) 1441 END IF 1442 1443 CALL CCSDT_L3AM_R(TBAR3AM,-FREQB,XINT1T0,XINT2T0,XIAJB,FOCK0, 1444 & WORK(KL1AM),WORK(KL2AM),SCR1,FOCKD,.FALSE.) 1445 1446 IF (LOCDBG) THEN 1447 WRITE(LUPRI,*) 'CCSDT_TBAR31_NODDY> FREQB:',FREQB 1448 WRITE(LUPRI,*) 'CCSDT_TBAR31_NODDY> norm^2(XIAJB):', 1449 & DDOT(NT1AMX*NT1AMX,XIAJB,1,XIAJB,1) 1450 WRITE(LUPRI,*) 'CCSDT_TBAR31_NODDY> norm^2(XINT1T0):', 1451 & DDOT(NT1AMX*NVIRT*NVIRT,XINT1T0,1,XINT1T0,1) 1452 WRITE(LUPRI,*) 'CCSDT_TBAR31_NODDY> norm^2(XINT2T0):', 1453 & DDOT(NT1AMX*NRHFT*NRHFT,XINT1T0,1,XINT1T0,1) 1454 WRITE(LUPRI,*) 'CCSDT_TBAR31_NODDY> norm^2(FOCK0):', 1455 & DDOT(NORBT*NORBT,FOCK0,1,FOCK0,1) 1456 WRITE(LUPRI,*) 1457 & 'CCSDT_TBAR31_NODDY> norm^2(F3AM+E3AM+L3AM(Tbar)):', 1458 & DDOT(NT1AMX*NT1AMX*NT1AMX,TBAR3AM,1,TBAR3AM,1) 1459 END IF 1460 END IF 1461 1462*=====================================================================* 1463* Divide by orbital energy differences and frequency: 1464* (for finite difference we solve the equations iteratively) 1465*=====================================================================* 1466 IF (RHS_ONLY) THEN 1467 CALL DSCAL(NT1AMX*NT1AMX*NT1AMX,-1.0D0,TBAR3AM,1) 1468 ELSE 1469 KT3SCR = KL3AM ! scratch array for NONHF case 1470 1471 CALL CCSDT_3AM(TBAR3AM,-FREQB,SCR1,FOCKD, 1472 & NONHF,WORK(KFIELD),.TRUE.,WORK(KT3SCR)) 1473 1474 END IF 1475*---------------------------------------------------------------------* 1476* clean forbidden elements, print, if required and return: 1477*---------------------------------------------------------------------* 1478 CALL CCSDT_CLEAN_T3(TBAR3AM,NT1AMX,NVIRT,NRHFT) 1479 1480 IF (PRINT_T3) THEN 1481 WRITE(LUPRI,*)'CCSDT_TBAR31_AM> first-order T3-bar vector:' 1482 WRITE(LUPRI,*)'CCSDT_TBAR31_AM> list,idlst:',LISTB,IDLSTB 1483 WRITE(LUPRI,*)'CCSDT_TBAR31_AM> freq,label:',FREQB,LABELB 1484 CALL PRINT_PT3_NODDY(TBAR3AM) 1485 END IF 1486 IF (LOCDBG) THEN 1487 WRITE(LUPRI,*) 'CCSDT_TBAR31_AM> LISTB,IDLSTB:',LISTB,IDLSTB 1488 WRITE(LUPRI,*) 'CCSDT_TBAR31_AM> FREQB,RHS_ONLY:',FREQB,RHS_ONLY 1489 WRITE(LUPRI,*) 'CCSDT_TBAR31_AM> LISTL,LISTR:',LISTL,LISTR 1490 WRITE(LUPRI,*) 'CCSDT_TBAR31_AM> norm^2(TBAR3AM):', 1491 & DDOT(NT1AMX**3,TBAR3AM,1,TBAR3AM,1) 1492 END IF 1493 1494 CALL QEXIT('CCTBAR31') 1495 1496 RETURN 1497 END 1498*---------------------------------------------------------------------* 1499* END OF SUBROUTINE CCSDT_TBAR31_NODDY * 1500*---------------------------------------------------------------------* 1501*=====================================================================* 1502 SUBROUTINE CCSDT_FBC_NODDY(NEW_NODDY_F_CONT, 1503 & IBTRAN,IDOTS,DOTPROD,NBTRAN,MXVEC, 1504 & LISTL,LISTB,LISTC,WORK,LWORK) 1505*---------------------------------------------------------------------* 1506* 1507* Purpose: compute a triples contribution to F/B matrix contraction 1508* 1509* BCON = <L2|[[H,R_1],R_3]|HF> 1510* 1511* Written by Christof Haettig, May 2002 1512* 1513*=====================================================================* 1514#if defined (IMPLICIT_NONE) 1515 IMPLICIT NONE 1516#else 1517# include "implicit.h" 1518#endif 1519#include "priunit.h" 1520#include "ccsdinp.h" 1521#include "maxorb.h" 1522#include "ccsdsym.h" 1523#include "ccfield.h" 1524#include "ccorb.h" 1525#include "ccr1rsp.h" 1526#include "ccexci.h" 1527#include "dummy.h" 1528#include "inftap.h" 1529 1530 LOGICAL LOCDBG, XI_ONLY 1531 PARAMETER (LOCDBG=.FALSE., XI_ONLY=.FALSE.) 1532 1533 INTEGER ISYM0 1534 PARAMETER (ISYM0 = 1) 1535 1536 LOGICAL NEW_NODDY_F_CONT 1537 CHARACTER*3 LISTL, LISTB, LISTC 1538 INTEGER LWORK, NBTRAN, MXVEC 1539 INTEGER IDOTS(MXVEC,NBTRAN), IBTRAN(3,NBTRAN) 1540#if defined (SYS_CRAY) 1541 REAL DOTPROD(MXVEC,NBTRAN), DDOT 1542 REAL WORK(LWORK) 1543 REAL FREQB, TCON 1544 REAL FREQLST 1545#else 1546 DOUBLE PRECISION DOTPROD(MXVEC,NBTRAN), DDOT 1547 DOUBLE PRECISION WORK(LWORK) 1548 DOUBLE PRECISION FREQB, TCON 1549 DOUBLE PRECISION FREQLST 1550#endif 1551 1552 CHARACTER MODEL*(10) 1553 INTEGER KEND1, ITRAN, IDLSTB, IDLSTC, ISYMRB, KSCR1, 1554 & KFOCKD, KFOCK0, KT1AMP0, KLAMP0, KLAMH0, KFOCKB, 1555 & KLAMPB, KLAMHB, KOME1, KLAM2, KINT1T0, KINT2T0, KDUM, 1556 & KXIAJB, KYIAJB, KINT1SB, KINT2SB, LWRK1, IOPT, KTB3AM, 1557 & KEND2, LWRK2, IVEC, IDLSTL, ISYML, ISYMRC, KL2AM, 1558 & LUFOCK, KRC1AM, KINT1S0, KINT2S0, KFIELD, KFIELDAO 1559 1560 INTEGER ILSTSYM 1561 1562 IF (.NOT. NEW_NODDY_F_CONT) THEN 1563 RETURN 1564 ENDIF 1565 1566 CALL QENTER('CCSDT_FBC_NODDY') 1567 1568 IF (DIRECT) CALL QUIT('CCSDT_FBC_NODDY: DIRECT NOT IMPLEMENTED') 1569 1570 IF (LOCDBG) WRITE(LUPRI,*) 'entered CCSDT_FBC_NODDY...' 1571 1572*---------------------------------------------------------------------* 1573* Memory allocation: 1574*---------------------------------------------------------------------* 1575 KSCR1 = 1 1576 KFOCKD = KSCR1 + NT1AMX 1577 KEND1 = KFOCKD + NORBT 1578 1579 KFOCK0 = KEND1 1580 KEND1 = KFOCK0 + NORBT*NORBT 1581 1582 IF (NONHF) THEN 1583 KFIELD = KEND1 1584 KFIELDAO = KFIELD + NBAST*NBAST 1585 KEND1 = KFIELDAO + NBAST*NBAST 1586 END IF 1587 1588 KT1AMP0 = KEND1 1589 KLAMP0 = KT1AMP0 + NT1AMX 1590 KLAMH0 = KLAMP0 + NLAMDT 1591 KEND1 = KLAMH0 + NLAMDT 1592 1593 KFOCKB = KEND1 1594 KLAMPB = KFOCKB + NORBT*NORBT 1595 KLAMHB = KLAMPB + NLAMDT 1596 KEND1 = KLAMHB + NLAMDT 1597 1598 KOME1 = KEND1 1599 KRC1AM = KOME1 + NT1AMX 1600 KEND1 = KRC1AM + NT1AMX 1601 1602 KL2AM = KEND1 1603 KEND1 = KL2AM + NT1AMX*NT1AMX 1604 1605 KINT1T0 = KEND1 1606 KINT2T0 = KINT1T0 + NT1AMX*NVIRT*NVIRT 1607 KEND1 = KINT2T0 + NRHFT*NRHFT*NT1AMX 1608 1609 KXIAJB = KEND1 1610 KYIAJB = KXIAJB + NT1AMX*NT1AMX 1611 KEND1 = KYIAJB + NT1AMX*NT1AMX 1612 1613 KINT1S0 = KEND1 1614 KINT2S0 = KINT1S0 + NT1AMX*NVIRT*NVIRT 1615 KEND1 = KINT2S0 + NRHFT*NRHFT*NT1AMX 1616 1617 KINT1SB = KEND1 1618 KINT2SB = KINT1SB + NT1AMX*NVIRT*NVIRT 1619 KEND1 = KINT2SB + NRHFT*NRHFT*NT1AMX 1620 1621 LWRK1 = LWORK - KEND1 1622 IF (LWRK1 .LT. 0) THEN 1623 CALL QUIT('Insufficient space in CCSDT_FBC_NODDY') 1624 ENDIF 1625 1626*---------------------------------------------------------------------* 1627* Get zeroth-order Lambda matrices: 1628*---------------------------------------------------------------------* 1629 IOPT = 1 1630 KDUM = KEND1 1631 Call CC_RDRSP('R0',0,ISYM0,IOPT,MODEL,WORK(KT1AMP0),WORK(KDUM)) 1632 1633 Call LAMMAT(WORK(KLAMP0),WORK(KLAMH0),WORK(KT1AMP0), 1634 & WORK(KEND1),LWRK1) 1635 1636*---------------------------------------------------------------------* 1637* Read precalculated integrals from file: 1638*---------------------------------------------------------------------* 1639 CALL CCSDT_READ_NODDY(.TRUE.,WORK(KFOCKD),WORK(KFOCK0), 1640 & WORK(KFIELD),WORK(KFIELDAO), 1641 & .TRUE.,WORK(KXIAJB),WORK(KYIAJB), 1642 & .TRUE.,WORK(KINT1S0),WORK(KINT2S0), 1643 & .TRUE.,WORK(KINT1T0),WORK(KINT2T0), 1644 & .FALSE.,DUMMY,DUMMY,DUMMY,DUMMY, 1645 & NORBT,NLAMDT,NRHFT,NVIRT,NT1AMX) 1646 1647*---------------------------------------------------------------------* 1648* Loop right amplitude vectors: 1649*---------------------------------------------------------------------* 1650 DO ITRAN = 1, NBTRAN 1651 IDLSTB = IBTRAN(1,ITRAN) 1652 IDLSTC = IBTRAN(2,ITRAN) 1653 1654 ISYMRB = ILSTSYM(LISTB,IDLSTB) 1655 ISYMRC = ILSTSYM(LISTC,IDLSTC) 1656 1657 FREQB = FREQLST(LISTB,IDLSTB) 1658 1659 IF (LOCDBG) THEN 1660 WRITE(LUPRI,*) 'IDLSTB:',IDLSTB 1661 WRITE(LUPRI,*) 'FREQB:',FREQB 1662 WRITE(LUPRI,*) 'ISYMRB:',ISYMRB 1663 WRITE(LUPRI,*) 'IDLSTC:',IDLSTC 1664 WRITE(LUPRI,*) 'ISYMRC:',ISYMRC 1665 END IF 1666 1667*---------------------------------------------------------------------* 1668* Compute response amplitudes: 1669*---------------------------------------------------------------------* 1670 KTB3AM = KEND1 1671 KEND2 = KTB3AM + NT1AMX*NT1AMX*NT1AMX 1672 1673 LWRK2 = LWORK - KEND2 1674 IF (LWRK2 .LT. 0) THEN 1675 CALL QUIT('Insufficient space in CCSDT_FBC_NODDY') 1676 ENDIF 1677 1678 IF (LISTB(1:3).EQ.'R1 ' .OR. LISTB(1:3).EQ.'RE ' .OR. 1679 & LISTB(1:3).EQ.'RC ' ) THEN 1680 KDUM = KEND2 1681 CALL CCSDT_T31_NODDY(WORK(KTB3AM),LISTB,IDLSTB,FREQB,XI_ONLY, 1682 & .FALSE.,WORK(KINT1S0),WORK(KINT2S0), 1683 & .FALSE.,WORK(KINT1T0),WORK(KINT2T0), 1684 & .FALSE.,WORK(KXIAJB),WORK(KYIAJB), 1685 & WORK(KINT1SB),WORK(KINT2SB), 1686 & WORK(KLAMPB),WORK(KLAMHB),WORK(KFOCKB), 1687 & WORK(KLAMP0),WORK(KLAMH0),WORK(KFOCK0), 1688 & WORK(KDUM),WORK(KFOCKD), 1689 & WORK(KEND2),LWRK2) 1690 ELSE IF (LISTB(1:3).EQ.'R2 ' .OR. LISTB(1:3).EQ.'ER1') THEN 1691 CALL CCSDT_T32_NODDY(WORK(KTB3AM),LISTB,IDLSTB,FREQB, 1692 & WORK(KINT1S0),WORK(KINT2S0), 1693 & WORK(KLAMP0),WORK(KLAMH0),WORK(KFOCK0), 1694 & WORK(KFOCKD),DUMMY,DUMMY, 1695 & WORK(KSCR1),WORK(KEND2),LWRK2) 1696 ELSE 1697 CALL QUIT('Unknown or illegal list in CCSDT_FBC_NODDY.') 1698 END IF 1699 1700*---------------------------------------------------------------------* 1701* Compute contribution from <L_2|[[H,T^B_3],\tau_nu_1|HF>: 1702*---------------------------------------------------------------------* 1703 IVEC = 1 1704 DO WHILE(IDOTS(IVEC,ITRAN).NE.0 .AND. IVEC.LE.MXVEC) 1705 IDLSTL = IDOTS(IVEC,ITRAN) 1706 1707 IF (LWRK2 .LT. NT2AMX) THEN 1708 CALL QUIT('Insufficient space in CCSDT_FBC_NODDY') 1709 ENDIF 1710 1711 ISYML = ILSTSYM(LISTL,IDLSTL) 1712 1713 ! read left amplitudes from file and square up the doubles part 1714 IOPT = 2 1715 Call CC_RDRSP(LISTL,IDLSTL,ISYML,IOPT,MODEL,DUMMY,WORK(KEND2)) 1716 CALL CC_T2SQ(WORK(KEND2),WORK(KL2AM),ISYML) 1717 1718 CALL DZERO(WORK(KOME1),NT1AMX) 1719 1720 CALL T3_LEFT1(WORK(KOME1),WORK(KYIAJB),WORK(KXIAJB), 1721 & WORK(KL2AM),WORK(KTB3AM)) 1722 1723 CALL T3_LEFT2(WORK(KOME1),WORK(KYIAJB),WORK(KXIAJB), 1724 & WORK(KL2AM),WORK(KTB3AM)) 1725 1726 CALL T3_LEFT3(WORK(KOME1),WORK(KXIAJB),WORK(KL2AM),WORK(KTB3AM)) 1727 1728 ! read RC_1 amplitudes from file 1729 ISYMRC = ILSTSYM(LISTC,IDLSTC) 1730 IOPT = 1 1731 Call CC_RDRSP(LISTC,IDLSTC,ISYMRC,IOPT,MODEL,WORK(KRC1AM),DUMMY) 1732 1733 IF (ISYML.NE.MULD2H(ISYMRB,ISYMRC)) 1734 & CALL QUIT('Symmetry mismatch in CCSDT_FBC_NODDY.') 1735 1736 TCON = DDOT(NT1AM(ISYMRC),WORK(KRC1AM),1,WORK(KOME1),1) 1737 DOTPROD(IVEC,ITRAN) = DOTPROD(IVEC,ITRAN) + TCON 1738 1739 IF (LOCDBG) THEN 1740 WRITE(LUPRI,*) 'ITRAN,IVEC,IDLSTL:',ITRAN,IVEC,IDLSTL 1741 WRITE(LUPRI,*) 'TCON:',TCON 1742 END IF 1743 1744 IVEC = IVEC + 1 1745 END DO ! WHILE 1746 1747 END DO ! ITRAN 1748 1749 CALL QEXIT('CCSDT_FBC_NODDY') 1750 1751 RETURN 1752 END 1753*---------------------------------------------------------------------* 1754* END OF SUBROUTINE CCSDT_FBC_NODDY * 1755*---------------------------------------------------------------------* 1756