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! 18 SUBROUTINE CCSDT_AAMAT_NODDY(IOPTRES,FREQRES,LABELB,ISYMB, 19 & LISTC,IDLSTC,OEFF_INIT, 20 & OMEGA1,OMEGA2, 21 & OMEGA1EFF,OMEGA2EFF, 22 & IDOTS,DOTPROD,LISTDP,ITRAN, 23 & NXTRAN,MXVEC,WORK,LWORK) 24*---------------------------------------------------------------------* 25* 26* Purpose: compute triples contribution to A{B} transformed vector 27* 28* (A{B} T^C)^eff_1,2 = (A{B} T^C)_1,2(CCSD) + (A{B} T^C)_1,2(T3) 29* - A_1,2;3 (w_3 - w)^1 (A{B} T^C)_3 30* 31* 32* Written by Christof Haettig, April 2002, based on CCSDT_XI_NODDY 33* 34*=====================================================================* 35 IMPLICIT NONE 36#include "priunit.h" 37#include "ccsdinp.h" 38#include "maxorb.h" 39#include "ccsdsym.h" 40#include "ccfield.h" 41#include "ccorb.h" 42 43 LOGICAL LOCDBG, PRINT_T3 44 PARAMETER (LOCDBG=.FALSE., PRINT_T3=.FALSE.) 45 INTEGER ISYM0 46 PARAMETER (ISYM0 = 1) 47 48 CHARACTER LISTDP*3, LABELB*8, LISTC*3 49 LOGICAL OEFF_INIT 50 INTEGER LWORK, IOPTRES, ITRAN, MXVEC, NXTRAN, ISYMB, IDLSTC 51 INTEGER IDOTS(MXVEC,NXTRAN) 52 53#if defined (SYS_CRAY) 54 REAL DOTPROD(MXVEC,NXTRAN), FREQC, FREQRES 55 REAL WORK(LWORK), ONE, TWO, DUMMY, FF, SIGN, DDOT 56 REAL OMEGA1(*),OMEGA2(*) 57 REAL OMEGA1EFF(*),OMEGA2EFF(*) 58#else 59 DOUBLE PRECISION DOTPROD(MXVEC,NXTRAN), FREQC, FREQRES 60 DOUBLE PRECISION WORK(LWORK), ONE, TWO, DUMMY, FF, SIGN, DDOT 61 DOUBLE PRECISION OMEGA1(*),OMEGA2(*) 62 DOUBLE PRECISION OMEGA1EFF(*),OMEGA2EFF(*) 63#endif 64 PARAMETER( ONE = 1.0D0, TWO = 2.0D0 ) 65 66 CHARACTER*10 MODEL 67 INTEGER KFOCKB, KT3AM, KT2AM0, KT2AMC, KEND2, KEND1, KOMEGA1, 68 & KOMEGA2, KOMEGA3, KSCR1, KFOCKD, KFOCK0, KT1AMC, 69 & KFOCKB_AO, KFOCKBC, KLAMPC, KLAMHC, KFOCKC, KINT1SC, 70 & KINT2SC, KINT1S, KINT2S, KXIAJB, KYIAJB, KINT1T, KINT2T, 71 & LWRK1, KDUM, LWRK2, KT1AMP0, KLAMP0, KLAMH0 72 INTEGER IJ, NIJ, LUSIFC, INDEX, IDUMMY, ILSTSYM, ISYMC, LUFOCK, 73 & IRREP, IERR, ILLL, IDEL, ISYDIS, IOPT, ISYMD, IVEC, 74 & IDLSTD, KFCKBUF, KFIELD, KFIELDAO 75 76 INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 77 78*---------------------------------------------------------------------* 79* Begin: 80*---------------------------------------------------------------------* 81 CALL QENTER('CCSDT_AAMAT_NODDY') 82 83 IF (LOCDBG) WRITE(LUPRI,*) 'Entered CCSDT_AAMAT_NODDY...' 84 85 IF(DIRECT)CALL QUIT('DIRECT NOT IMPLEMENTED IN CCSDT_AAMAT_NODDY') 86 87*---------------------------------------------------------------------* 88* Memory allocation: 89*---------------------------------------------------------------------* 90 KEND1 = 1 91 92 KOMEGA1 = KEND1 93 KOMEGA2 = KOMEGA1 + NT1AMX 94 KOMEGA3 = KOMEGA2 + NT1AMX*NT1AMX 95 KEND1 = KOMEGA3 + NT1AMX*NT1AMX*NT1AMX 96 97 KSCR1 = KEND1 98 KFOCKD = KSCR1 + NT1AMX 99 KLAMP0 = KFOCKD + NORBT 100 KLAMH0 = KLAMP0 + NLAMDT 101 KFOCK0 = KLAMH0 + NLAMDT 102 KT1AMP0 = KFOCK0 + NORBT*NORBT 103 KEND1 = KT1AMP0+ NT1AMX 104 105 106 KFOCKB = KEND1 107 KFOCKB_AO = KFOCKB + NORBT*NORBT 108 KFOCKBC = KFOCKB_AO + NORBT*NORBT 109 KFCKBUF = KFOCKBC + NORBT*NORBT 110 KEND1 = KFCKBUF + NORBT*NORBT 111 112 IF (NONHF) THEN 113 KFIELD = KEND1 114 KFIELDAO = KFIELD + NORBT*NORBT 115 KEND1 = KFIELDAO + NORBT*NORBT 116 END IF 117 118 KLAMPC = KEND1 119 KLAMHC = KLAMPC + NLAMDT 120 KFOCKC = KLAMHC + NLAMDT 121 KEND1 = KFOCKC + NORBT*NORBT 122 123 KINT1S = KEND1 124 KINT2S = KINT1S + NT1AMX*NVIRT*NVIRT 125 KEND1 = KINT2S + NRHFT*NRHFT*NT1AMX 126 127 KINT1SC = KEND1 128 KINT2SC = KINT1SC + NT1AMX*NVIRT*NVIRT 129 KEND1 = KINT2SC + NRHFT*NRHFT*NT1AMX 130 131 KXIAJB = KEND1 132 KYIAJB = KXIAJB + NT1AMX*NT1AMX 133 KEND1 = KYIAJB + NT1AMX*NT1AMX 134 135 KINT1T = KEND1 136 KINT2T = KINT1T + NT1AMX*NVIRT*NVIRT 137 KEND1 = KINT2T + NRHFT*NRHFT*NT1AMX 138 139 LWRK1 = LWORK - KEND1 140 IF (LWRK1 .LT. 0) THEN 141 CALL QUIT('Insufficient space in CCSDT_AAMAT_NODDY') 142 ENDIF 143 144 CALL DZERO(WORK(KOMEGA1),NT1AMX) 145 146*---------------------------------------------------------------------* 147* Get zeroth-order Lambda matrices: 148*---------------------------------------------------------------------* 149 IOPT = 1 150 KDUM = KEND1 151 Call CC_RDRSP('R0',0,ISYM0,IOPT,MODEL,WORK(KT1AMP0),WORK(KDUM)) 152 153 Call LAMMAT(WORK(KLAMP0),WORK(KLAMH0),WORK(KT1AMP0), 154 & WORK(KEND1),LWRK1) 155 156*---------------------------------------------------------------------* 157* Read SCF orbital energies from file: 158*---------------------------------------------------------------------* 159 LUSIFC = -1 160 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY, 161 & .FALSE.) 162 REWIND LUSIFC 163 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 164 READ (LUSIFC) 165 READ (LUSIFC) (WORK(KFOCKD+I-1), I=1,NORBT) 166 CALL GPCLOSE(LUSIFC,'KEEP') 167 168*---------------------------------------------------------------------* 169* read zeroth-order AO Fock matrix from file: 170*---------------------------------------------------------------------* 171 LUFOCK = -1 172 CALL GPOPEN(LUFOCK,'CC_FCKH','OLD',' ','UNFORMATTED',IDUMMY, 173 & .FALSE.) 174 REWIND(LUFOCK) 175 READ(LUFOCK) (WORK(KFOCK0-1+I),I=1,N2BST(ISYM0)) 176 CALL GPCLOSE(LUFOCK,'KEEP') 177 178 CALL CC_FCKMO(WORK(KFOCK0),WORK(KLAMP0),WORK(KLAMH0), 179 & WORK(KEND1),LWRK1,ISYM0,ISYM0,ISYM0) 180 181*---------------------------------------------------------------------* 182* If needed get external field: 183*---------------------------------------------------------------------* 184 IF ((NONHF) .AND. NFIELD .GT. 0) THEN 185 CALL DZERO(WORK(KFIELDAO),NORBT*NORBT) 186 DO I = 1, NFIELD 187 FF = EFIELD(I) 188 CALL CC_ONEP(WORK(KFIELDAO),WORK(KEND1),LWRK1,FF,1,LFIELD(I)) 189 ENDDO 190 191 CALL DCOPY(NORBT*NORBT,WORK(KFIELDAO),1,WORK(KFIELD),1) 192 CALL CC_FCKMO(WORK(KFIELD),WORK(KLAMP0),WORK(KLAMH0), 193 * WORK(KEND1),LWRK1,1,1,1) 194 ENDIF 195 196*---------------------------------------------------------------------* 197* Matrix with property integrals in MO basis: 198*---------------------------------------------------------------------* 199 ! read property integrals from file: 200 CALL CCPRPAO(LABELB,.TRUE.,WORK(KFOCKB_AO),IRREP,ISYMB,IERR, 201 & WORK(KEND1),LWRK1) 202 CALL DCOPY(NORBT*NORBT,WORK(KFOCKB_AO),1,WORK(KFOCKB),1) 203 IF ((IERR.GT.0) .OR. (IERR.EQ.0 .AND. IRREP.NE.ISYMB)) THEN 204 CALL QUIT('CCSDT_AAMAT_NODDY: error reading operator '//LABELB) 205 ELSE IF (IERR.LT.0) THEN 206 CALL DZERO(WORK(KFOCKB),N2BST(ISYMB)) 207 END IF 208 209 ! transform property integrals to Lambda-MO basis 210 CALL CC_FCKMO(WORK(KFOCKB),WORK(KLAMP0),WORK(KLAMH0), 211 & WORK(KEND1),LWRK1,ISYMB,1,1) 212 213*---------------------------------------------------------------------* 214* Compute some integrals: 215* XINT1S = (CK|BD) 216* XINT2S = (CK|LJ) 217* XINT1T = (KC|BD) 218* XINT2T = (KC|LJ) 219* XIAJB = 2(IA|JB) - (IB|JA) 220* YIAJB = (IA|JB) 221*---------------------------------------------------------------------* 222 CALL CCSDT_INTS0_NODDY(.TRUE.,WORK(KXIAJB),WORK(KYIAJB), 223 & .TRUE.,WORK(KINT1S),WORK(KINT2S), 224 & .TRUE.,WORK(KINT1T),WORK(KINT2T), 225 & WORK(KLAMP0),WORK(KLAMH0), 226 & WORK(KEND1),LWRK1) 227 228*---------------------------------------------------------------------* 229* allocate work space for a triples and two doubles vectors and 230* read zeroth-order and response doubles amplitudes from file, 231* compute response Lambda matrices: 232*---------------------------------------------------------------------* 233 234 KT3AM = KEND1 235 KT2AM0 = KT3AM + NT1AMX*NT1AMX*NT1AMX 236 KT2AMC = KT2AM0 + NT1AMX*NT1AMX 237 KT1AMC = KT2AMC + NT1AMX*NT1AMX 238 KEND2 = KT1AMC + NT1AMX 239 LWRK2 = LWORK - KEND2 240 IF (LWRK2 .LT. 0) THEN 241 CALL QUIT('Insufficient space in CCSDT_AAMAT_NODDY') 242 ENDIF 243 244 IOPT = 2 245 CALL CC_RDRSP('R0',0,ISYMOP,IOPT,MODEL,DUMMY,WORK(KT3AM)) 246 CALL CC_T2SQ(WORK(KT3AM),WORK(KT2AM0),ISYMOP) 247 248 ISYMC = ILSTSYM(LISTC,IDLSTC) 249 IOPT = 3 250 CALL CC_RDRSP(LISTC,IDLSTC,ISYMC,IOPT,MODEL, 251 & WORK(KT1AMC),WORK(KT3AM)) 252 CALL CCLR_DIASCL(WORK(KT3AM),TWO,ISYMC) 253 CALL CC_T2SQ(WORK(KT3AM),WORK(KT2AMC),ISYMC) 254 255 CALL CCLR_LAMTRA(WORK(KLAMP0),WORK(KLAMPC), 256 & WORK(KLAMH0),WORK(KLAMHC),WORK(KT1AMC),ISYMC) 257 258*---------------------------------------------------------------------* 259* compute zeroth-order triples response amplitudes and its 260* contributions to the result vector 261*---------------------------------------------------------------------* 262 263C -------------------------------------------------------------- 264C compute zeroth-order triples; note: CCSDT_T3AM uses in non-hf 265C case WORK(KOMEGA3) as scratch vector for the iterative 266C solution of the triples equations! 267C -------------------------------------------------------------- 268 CALL DZERO(WORK(KT3AM),NT1AMX*NT1AMX*NT1AMX) 269 CALL CCSDT_T03AM(WORK(KT3AM),WORK(KINT1S),WORK(KINT2S), 270 & WORK(KT2AM0),WORK(KSCR1),WORK(KFOCKD), 271 & WORK(KFIELD),WORK(KOMEGA3)) 272 273C -------------------------------------------------------------- 274C calculate one-index transf. property integrals: B^C = [B,T1^C] 275C -------------------------------------------------------------- 276 CALL DCOPY(NORBT*NORBT,WORK(KFOCKB_AO),1,WORK(KFCKBUF),1) 277 CALL DCOPY(NORBT*NORBT,WORK(KFOCKB_AO),1,WORK(KFOCKBC),1) 278 CALL CC_FCKMO(WORK(KFCKBUF),WORK(KLAMPC),WORK(KLAMH0), 279 & WORK(KEND2),LWRK2,ISYMB,ISYMC,ISYM0) 280 CALL CC_FCKMO(WORK(KFOCKBC),WORK(KLAMP0),WORK(KLAMHC), 281 & WORK(KEND2),LWRK2,ISYMB,ISYM0,ISYMC) 282 CALL DAXPY(NORBT*NORBT,ONE,WORK(KFCKBUF),1,WORK(KFOCKBC),1) 283 284C --------------------------------- 285C Initialize triples result vector: 286C --------------------------------- 287 CALL DZERO(WORK(KOMEGA3),NT1AMX*NT1AMX*NT1AMX) 288 289C ---------------------------------------------- 290C add triples contribution: <mu_3|[B^C,T3^0]|HF> 291C ---------------------------------------------- 292 CALL CCSDT_XKSI3_2(WORK(KOMEGA3),WORK(KFOCKBC),WORK(KT3AM)) 293 294*---------------------------------------------------------------------* 295* compute contribution from doubles amplitudes: 296*---------------------------------------------------------------------* 297 298C --------------------------------------------------- 299C add triples contribution: <mu_3|[[B,T2^C],T2^0]|HF> 300C --------------------------------------------------- 301 CALL CCSDT_XKSI3_1(WORK(KOMEGA3),WORK(KFOCKB), 302 & WORK(KT2AM0),WORK(KT2AMC),ONE) 303 CALL CCSDT_XKSI3_1(WORK(KOMEGA3),WORK(KFOCKB), 304 & WORK(KT2AMC),WORK(KT2AM0),ONE) 305 306*---------------------------------------------------------------------* 307* compute triples response amplitudes and its contributions 308* to the result vector 309*---------------------------------------------------------------------* 310 IF (LISTC(1:3).EQ.'R1 ' .OR. LISTC(1:3).EQ.'RE ' .OR. 311 & LISTC(1:3).EQ.'RC ' ) THEN 312 313C ----------------------------------------- 314C calculate first-order triples amplitudes: 315C ----------------------------------------- 316 KDUM = KEND2 317 CALL CCSDT_T31_NODDY(WORK(KT3AM),LISTC,IDLSTC,FREQC,.FALSE., 318 & .FALSE.,WORK(KINT1S),WORK(KINT2S), 319 & .FALSE.,WORK(KDUM),WORK(KDUM), 320 & .FALSE.,WORK(KDUM),WORK(KDUM), 321 & WORK(KINT1SC),WORK(KINT2SC), 322 & WORK(KLAMPC),WORK(KLAMHC),WORK(KFOCKC), 323 & WORK(KLAMP0),WORK(KLAMH0),WORK(KFOCK0), 324 & WORK(KDUM),WORK(KFOCKD), 325 & WORK(KEND2),LWRK2) 326 327 ELSE IF (LISTC(1:3).EQ.'R2 ' .OR. LISTC(1:3).EQ.'ER1') THEN 328 329C ------------------------------------------ 330C calculate second-order triples amplitudes: 331C ------------------------------------------ 332 CALL CCSDT_T32_NODDY(WORK(KT3AM),LISTC,IDLSTC,FREQC, 333 & WORK(KINT1S),WORK(KINT2S), 334 & WORK(KLAMP0),WORK(KLAMH0),WORK(KFOCK0), 335 & WORK(KFOCKD),WORK(KFIELDAO),WORK(KFIELD), 336 & WORK(KSCR1),WORK(KEND2),LWRK2) 337 338 ELSE 339 340 CALL QUIT('Unknown list '//LISTC//' in CCSDT_AA_NODDY.') 341 342 END IF 343 344 CALL DSCAL(NT1AMX*NT1AMX*NT1AMX,-1.0D0,WORK(KT3AM),1) 345 346 IF (LOCDBG) THEN 347 WRITE(LUPRI,*) 'CCSDT_AAMAT_NODDY> T^C:' 348 WRITE(LUPRI,*) 'CCSDT_AAMAT_NODDY> LISTC,IDLSTC:',LISTC,IDLSTC 349 CALL PRINT_PT3_NODDY(WORK(KT3AM)) 350 END IF 351 352C ----------------------------------------------- 353C add contribution to doubles: <mu_2|[B,T3^C]|HF> 354C ----------------------------------------------- 355 CALL DZERO(WORK(KOMEGA2),NT1AMX*NT1AMX) 356 CALL CCSDT_XKSI2_2(WORK(KOMEGA2),WORK(KFOCKB),WORK(KT3AM)) 357 358 DO I = 1,NT1AMX 359 DO J = 1,I 360 IJ = NT1AMX*(I-1) + J 361 NIJ = INDEX(I,J) 362 OMEGA2(NIJ) = OMEGA2(NIJ) + WORK(KOMEGA2+IJ-1) 363 IF (.NOT. OEFF_INIT) 364 & OMEGA2EFF(NIJ) = OMEGA2EFF(NIJ) + WORK(KOMEGA2+IJ-1) 365 END DO 366 END DO 367 368C ------------------------------------------- 369C contribution to triples: <mu_3|[B,T3^C]|HF> 370C ------------------------------------------- 371 CALL CCSDT_XKSI3_2(WORK(KOMEGA3),WORK(KFOCKB),WORK(KT3AM)) 372 373 if (print_t3) then 374 write(lupri,*)'CCSDT_AA_NODDY> vector type:',listc 375 write(lupri,*)'CCSDT_AA_NODDY> idlst:',idlstc 376 write(lupri,*)'CCSDT_AA_NODDY> freq:',freqres 377 call ccsdt_clean_t3(work(komega3),nt1amx,nvirt,nrhft) 378 call print_pt3_noddy(work(komega3)) 379 end if 380 381*---------------------------------------------------------------------* 382* Now we split: 383* for IOPTRES < 5 we compute the effective result vector 384* for IOPTRES = 5 we compute the contractions Tbar^D A{B} T^C 385*---------------------------------------------------------------------* 386 IF (IOPTRES.GE.1 .AND. IOPTRES.LE.4) THEN 387 388 IF (OEFF_INIT) THEN 389 CALL DCOPY(NT1AMX,OMEGA1,1,OMEGA1EFF,1) 390 CALL DCOPY(NT2AMX,OMEGA2,1,OMEGA2EFF,1) 391 END IF 392 393 CALL CC_RHPART_NODDY(OMEGA1EFF,OMEGA2EFF,WORK(KOMEGA3),FREQRES, 394 & WORK(KFOCKD),WORK(KFOCK0),WORK(KFIELD), 395 & WORK(KXIAJB),WORK(KINT1T),WORK(KINT2T), 396 & WORK(KEND1),LWRK1) 397 398 ELSE IF (IOPTRES.EQ.5) THEN 399 400 SIGN = -1.0D0 401 CALL CCDOTRSP_NODDY(WORK(KOMEGA1),WORK(KOMEGA2), 402 & WORK(KOMEGA3),SIGN, 403 & ITRAN,LISTDP,IDOTS,DOTPROD,MXVEC, 404 & WORK(KLAMP0),WORK(KLAMH0), 405 & WORK(KFOCK0),WORK(KFOCKD), 406 & WORK(KXIAJB),WORK(KYIAJB), 407 & WORK(KINT1T),WORK(KINT2T), 408 & WORK(KINT1S),WORK(KINT2S), 409 & 'CCSDT_AAMAT_NODDY',LOCDBG,LOCDBG,.FALSE., 410 & WORK(KEND1),LWRK1) 411 412 ELSE 413 CALL QUIT('Illegal value for IOPTRES IN CCSDT_AAMAT_NODDY') 414 END IF 415 416 CALL QEXIT('CCSDT_AAMAT_NODDY') 417 RETURN 418 END 419 420*---------------------------------------------------------------------* 421* END OF SUBROUTINE CCSDT_AAMAT_NODDY * 422*---------------------------------------------------------------------* 423