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_BMAT_NODDY(LISTA,IDLSTA,LISTB,IDLSTB,IOPTRES, 21 & XLAMDP0,XLAMDH0, 22 & OMEGA1,OMEGA2,OMEGA1EFF,OMEGA2EFF, 23 & IDOTS,DOTPROD,LISTDP,ITRAN, 24 & NBTRAN,MXVEC,WORK,LWORK) 25*---------------------------------------------------------------------* 26* 27* Purpose: compute triples contribution to B matrix transformation 28* 29* (B T^A T^B)^eff_1,2 = (B T^A T^B)_1,2(CCSD) + (B T^A T^B)_1,2(T3) 30* - A_1,2;3 (w_3 - w)^1 (B T^A T^B)_3 31* 32* 33* Written by Christof Haettig, April 2002, based on CCSDT_TRIPLE. 34* 35*=====================================================================* 36#if defined (IMPLICIT_NONE) 37 IMPLICIT NONE 38#else 39# include "implicit.h" 40#endif 41#include "dummy.h" 42#include "priunit.h" 43#include "ccsdinp.h" 44#include "maxorb.h" 45#include "ccsdsym.h" 46#include "ccfield.h" 47#include "ccorb.h" 48 49 LOGICAL LOCDBG, PRINT_T3 50 PARAMETER (LOCDBG=.FALSE., PRINT_T3=.FALSE.) 51 INTEGER ISYM0 52 PARAMETER (ISYM0=1) 53 54 CHARACTER*3 LISTDP, LISTA, LISTB 55 INTEGER LWORK, IOPTRES, ITRAN, MXVEC, NBTRAN, IDLSTA, IDLSTB 56 INTEGER IDOTS(MXVEC,NBTRAN) 57 58#if defined (SYS_CRAY) 59 REAL FREQA, FREQB, FREQ, DDOT, FF 60 REAL DOTPROD(MXVEC,NBTRAN), WORK(LWORK), TWO, ONE 61 REAL XLAMDP0(*), XLAMDH0(*), SIGN 62 REAL OMEGA1(*), OMEGA2(*) 63 REAL OMEGA1EFF(*), OMEGA2EFF(*) 64#else 65 DOUBLE PRECISION FREQA, FREQB, FREQ, DDOT, FF 66 DOUBLE PRECISION DOTPROD(MXVEC,NBTRAN), WORK(LWORK), TWO, ONE 67 DOUBLE PRECISION XLAMDP0(*), XLAMDH0(*), SIGN 68 DOUBLE PRECISION OMEGA1(*), OMEGA2(*) 69 DOUBLE PRECISION OMEGA1EFF(*), OMEGA2EFF(*) 70#endif 71 PARAMETER(ONE = 1.0D0, TWO = 2.0D0 ) 72 73 CHARACTER*10 MODEL 74 INTEGER KSCR1, KFOCKD, KEND1, KINT1T0, KINT2T0, KINT1TA, KINT2TA, 75 & KINT1TB, KINT2TB, KINT1SA, KINT2SA, KINT1SB, KINT2SB, 76 & KINT1SAB, KINT2SAB, KTA3AM, KTB3AM, KFOCKA, KFOCKB, 77 & KLAMHB, KLAMPB, KLAMHA, KLAMPA, KT2AM0, KDUM, KXIAJB, 78 & KYIAJB, KFOCK0, KOMEGA2, LWRK1, KXINT, KEND2, IDX, 79 & LWRK2, KB3AM, IRECNR, KFIELDAO, KFLDA1, KFLDB1, KFLDBUF 80 INTEGER IJ, NIJ, LUSIFC, INDEX, ILSTSYM, IOPT, LUFOCK, LUTEMP, 81 & ILLL, IDEL, ISYDIS, ISYMA, ISYMB, IVEC, ISYMC, IDLSTC, 82 & ISYMD, KFIELD, KINT1S0, KINT2S0, KFCKAR, KFCKBR 83 INTEGER KTA1AM, KTB1AM, KEND1A, LWRK1A, KT2AMB, KT2AMA, KT3AM 84 85 INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 86 87 CALL QENTER('CCSDT_BMAT_NODDY') 88 89 IF(DIRECT)CALL QUIT('DIRECT NOT IMPLEMENTED IN CCSDT_BMAT_NODDY') 90 91*---------------------------------------------------------------------* 92* Memory allocation: 93*---------------------------------------------------------------------* 94 KSCR1 = 1 95 KFOCKD = KSCR1 + NT1AMX 96 KFOCK0 = KFOCKD + NORBT 97 KEND1 = KFOCK0 + NORBT*NORBT 98 99 IF (NONHF) THEN 100 KFIELD = KEND1 101 KFIELDAO = KFIELD + NORBT*NORBT 102 KFLDB1 = KFIELDAO + NORBT*NORBT 103 KFLDA1 = KFLDB1 + NORBT*NORBT 104 KFLDBUF = KFLDA1 + NORBT*NORBT 105 KEND1 = KFLDBUF + NORBT*NORBT 106 END IF 107 108 KTA1AM = KEND1 109 KFOCKA = KTA1AM + NT1AMX 110 KFCKAR = KFOCKA + NORBT*NORBT 111 KLAMPA = KFCKAR + NORBT*NORBT 112 KLAMHA = KLAMPA + NLAMDT 113 KEND1 = KLAMHA + NLAMDT 114 115 KTB1AM = KEND1 116 KFOCKB = KTB1AM + NT1AMX 117 KFCKBR = KFOCKB + NORBT*NORBT 118 KLAMPB = KFCKBR + NORBT*NORBT 119 KLAMHB = KLAMPB + NLAMDT 120 KEND1 = KLAMHB + NLAMDT 121 122 KXIAJB = KEND1 123 KYIAJB = KXIAJB + NT1AMX*NT1AMX 124 KEND1 = KYIAJB + NT1AMX*NT1AMX 125 126 KINT1T0 = KEND1 127 KINT2T0 = KINT1T0 + NT1AMX*NVIRT*NVIRT 128 KEND1 = KINT2T0 + NRHFT*NRHFT*NT1AMX 129 130 KINT1S0 = KEND1 131 KINT2S0 = KINT1S0 + NT1AMX*NVIRT*NVIRT 132 KEND1 = KINT2S0 + NRHFT*NRHFT*NT1AMX 133 134 KB3AM = KEND1 135 KEND1 = KB3AM + NT1AMX*NT1AMX*NT1AMX 136 137 ! what is above has to be kept until the end... 138 ! everything below might be overwritten 139 ! in CC_RHPART_NODDY / CCDOTRSP_NODDY 140 KEND1A = KEND1 141 LWRK1A = LWORK - KEND1A 142 143 KINT1TA = KEND1 144 KINT2TA = KINT1TA + NT1AMX*NVIRT*NVIRT 145 KEND1 = KINT2TA + NRHFT*NRHFT*NT1AMX 146 147 KINT1TB = KEND1 148 KINT2TB = KINT1TB + NT1AMX*NVIRT*NVIRT 149 KEND1 = KINT2TB + NRHFT*NRHFT*NT1AMX 150 151 KINT1SA = KEND1 152 KINT2SA = KINT1SA + NT1AMX*NVIRT*NVIRT 153 KEND1 = KINT2SA + NRHFT*NRHFT*NT1AMX 154 155 KINT1SB = KEND1 156 KINT2SB = KINT1SB + NT1AMX*NVIRT*NVIRT 157 KEND1 = KINT2SB + NRHFT*NRHFT*NT1AMX 158 159 KINT1SAB = KEND1 160 KINT2SAB = KINT1SAB + NT1AMX*NVIRT*NVIRT 161 KEND1 = KINT2SAB + NRHFT*NRHFT*NT1AMX 162 163 KT3AM = KEND1 164 KEND1 = KT3AM + NT1AMX*NT1AMX*NT1AMX 165 166 KOMEGA2 = KEND1 167 KEND1 = KOMEGA2 + NT1AMX*NT1AMX 168 169 KT2AMB = KEND1 170 KT2AMA = KT2AMB + NT1AMX*NT1AMX 171 KEND1 = KT2AMA + NT1AMX*NT1AMX 172 173 KT2AM0 = KOMEGA2 174 KTA3AM = KB3AM 175 KTB3AM = KT3AM 176 177 LWRK1 = LWORK - KEND1 178 IF (LWRK1 .LT. 0) THEN 179 CALL QUIT('Insufficient space in CCSDT_BMAT_NODDY') 180 ENDIF 181 182*---------------------------------------------------------------------* 183* Read SCF orbital energies from file: 184*---------------------------------------------------------------------* 185 LUSIFC = -1 186 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY, 187 & .FALSE.) 188 REWIND LUSIFC 189 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 190 READ (LUSIFC) 191 READ (LUSIFC) (WORK(KFOCKD+I-1), I=1,NORBT) 192 CALL GPCLOSE(LUSIFC,'KEEP') 193 194*---------------------------------------------------------------------* 195* read zeroth-order AO Fock matrix from file: 196*---------------------------------------------------------------------* 197 LUFOCK = -1 198 CALL GPOPEN(LUFOCK,'CC_FCKH','OLD',' ','UNFORMATTED',IDUMMY, 199 & .FALSE.) 200 REWIND(LUFOCK) 201 READ(LUFOCK) (WORK(KFOCK0-1+I),I=1,N2BST(ISYM0)) 202 CALL GPCLOSE(LUFOCK,'KEEP') 203 204 CALL CC_FCKMO(WORK(KFOCK0),XLAMDP0,XLAMDH0, 205 & WORK(KEND1),LWRK1,ISYM0,ISYM0,ISYM0) 206 207*---------------------------------------------------------------------* 208* If needed get external field: 209*---------------------------------------------------------------------* 210 IF (NONHF .AND. NFIELD.GT.0) THEN 211 CALL DZERO(WORK(KFIELDAO),NORBT*NORBT) 212 DO I = 1, NFIELD 213 FF = EFIELD(I) 214 CALL CC_ONEP(WORK(KFIELDAO),WORK(KEND1),LWRK1,FF,1,LFIELD(I)) 215 ENDDO 216 217 CALL DCOPY(NORBT*NORBT,WORK(KFIELDAO),1,WORK(KFIELD),1) 218 CALL CC_FCKMO(WORK(KFIELD),XLAMDP0,XLAMDH0, 219 * WORK(KEND1),LWRK1,1,1,1) 220 ENDIF 221 222*---------------------------------------------------------------------* 223* Compute some integrals: 224* XINT1S0 = (CK|BD) 225* XINT2S0 = (CK|LJ) 226* XINT1T0 = (KC|BD) 227* XINT2T0 = (KC|LJ) 228* XIAJB = 2(IA|JB) - (IB|JA) 229* YIAJB = (IA|JB) 230*---------------------------------------------------------------------* 231 CALL CCSDT_INTS0_NODDY(.TRUE.,WORK(KXIAJB), WORK(KYIAJB), 232 & .TRUE.,WORK(KINT1S0),WORK(KINT2S0), 233 & .TRUE.,WORK(KINT1T0),WORK(KINT2T0), 234 & XLAMDP0,XLAMDH0, 235 & WORK(KEND1),LWRK1) 236 237*---------------------------------------------------------------------* 238* Compute corrections to triples vector T^A_3, and corresponding 239* lambda matrices and the XINT1SA,XINT2SA integrals and set FREQA: 240*---------------------------------------------------------------------* 241 IF (LISTA(1:3).EQ.'R1 ' .or. LISTA(1:3).EQ.'RE ' .or. 242 & LISTA(1:3).EQ.'RC ' ) THEN 243 244 KDUM = KEND1 245 CALL CCSDT_T31_NODDY(WORK(KTA3AM),LISTA,IDLSTA,FREQA,.FALSE., 246 & .FALSE.,WORK(KINT1S0),WORK(KINT2S0), 247 & .FALSE.,WORK(KINT1T0),WORK(KINT2T0), 248 & .FALSE.,WORK(KXIAJB), WORK(KYIAJB), 249 & WORK(KINT1SA),WORK(KINT2SA), 250 & WORK(KLAMPA),WORK(KLAMHA),WORK(KFOCKA), 251 & XLAMDP0,XLAMDH0,WORK(KFOCK0), 252 & WORK(KDUM),WORK(KFOCKD), 253 & WORK(KEND1),LWRK1) 254 255 ELSE IF (LISTA(1:3).EQ.'R2 ' .or. LISTA(1:3).EQ.'ER1') THEN 256 257 CALL CCSDT_T32_NODDY(WORK(KTA3AM),LISTA,IDLSTA,FREQA, 258 & WORK(KINT1S0),WORK(KINT2S0), 259 & XLAMDP0,XLAMDH0,WORK(KFOCK0), 260 & WORK(KFOCKD),WORK(KFIELDAO),WORK(KFIELD), 261 & WORK(KSCR1),WORK(KEND1),LWRK1) 262 263 CALL CCLR_LAMTRA(XLAMDP0,WORK(KLAMPA),XLAMDH0,WORK(KLAMHA), 264 & WORK(KTA1AM),ISYMA) 265 266 CALL CCSDT_INTS1_NODDY(.TRUE.,WORK(KINT1SA),WORK(KINT2SA), 267 & .FALSE.,DUMMY,DUMMY, 268 & XLAMDP0,XLAMDH0, 269 & WORK(KLAMPA),WORK(KLAMHA), 270 & WORK(KEND1),LWRK1) 271 272 ELSE 273 CALL QUIT('Unknown LISTB in CCSDT_BMAT_NODDY.') 274 END IF 275 276 CALL DSCAL(NT1AMX*NT1AMX*NT1AMX,-1.0D0,WORK(KTA3AM),1) 277 278 IF (NONHF) THEN 279 LUTEMP = -1 280 CALL GPOPEN(LUTEMP,'T3AMPA','UNKNOWN',' ','UNFORMATTED', 281 & IDUMMY,.FALSE.) 282 REWIND LUTEMP 283 WRITE (LUTEMP) (WORK(KTA3AM-1+IDX),IDX=1,NT1AMX*NT1AMX*NT1AMX) 284 CALL GPCLOSE(LUTEMP,'KEEP') 285 END IF 286 287*---------------------------------------------------------------------* 288* Compute corrections to triples vector T^B_3, and corresponding 289* lambda matrices and the XINT1SB,XINT2SB integrals and set FREQB: 290*---------------------------------------------------------------------* 291 IF (LISTB(1:3).EQ.'R1 ' .or. LISTB(1:3).EQ.'RE ' .or. 292 & LISTB(1:3).EQ.'RC ' ) THEN 293 294 KDUM = KEND1 295 CALL CCSDT_T31_NODDY(WORK(KTB3AM),LISTB,IDLSTB,FREQB,.FALSE., 296 & .FALSE.,WORK(KINT1S0),WORK(KINT2S0), 297 & .FALSE.,WORK(KINT1T0),WORK(KINT1T0), 298 & .FALSE.,WORK(KXIAJB), WORK(KYIAJB), 299 & WORK(KINT1SB),WORK(KINT2SB), 300 & WORK(KLAMPB),WORK(KLAMHB),WORK(KFOCKB), 301 & XLAMDP0,XLAMDH0,WORK(KFOCK0), 302 & WORK(KDUM),WORK(KFOCKD), 303 & WORK(KEND1),LWRK1) 304 305 ELSE IF (LISTB(1:3).EQ.'R2 ' .or. LISTB(1:3).EQ.'ER1') THEN 306 307 CALL CCSDT_T32_NODDY(WORK(KTB3AM),LISTB,IDLSTB,FREQB, 308 & WORK(KINT1S0),WORK(KINT2S0), 309 & XLAMDP0,XLAMDH0,WORK(KFOCK0), 310 & WORK(KFOCKD),WORK(KFIELDAO),WORK(KFIELD), 311 & WORK(KSCR1),WORK(KEND1),LWRK1) 312 313 CALL CCLR_LAMTRA(XLAMDP0,WORK(KLAMPB),XLAMDH0,WORK(KLAMHB), 314 & WORK(KTB1AM),ISYMB) 315 316 CALL CCSDT_INTS1_NODDY(.TRUE.,WORK(KINT1SB),WORK(KINT2SB), 317 & .FALSE.,DUMMY,DUMMY, 318 & XLAMDP0,XLAMDH0, 319 & WORK(KLAMPB),WORK(KLAMHB), 320 & WORK(KEND1),LWRK1) 321 322 ELSE 323 CALL QUIT('Unknown LISTB in CCSDT_BMAT_NODDY.') 324 END IF 325 326 CALL DSCAL(NT1AMX*NT1AMX*NT1AMX,-1.0D0,WORK(KTB3AM),1) 327 328 IF (NONHF) THEN 329 LUTEMP = -1 330 CALL GPOPEN(LUTEMP,'T3AMPB','UNKNOWN',' ','UNFORMATTED', 331 & IDUMMY,.FALSE.) 332 REWIND LUTEMP 333 WRITE (LUTEMP) (WORK(KTB3AM-1+IDX),IDX=1,NT1AMX*NT1AMX*NT1AMX) 334 CALL GPCLOSE(LUTEMP,'KEEP') 335 END IF 336 337*---------------------------------------------------------------------* 338* Compute required integrals: 339*---------------------------------------------------------------------* 340 CALL DZERO(WORK(KINT1SAB),NT1AMX*NVIRT*NVIRT) 341 CALL DZERO(WORK(KINT2SAB),NT1AMX*NRHFT*NRHFT) 342 343 CALL DZERO(WORK(KINT1TA),NT1AMX*NVIRT*NVIRT) 344 CALL DZERO(WORK(KINT2TA),NT1AMX*NRHFT*NRHFT) 345 346 CALL DZERO(WORK(KINT1TB),NT1AMX*NVIRT*NVIRT) 347 CALL DZERO(WORK(KINT2TB),NT1AMX*NRHFT*NRHFT) 348 349 DO ISYMD = 1, NSYM 350 DO ILLL = 1,NBAS(ISYMD) 351 IDEL = IBAS(ISYMD) + ILLL 352 ISYDIS = MULD2H(ISYMD,ISYMOP) 353 354C ---------------------------- 355C Work space allocation no. 2. 356C ---------------------------- 357 KXINT = KEND1 358 KEND2 = KXINT + NDISAO(ISYDIS) 359 LWRK2 = LWORK - KEND2 360 IF (LWRK2 .LT. 0) THEN 361 WRITE(LUPRI,*) 'Need : ',KEND2,'Available : ',LWORK 362 CALL QUIT('Insufficient space in CCSDT_BMAT_NODDY') 363 ENDIF 364 365C --------------------------- 366C Read in batch of integrals. 367C --------------------------- 368 CALL CCRDAO(WORK(KXINT),IDEL,1,WORK(KEND2),LWRK2, 369 * IRECNR,DIRECT) 370 371C ---------------------------------- 372C Calculate integrals needed in CC3: 373C ---------------------------------- 374 CALL CCSDT_TRAN1_R(WORK(KINT1TA),WORK(KINT2TA), 375 & XLAMDP0,XLAMDH0, 376 & WORK(KLAMPA),WORK(KLAMHA), 377 & WORK(KXINT),IDEL) 378 379 CALL CCSDT_TRAN1_R(WORK(KINT1TB),WORK(KINT2TB), 380 & XLAMDP0,XLAMDH0, 381 & WORK(KLAMPB),WORK(KLAMHB), 382 & WORK(KXINT),IDEL) 383 384 385 ! XINT1SAB = XINT1SAB + (C-barB K-barA|B D) 386 ! XINT2SAB = XINT2SAB + (C-barB K-barA|L J) 387 CALL CCSDT_TRAN3_R(WORK(KINT1SAB),WORK(KINT2SAB), 388 & XLAMDP0,XLAMDH0, 389 & WORK(KLAMPB),WORK(KLAMHA), 390 & XLAMDP0,XLAMDH0, 391 & WORK(KXINT),IDEL) 392 393 ! XINT1SAB = XINT1SAB + (C-barA K-barB|B D) 394 ! XINT2SAB = XINT2SAB + (C-barA K-barB|L J) 395 CALL CCSDT_TRAN3_R(WORK(KINT1SAB),WORK(KINT2SAB), 396 & XLAMDP0,XLAMDH0, 397 & WORK(KLAMPA),WORK(KLAMHB), 398 & XLAMDP0,XLAMDH0, 399 & WORK(KXINT),IDEL) 400 401 ! XINT1SAB = XINT1SAB + (C-barB K|B-barA D) 402 ! XINT2SAB = XINT2SAB + (C-barB K|L J-barA) 403 CALL CCSDT_TRAN3_R(WORK(KINT1SAB),WORK(KINT2SAB), 404 & XLAMDP0,XLAMDH0, 405 & WORK(KLAMPB),XLAMDH0, 406 & WORK(KLAMPA),WORK(KLAMHA), 407 & WORK(KXINT),IDEL) 408 409 ! XINT1SAB = XINT1SAB + (C-barA K|B-barB D) 410 ! XINT2SAB = XINT2SAB + (C-barA K|L J-barB) 411 CALL CCSDT_TRAN3_R(WORK(KINT1SAB),WORK(KINT2SAB), 412 & XLAMDP0,XLAMDH0, 413 & WORK(KLAMPA),XLAMDH0, 414 & WORK(KLAMPB),WORK(KLAMHB), 415 & WORK(KXINT),IDEL) 416 417 ! XINT1SAB = XINT1SAB + (C K-barB|B-barA D) 418 ! XINT2SAB = XINT2SAB + (C K-barB|L J-barA) 419 CALL CCSDT_TRAN3_R(WORK(KINT1SAB),WORK(KINT2SAB), 420 & XLAMDP0,XLAMDH0, 421 & XLAMDP0,WORK(KLAMHB), 422 & WORK(KLAMPA),WORK(KLAMHA), 423 & WORK(KXINT),IDEL) 424 425 ! XINT1SAB = XINT1SAB + (C K-barA|B-barB D) 426 ! XINT2SAB = XINT2SAB + (C K-barA|L J-barB) 427 CALL CCSDT_TRAN3_R(WORK(KINT1SAB),WORK(KINT2SAB), 428 & XLAMDP0,XLAMDH0, 429 & XLAMDP0,WORK(KLAMHA), 430 & WORK(KLAMPB),WORK(KLAMHB), 431 & WORK(KXINT),IDEL) 432 433 END DO 434 END DO 435 436C CALL CCSDT_INTS2_NODDY(WORK(KXINT1SAB),WORK(KXINT2SAB), 437C & XLAMDP0,XLAMDH0, 438C & WORK(KLAMPB),WORK(KLAMHB), 439C & WORK(KLAMPA),WORK(KLAMHA), 440C & WORK(KEND1),LWRK1) 441 442C CALL CCSDT_INTS1_NODDY(.FALSE.,WORK(KINT1SA),WORK(KINT2SA), 443C & .TRUE.,WORK(KINT1TA),WORK(KINT2TA), 444C & XLAMDP0,XLAMDH0, 445C & WORK(KLAMPA),WORK(KLAMHA), 446C & WORK(KEND1),LWRK1) 447 448C CALL CCSDT_INTS1_NODDY(.FALSE.,WORK(KINT1SB),WORK(KINT2SB), 449C & .TRUE.,WORK(KINT1TB),WORK(KINT2TB), 450C & XLAMDP0,XLAMDH0, 451C & WORK(KLAMPB),WORK(KLAMHB), 452C & WORK(KEND1),LWRK1) 453 454*---------------------------------------------------------------------* 455* Compute corrections to doubles result vector from T^A_3, T^B_3: 456* omega2 += P^AB <mu_2|[[H,T^A_1],T^B_3]|HF> 457*---------------------------------------------------------------------* 458 ! read RA_1 amplitudes from file 459 ISYMA = ILSTSYM(LISTA,IDLSTA) 460 IOPT = 1 461 Call CC_RDRSP(LISTA,IDLSTA,ISYMA,IOPT,MODEL,WORK(KTA1AM),DUMMY) 462 463 ! read RB_1 amplitudes from file 464 ISYMB = ILSTSYM(LISTB,IDLSTB) 465 IOPT = 1 466 Call CC_RDRSP(LISTB,IDLSTB,ISYMB,IOPT,MODEL,WORK(KTB1AM),DUMMY) 467 468 ! compute fock matrix like contribution to [H,T^A_1] 469 CALL DZERO(WORK(KFCKAR),NORBT*NORBT) 470 CALL CCSDT_FCK_R(WORK(KFCKAR),WORK(KXIAJB),WORK(KTA1AM)) 471 472 ! compute fock matrix like contribution to [H,T^B_1] 473 CALL DZERO(WORK(KFCKBR),NORBT*NORBT) 474 CALL CCSDT_FCK_R(WORK(KFCKBR),WORK(KXIAJB),WORK(KTB1AM)) 475 476 477 ! <mu_2|[[H,T^A_1],T^B_3]|HF> 478 CALL DZERO(WORK(KOMEGA2),NT1AMX*NT1AMX) 479 CALL CCSDT_OMEGA2(WORK(KOMEGA2),WORK(KINT1TA),WORK(KINT2TA), 480 & WORK(KTB3AM),WORK(KFCKAR)) 481 482 DO I = 1,NT1AMX 483 DO J = 1,I 484 IJ = NT1AMX*(I-1) + J 485 NIJ = INDEX(I,J) 486 OMEGA2(NIJ) = OMEGA2(NIJ) + WORK(KOMEGA2+IJ-1) 487 END DO 488 END DO 489 490 491 ! <mu_2|[[H,T^B_1],T^A_3]|HF> 492 CALL DZERO(WORK(KOMEGA2),NT1AMX*NT1AMX) 493 494 CALL CCSDT_OMEGA2(WORK(KOMEGA2),WORK(KINT1TB),WORK(KINT2TB), 495 & WORK(KTA3AM),WORK(KFCKBR)) 496 497 DO I = 1,NT1AMX 498 DO J = 1,I 499 IJ = NT1AMX*(I-1) + J 500 NIJ = INDEX(I,J) 501 OMEGA2(NIJ) = OMEGA2(NIJ) + WORK(KOMEGA2+IJ-1) 502 END DO 503 END DO 504 505*---------------------------------------------------------------------* 506* Compute triples vector (B T^A T^B)_3: 507*---------------------------------------------------------------------* 508 509 CALL DZERO(WORK(KB3AM),NT1AMX*NT1AMX*NT1AMX) 510 511 if (.true. ) then 512 513 IF (LWRK1 .LT. NT2AMX) THEN 514 CALL QUIT('Insufficient space in CCSDT_BMAT_NODDY') 515 ENDIF 516 517 ! add <nu_3|[H^AB,T^0_2]|HF> 518 IOPT = 2 519 KDUM = KEND1 520 CALL CC_RDRSP('R0',0,ISYM0,IOPT,MODEL,WORK(KDUM),WORK(KEND1)) 521 CALL CC_T2SQ(WORK(KEND1),WORK(KT2AM0),ISYM0) 522 523 CALL CCSDT_T3AM_R(WORK(KB3AM),0.0D0, 524 & WORK(KINT1SAB),WORK(KINT2SAB),WORK(KT2AM0), 525 & WORK(KSCR1),WORK(KFOCKD), 526 & .FALSE.,DUMMY,.FALSE.) 527 528 529 ! add <nu_3|[H^A,T^B_2]|HF> 530 ISYMB = ILSTSYM(LISTB,IDLSTB) 531 IOPT = 2 532 KDUM = KEND1 533 CALL CC_RDRSP(LISTB,IDLSTB,ISYMB,IOPT,MODEL, 534 & WORK(KDUM),WORK(KEND1)) 535 CALL CCLR_DIASCL(WORK(KEND1),TWO,ISYMB) 536 CALL CC_T2SQ(WORK(KEND1),WORK(KT2AMB),ISYMB) 537 538 CALL CCSDT_T3AM_R(WORK(KB3AM),0.0D0, 539 & WORK(KINT1SA),WORK(KINT2SA),WORK(KT2AMB), 540 & WORK(KSCR1),WORK(KFOCKD), 541 & .FALSE.,DUMMY,.FALSE.) 542 543 ! add <nu_3|[H^B,T^A_2]|HF> 544 ISYMA = ILSTSYM(LISTA,IDLSTA) 545 IOPT = 2 546 KDUM = KEND1 547 CALL CC_RDRSP(LISTA,IDLSTA,ISYMA,IOPT,MODEL, 548 & WORK(KDUM),WORK(KEND1)) 549 CALL CCLR_DIASCL(WORK(KEND1),TWO,ISYMA) 550 CALL CC_T2SQ(WORK(KEND1),WORK(KT2AMA),ISYMA) 551 552 CALL CCSDT_T3AM_R(WORK(KB3AM),0.0D0, 553 & WORK(KINT1SB),WORK(KINT2SB),WORK(KT2AMA), 554 & WORK(KSCR1),WORK(KFOCKD), 555 & .FALSE.,DUMMY,.FALSE.) 556 557 IF (.TRUE. .AND. NONHF .AND. NFIELD.GT.0) THEN 558cch 559 write(lupri,*) 'norm^2(field) before FF:', 560 & ddot(norbt*norbt,work(kfield),1,work(kfield),1) 561 write(lupri,*) 'norm^2(b3am) before FF:', 562 & ddot(nt1amx**3,work(kb3am),1,work(kb3am),1) 563cch 564 CALL CCSDT_XKSI3_1(WORK(KB3AM),WORK(KFIELD), 565 & WORK(KT2AMB),WORK(KT2AMA),ONE) 566 CALL CCSDT_XKSI3_1(WORK(KB3AM),WORK(KFIELD), 567 & WORK(KT2AMA),WORK(KT2AMB),ONE) 568cch 569 write(lupri,*) 'norm^2(b3am) after FF-1:', 570 & ddot(nt1amx**3,work(kb3am),1,work(kb3am),1) 571cch 572 573 ! calculate one-index transf. field integrals: V^A = [V,T1^A] 574 CALL DCOPY(NORBT*NORBT,WORK(KFIELDAO),1,WORK(KFLDBUF),1) 575 CALL DCOPY(NORBT*NORBT,WORK(KFIELDAO),1,WORK(KFLDA1),1) 576 CALL CC_FCKMO(WORK(KFLDBUF),WORK(KLAMPA),XLAMDH0, 577 & WORK(KEND1),LWRK1,ISYMA,ISYM0,ISYMA) 578 CALL CC_FCKMO(WORK(KFLDA1),XLAMDP0,WORK(KLAMHA), 579 & WORK(KEND1),LWRK1,ISYM0,ISYMA,ISYMA) 580 CALL DAXPY(NORBT*NORBT,ONE,WORK(KFLDBUF),1,WORK(KFLDA1),1) 581ctest 582 write(lupri,*) 'norm^2(flda1) before FF:', 583 & ddot(norbt*norbt,work(kflda1),1,work(kflda1),1) 584 585C CALL DSCAL(NORBT*NORBT,-1.0D0,WORK(KFLDA1),1) 586ctest 587 588 LUTEMP = -1 589 CALL GPOPEN(LUTEMP,'T3AMPB','UNKNOWN',' ','UNFORMATTED', 590 & IDUMMY,.FALSE.) 591 REWIND LUTEMP 592 READ (LUTEMP) (WORK(KT3AM-1+IDX),IDX=1,NT1AMX*NT1AMX*NT1AMX) 593 CALL GPCLOSE(LUTEMP,'DELETE') 594 595 CALL CCSDT_XKSI3_2(WORK(KB3AM),WORK(KFLDA1),WORK(KT3AM)) 596 597 ! calculate one-index transf. field integrals: V^B = [V,T1^B] 598 CALL DCOPY(NORBT*NORBT,WORK(KFIELDAO),1,WORK(KFLDBUF),1) 599 CALL DCOPY(NORBT*NORBT,WORK(KFIELDAO),1,WORK(KFLDB1),1) 600 CALL CC_FCKMO(WORK(KFLDBUF),WORK(KLAMPB),XLAMDH0, 601 & WORK(KEND1),LWRK1,ISYMB,ISYM0,ISYMB) 602 CALL CC_FCKMO(WORK(KFLDB1),XLAMDP0,WORK(KLAMHB), 603 & WORK(KEND1),LWRK1,ISYM0,ISYMB,ISYMB) 604 CALL DAXPY(NORBT*NORBT,ONE,WORK(KFLDBUF),1,WORK(KFLDB1),1) 605ctest 606 write(lupri,*) 'norm^2(fldb1) before FF:', 607 & ddot(norbt*norbt,work(kfldb1),1,work(kfldb1),1) 608 609C CALL DSCAL(NORBT*NORBT,-1.0D0,WORK(KFLDB1),1) 610ctest 611 612 LUTEMP = -1 613 CALL GPOPEN(LUTEMP,'T3AMPA','UNKNOWN',' ','UNFORMATTED', 614 & IDUMMY,.FALSE.) 615 REWIND LUTEMP 616 READ (LUTEMP) (WORK(KT3AM-1+IDX),IDX=1,NT1AMX*NT1AMX*NT1AMX) 617 CALL GPCLOSE(LUTEMP,'DELETE') 618 619 CALL CCSDT_XKSI3_2(WORK(KB3AM),WORK(KFLDB1),WORK(KT3AM)) 620cch 621 write(lupri,*) 'norm^2(b3am) after FF:', 622 & ddot(nt1amx**3,work(kb3am),1,work(kb3am),1) 623cch 624 END IF 625 626 else 627 628 CALL CCSDT_B3AM(WORK(KB3AM), 629 & WORK(KINT1SAB),WORK(KINT2SAB),WORK(KFOCKD), 630 & LISTA,IDLSTA,WORK(KINT1SA),WORK(KINT2SA), 631 & LISTB,IDLSTB,WORK(KINT1SB),WORK(KINT2SB), 632 & WORK(KSCR1),WORK(KT2AM0),WORK(KEND1),LWRK1) 633 634 end if 635 636C if (ioptres.eq.5) then 637C write(lupri,*) 'IOPTRES=5... set B3AM to zero!!!' 638C call dzero(work(kb3am),nt1amx*nt1amx*nt1amx) 639C end if 640 641 if (print_t3) then 642 write(lupri,*)'CCSDT_B_NODDY> vector types:',lista,listb 643 write(lupri,*)'CCSDT_B_NODDY> idlsts:',idlsta,idlstb 644 write(lupri,*)'CCSDT_B_NODDY> freq:',freqa+freqb 645 call print_pt3_noddy(WORK(KB3AM)) 646 end if 647 648*---------------------------------------------------------------------* 649* Now we split: 650* for IOPTRES < 5 we compute an effective rhs vector 651* for IOPTRES = 5 we compute contractions L B T^A T^B 652*---------------------------------------------------------------------* 653 IF (IOPTRES.GE.1 .AND. IOPTRES.LE.4) THEN 654 655 CALL DCOPY(NT1AMX,OMEGA1,1,OMEGA1EFF,1) 656 CALL DCOPY(NT2AMX,OMEGA2,1,OMEGA2EFF,1) 657 658 FREQ = FREQA + FREQB 659 660 CALL CC_RHPART_NODDY(OMEGA1EFF,OMEGA2EFF,WORK(KB3AM),FREQ, 661 & WORK(KFOCKD),WORK(KFOCK0),WORK(KFIELD), 662 & WORK(KXIAJB),WORK(KINT1T0),WORK(KINT2T0), 663 & WORK(KEND1A),LWRK1A) 664 665 ELSE IF (IOPTRES.EQ.5) THEN 666 667 SIGN = -1.0D0 668 CALL CCDOTRSP_NODDY(DUMMY,DUMMY,WORK(KB3AM),SIGN, 669 & ITRAN,LISTDP,IDOTS,DOTPROD,MXVEC, 670 & XLAMDP0,XLAMDH0,WORK(KFOCK0),WORK(KFOCKD), 671 & WORK(KXIAJB), WORK(KYIAJB), 672 & WORK(KINT1T0),WORK(KINT2T0), 673 & WORK(KINT1S0),WORK(KINT2S0), 674 & 'CCSDT_B_NODDY',LOCDBG,.FALSE.,.FALSE., 675 & WORK(KEND1A),LWRK1A) 676 677 ELSE 678 CALL QUIT('Illegal value for IOPTRES IN CCSDT_BMAT_NODDY') 679 END IF 680 681 CALL QEXIT('CCSDT_BMAT_NODDY') 682 RETURN 683 END 684 685*---------------------------------------------------------------------* 686* END OF SUBROUTINE CCSDT_BMAT_NODDY * 687*---------------------------------------------------------------------* 688*=====================================================================* 689 SUBROUTINE CCSDT_B3AM(B3AM,XINT1SAB,XINT2SAB,FOCKD, 690 & LISTA,IDLSTA,XINT1SA,XINT2SA, 691 & LISTB,IDLSTB,XINT1SB,XINT2SB, 692 & SCR1,T2AM,WORK,LWORK) 693*---------------------------------------------------------------------* 694* Purpose: compute triples result of B matrix transformation 695*---------------------------------------------------------------------* 696 IMPLICIT NONE 697#include "dummy.h" 698#include "ccsdsym.h" 699 700 INTEGER ISYM0 701 PARAMETER (ISYM0=1) 702 703 CHARACTER*3 LISTA, LISTB 704 INTEGER LWORK, IDLSTB, IDLSTA 705 706#if defined (SYS_CRAY) 707 REAL B3AM(*), SCR1(*), T2AM(*), WORK(*) 708 REAL XINT1SA(*), XINT2SA(*), XINT1SB(*), XINT2SB(*) 709 REAL XINT1SAB(*), XINT2SAB(*), FOCKD(*) 710 REAL TWO, ZERO 711#else 712 DOUBLE PRECISION B3AM(*), SCR1(*), T2AM(*), WORK(*) 713 DOUBLE PRECISION XINT1SA(*), XINT2SA(*), XINT1SB(*), XINT2SB(*) 714 DOUBLE PRECISION XINT1SAB(*), XINT2SAB(*), FOCKD(*) 715 DOUBLE PRECISION TWO, ZERO 716#endif 717 PARAMETER( TWO = 2.0D0, ZERO = 0.0D0 ) 718 719 CHARACTER*10 MODEL 720 INTEGER IOPT, ISYMA, ISYMB, ILSTSYM 721 722 IF (LWORK .LT. NT2AMX) THEN 723 CALL QUIT('Insufficient space in CCSDT_B3AM') 724 ENDIF 725 726 ! add <nu_3|[H^AB,T^0_2]|HF> 727 IOPT = 2 728 CALL CC_RDRSP('R0',0,ISYM0,IOPT,MODEL,DUMMY,WORK) 729 CALL CC_T2SQ(WORK,T2AM,ISYM0) 730 731 CALL CCSDT_T3AM_R(B3AM,ZERO,XINT1SAB,XINT2SAB,T2AM, 732 & SCR1,FOCKD,.FALSE.,DUMMY,.FALSE.) 733 734 735 ! add <nu_3|[H^A,T^B_2]|HF> 736 ISYMB = ILSTSYM(LISTB,IDLSTB) 737 IOPT = 2 738 CALL CC_RDRSP(LISTB,IDLSTB,ISYMB,IOPT,MODEL,DUMMY,WORK) 739 CALL CCLR_DIASCL(WORK,TWO,ISYMB) 740 CALL CC_T2SQ(WORK,T2AM,ISYMB) 741 742 CALL CCSDT_T3AM_R(B3AM,ZERO,XINT1SA,XINT2SA,T2AM, 743 & SCR1,FOCKD,.FALSE.,DUMMY,.FALSE.) 744 745 ! add <nu_3|[H^B,T^A_2]|HF> 746 ISYMA = ILSTSYM(LISTA,IDLSTA) 747 IOPT = 2 748 CALL CC_RDRSP(LISTA,IDLSTA,ISYMA,IOPT,MODEL,DUMMY,WORK) 749 CALL CCLR_DIASCL(WORK,TWO,ISYMA) 750 CALL CC_T2SQ(WORK,T2AM,ISYMA) 751 752 CALL CCSDT_T3AM_R(B3AM,ZERO,XINT1SB,XINT2SB,T2AM, 753 & SCR1,FOCKD,.FALSE.,DUMMY,.FALSE.) 754 755 RETURN 756 END 757*---------------------------------------------------------------------* 758* END OF SUBROUTINE CCSDT_B3AM 759*---------------------------------------------------------------------* 760*=====================================================================* 761 SUBROUTINE CCSDT_INTS2_NODDY(XINT1SAB,XINT2SAB,XLAMDP0,XLAMDH0, 762 & XLAMDPB,XLAMDHB,XLAMDPA,XLAMDHA, 763 & WORK,LWORK) 764*---------------------------------------------------------------------* 765* Purpose: 766* Loop over distributions of integrals and compute second-order 767* response versions XINT1SAB=(ck|db)-AB and XINT2SAB=(ck|lj)-AB 768*---------------------------------------------------------------------* 769 IMPLICIT NONE 770#include "priunit.h" 771#include "ccsdinp.h" 772#include "maxorb.h" 773#include "ccsdsym.h" 774#include "ccorb.h" 775 776 INTEGER LWORK 777#if defined (SYS_CRAY) 778 REAL WORK(*), XINT1SAB(*), XINT2SAB(*) 779 REAL XLAMDP0(*), XLAMDH0(*) 780 REAL XLAMDPA(*), XLAMDHA(*) 781 REAL XLAMDPB(*), XLAMDHB(*) 782#else 783 DOUBLE PRECISION WORK(*), XINT1SAB(*), XINT2SAB(*) 784 DOUBLE PRECISION XLAMDP0(*), XLAMDH0(*) 785 DOUBLE PRECISION XLAMDPA(*), XLAMDHA(*) 786 DOUBLE PRECISION XLAMDPB(*), XLAMDHB(*) 787#endif 788 789 INTEGER ISYMD, ILLL, IDEL, ISYDIS, KXINT, KEND2, LWRK2, IRECNR 790 791 CALL DZERO(XINT1SAB,NT1AMX*NVIRT*NVIRT) 792 CALL DZERO(XINT2SAB,NT1AMX*NRHFT*NRHFT) 793 794 795 DO ISYMD = 1, NSYM 796 DO ILLL = 1,NBAS(ISYMD) 797 IDEL = IBAS(ISYMD) + ILLL 798 ISYDIS = ISYMD 799 800C ---------------------------- 801C Work space allocation no. 2. 802C ---------------------------- 803 KXINT = 1 804 KEND2 = KXINT + NDISAO(ISYDIS) 805 LWRK2 = LWORK - KEND2 806 IF (LWRK2 .LT. 0) THEN 807 WRITE(LUPRI,*) 'Need : ',KEND2,'Available : ',LWORK 808 CALL QUIT('Insufficient space in CCSDT_INTS2_NODDY') 809 ENDIF 810 811C --------------------------- 812C Read in batch of integrals. 813C --------------------------- 814 CALL CCRDAO(WORK(KXINT),IDEL,1,WORK(KEND2),LWRK2, 815 * IRECNR,DIRECT) 816 817C ---------------------------------- 818C Calculate integrals needed in CC3: 819C ---------------------------------- 820 ! XINT1SAB = XINT1SAB + (C-barB K-barA|B D) 821 ! XINT2SAB = XINT2SAB + (C-barB K-barA|L J) 822 CALL CCSDT_TRAN3_R(XINT1SAB,XINT2SAB,XLAMDP0,XLAMDH0, 823 & XLAMDPB,XLAMDHA, XLAMDP0,XLAMDH0, 824 & WORK(KXINT),IDEL) 825 826 ! XINT1SAB = XINT1SAB + (C-barA K-barB|B D) 827 ! XINT2SAB = XINT2SAB + (C-barA K-barB|L J) 828 CALL CCSDT_TRAN3_R(XINT1SAB,XINT2SAB,XLAMDP0,XLAMDH0, 829 & XLAMDPA,XLAMDHB, XLAMDP0,XLAMDH0, 830 & WORK(KXINT),IDEL) 831 832 ! XINT1SAB = XINT1SAB + (C-barB K|B-barA D) 833 ! XINT2SAB = XINT2SAB + (C-barB K|L J-barA) 834 CALL CCSDT_TRAN3_R(XINT1SAB,XINT2SAB,XLAMDP0,XLAMDH0, 835 & XLAMDPB,XLAMDH0, XLAMDPA,XLAMDHA, 836 & WORK(KXINT),IDEL) 837 838 ! XINT1SAB = XINT1SAB + (C-barA K|B-barB D) 839 ! XINT2SAB = XINT2SAB + (C-barA K|L J-barB) 840 CALL CCSDT_TRAN3_R(XINT1SAB,XINT2SAB,XLAMDP0,XLAMDH0, 841 & XLAMDPA,XLAMDH0, XLAMDPB,XLAMDHB, 842 & WORK(KXINT),IDEL) 843 844 ! XINT1SAB = XINT1SAB + (C K-barB|B-barA D) 845 ! XINT2SAB = XINT2SAB + (C K-barB|L J-barA) 846 CALL CCSDT_TRAN3_R(XINT1SAB,XINT2SAB,XLAMDP0,XLAMDH0, 847 & XLAMDP0,XLAMDHB, XLAMDPA,XLAMDHA, 848 & WORK(KXINT),IDEL) 849 850 ! XINT1SAB = XINT1SAB + (C K-barA|B-barB D) 851 ! XINT2SAB = XINT2SAB + (C K-barA|L J-barB) 852 CALL CCSDT_TRAN3_R(XINT1SAB,XINT2SAB,XLAMDP0,XLAMDH0, 853 & XLAMDP0,XLAMDHA, XLAMDPB,XLAMDHB, 854 & WORK(KXINT),IDEL) 855 856 END DO 857 END DO 858 859 RETURN 860 END 861*---------------------------------------------------------------------* 862* END OF SUBROUTINE CCSDT_INTS2_NODDY * 863*---------------------------------------------------------------------* 864