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_ETA_NODDY(LISTL,IDLSTL,IOPTRES, 21 & FOCKA,FOCKA_AO,FREQE,FOCK0, 22 & OMEGA1,OMEGA2, 23 & OMEGA1EFF,OMEGA2EFF, 24 & IDOTS,DOTPROD,LISTDP,ITRAN, 25 & NEATRAN,MXVEC,WORK,LWORK ) 26*---------------------------------------------------------------------* 27* 28* Purpose: compute triples contribution to Eta vector 29* 30* Eta^eff_1,2 = Eta_1,2(CCSD) + Eta_1,2(T^0_3) 31* - Eta_3 A_3;1,2 (w_3 - w)^1 32* 33* 34* Written by Christof Haettig, April 2002 35* based on different other noddy codes 36* 37*=====================================================================* 38 IMPLICIT NONE 39#include "dummy.h" 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 "ccnoddy.h" 47 48 LOGICAL LOCDBG, FD_TEST, XI_ONLY 49 PARAMETER (LOCDBG=.FALSE., FD_TEST=.false., 50 & XI_ONLY = .FALSE.) 51 52 INTEGER ISYM0 53 PARAMETER (ISYM0 = 1) 54 55 CHARACTER*3 LISTL, LISTDP 56 INTEGER LWORK, IDLSTL, IOPTRES, ITRAN, MXVEC, NEATRAN 57 INTEGER IDOTS(MXVEC,NEATRAN) 58 59#if defined (SYS_CRAY) 60 REAL DOTPROD(MXVEC,NEATRAN), DDOT 61 REAL WORK(LWORK), FREQL, FREQE, FREQC, SIGN 62 REAL FOCKA_AO(NORBT,NORBT) 63 REAL OMEGA1(*), OMEGA2(*), FOCKA(NORBT,NORBT) 64 REAL OMEGA1EFF(*), OMEGA2EFF(*), FOCK0(NORBT,NORBT) 65 REAL SIXTH, TWO, TCON, DCON, SCON, FF 66#else 67 DOUBLE PRECISION DOTPROD(MXVEC,NEATRAN), DDOT 68 DOUBLE PRECISION FOCKA_AO(NORBT,NORBT) 69 DOUBLE PRECISION WORK(LWORK), FREQL, FREQE, FREQC, SIGN 70 DOUBLE PRECISION OMEGA1(*), OMEGA2(*), FOCKA(NORBT,NORBT) 71 DOUBLE PRECISION OMEGA1EFF(*), OMEGA2EFF(*), FOCK0(NORBT,NORBT) 72 DOUBLE PRECISION SIXTH, TWO, TCON, DCON, SCON, FF 73#endif 74 PARAMETER(SIXTH=1.0D0/6.0D0, TWO=2.0D0) 75 76 CHARACTER*10 MODEL 77 LOGICAL L2INCL 78 INTEGER INDEX, LUSIFC, IOPT, ISYMD, ILLL, IDEL, ISYDIS, NIJ, IJ, 79 & IVEC, IDLSTC, ISYMC, LUFOCK, ILSTSYM, ISYML, LUTEMP, IDX 80 INTEGER KSCR1, KFOCKD, KEND1, KT1AMP0, KLAMP0, KLAMH0, KFIELD, 81 & KINT1T0, KINT2T0, KINT1S0, KINT2S0, KXIAJB, KYIAJB, 82 & K0IOVVO, K0IOOVV, K0IOOOO, K0IVVVV, KOME1, KOME2, KDUM, 83 & KXINT, KEND2, LWRK2, KL1AM, KL2AM, KL3AM, KT3AM, KT2AM, 84 & K3AM, KEND3, LWRK3, KINT1SC, KINT2SC, KLAMPC, KLAMHC, 85 & KFOCKC, LWRK1, KE3AM, KTC3AM, KTC1AM, KTC2AM, 86 & KFDETA1, KFDETA2, KFDETA3, KEND4, LWRK4, KMMAT, KDIAM, 87 & KT3SCR, KFIELDAO, KFOCK0 88 89 INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 90 91 CALL QENTER('CCSDT_ETA_NODDY') 92 KDUM = 1 93 94 IF (DIRECT) CALL QUIT('CCSDT_ETA_NODDY: DIRECT NOT IMPLEMENTED') 95 96 IF (LOCDBG) THEN 97 WRITE(LUPRI,*) 'CCSDT_ETA_NODDY> Eta vector on entry:' 98 CALL CC_PRP(OMEGA1,OMEGA2,1,1,1) 99 END IF 100 101*---------------------------------------------------------------------* 102* Memory allocation: 103*---------------------------------------------------------------------* 104 KSCR1 = 1 105 KFOCKD = KSCR1 + NT1AMX 106 KFOCK0 = KFOCKD + NORBT 107 KEND1 = KFOCK0 + NORBT*NORBT 108 109 IF (NONHF) THEN 110 KFIELD = KEND1 111 KFIELDAO = KFIELD + NBAST*NBAST 112 KEND1 = KFIELDAO + NBAST*NBAST 113 END IF 114 115 KT1AMP0 = KEND1 116 KLAMP0 = KT1AMP0 + NT1AMX 117 KLAMH0 = KLAMP0 + NLAMDT 118 KEND1 = KLAMH0 + NLAMDT 119 120 KINT1T0 = KEND1 121 KINT2T0 = KINT1T0 + NT1AMX*NVIRT*NVIRT 122 KEND1 = KINT2T0 + NRHFT*NRHFT*NT1AMX 123 124 KINT1S0 = KEND1 125 KINT2S0 = KINT1S0 + NT1AMX*NVIRT*NVIRT 126 KEND1 = KINT2S0 + NRHFT*NRHFT*NT1AMX 127 128 KXIAJB = KEND1 129 KYIAJB = KXIAJB + NT1AMX*NT1AMX 130 KEND1 = KYIAJB + NT1AMX*NT1AMX 131 132 K0IOVVO = KEND1 133 K0IOOVV = K0IOVVO + NRHFT*NVIRT*NVIRT*NRHFT 134 K0IOOOO = K0IOOVV + NRHFT*NVIRT*NVIRT*NRHFT 135 K0IVVVV = K0IOOOO + NRHFT*NRHFT*NRHFT*NRHFT 136 KEND1 = K0IVVVV + NVIRT*NVIRT*NVIRT*NVIRT 137 138 KOME1 = KEND1 139 KOME2 = KOME1 + NT1AMX 140 KEND1 = KOME2 + NT1AMX*NT1AMX 141 142 LWRK1 = LWORK - KEND1 143 IF (LWRK1 .LT. 0) THEN 144 CALL QUIT('Insufficient space in CCSDT_ETA_NODDY') 145 ENDIF 146 147*---------------------------------------------------------------------* 148* Get zeroth-order Lambda matrices: 149*---------------------------------------------------------------------* 150 IOPT = 1 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 precalculated integrals from file: 158*---------------------------------------------------------------------* 159 CALL CCSDT_READ_NODDY(.TRUE.,WORK(KFOCKD),WORK(KFOCK0), 160 & WORK(KFIELD),WORK(KFIELDAO), 161 & .TRUE.,WORK(KXIAJB),WORK(KYIAJB), 162 & .TRUE.,WORK(KINT1S0),WORK(KINT2S0), 163 & .TRUE.,WORK(KINT1T0),WORK(KINT2T0), 164 & .TRUE.,WORK(K0IOVVO),WORK(K0IOOVV), 165 & WORK(K0IOOOO),WORK(K0IVVVV), 166 & NORBT,NLAMDT,NRHFT,NVIRT,NT1AMX) 167 168*---------------------------------------------------------------------* 169* Some more memory allocations: 170*---------------------------------------------------------------------* 171 KL1AM = KEND1 172 KL2AM = KL1AM + NT1AMX 173 KL3AM = KL2AM + NT1AMX*NT1AMX 174 KEND2 = KL3AM + NT1AMX*NT1AMX*NT1AMX 175 LWRK2 = LWORK - KEND2 176 177 KT3AM = KEND2 178 KT2AM = KT3AM + NT1AMX*NT1AMX*NT1AMX 179 KEND3 = KT2AM + NT1AMX*NT1AMX 180 LWRK3 = LWORK - KEND3 181 182 KEND4 = KEND3 183 IF (NONHF) THEN 184 ! allocate scratch array needed in CCSDT_T03AM and 185 ! also in CCSDT_L03AM in case of finite difference runs 186 KT3SCR = KEND4 187 KEND4 = KT3SCR + NT1AMX*NT1AMX*NT1AMX 188 END IF 189 LWRK4 = LWORK - KEND4 190 191 IF (LWRK4 .LT. 0) THEN 192 CALL QUIT('Insufficient space in CCSDT_ETA_NODDY') 193 ENDIF 194 195 ! read T^0 doubles amplitudes from file and square up 196 IOPT = 2 197 Call CC_RDRSP('R0',0,ISYM0,IOPT,MODEL,WORK(KDUM),WORK(KT3AM)) 198 CALL CC_T2SQ(WORK(KT3AM),WORK(KT2AM),ISYM0) 199 200 ! read L^0 multipliers from file and square up doubles part 201 ISYML = ILSTSYM(LISTL,IDLSTL) 202 IOPT = 3 203 Call CC_RDRSP(LISTL,IDLSTL,ISYML,IOPT,MODEL, 204 & WORK(KL1AM),WORK(KT3AM)) 205 CALL CC_T2SQ(WORK(KT3AM),WORK(KL2AM),ISYM0) 206 207*---------------------------------------------------------------------* 208* Compute zeroth-order triples amplitudes: 209*---------------------------------------------------------------------* 210 LUTEMP = -1 211 CALL GPOPEN(LUTEMP,FILNODT30,'UNKNOWN',' ','UNFORMATTED', 212 & IDUMMY,.FALSE.) 213 READ(LUTEMP) (WORK(KT3AM+I-1), I=1,NT1AMX*NT1AMX*NT1AMX) 214 CALL GPCLOSE(LUTEMP,'KEEP') 215 216 IF (NONHF) THEN 217 LUTEMP = -1 218 CALL GPOPEN(LUTEMP,'T3AMPL','UNKNOWN',' ','UNFORMATTED', 219 & IDUMMY,.FALSE.) 220 REWIND LUTEMP 221 WRITE (LUTEMP) (WORK(KT3AM-1+IDX),IDX=1,NT1AMX*NT1AMX*NT1AMX) 222 CALL GPCLOSE(LUTEMP,'KEEP') 223 END IF 224*---------------------------------------------------------------------* 225* Compute L^0_3 multipliers: 226*---------------------------------------------------------------------* 227 IF (LISTL(1:3).EQ.'L0 ') THEN 228 229 ! remember that CCSDT_L03AM returns -L3 !! 230 CALL CCSDT_L03AM(WORK(KL3AM),WORK(KINT1T0),WORK(KINT2T0), 231 * WORK(KXIAJB),FOCK0,WORK(KL1AM), 232 * WORK(KL2AM),WORK(KSCR1),WORK(KFOCKD), 233 * WORK(KFIELD),WORK(KT3SCR)) 234 235 CALL DSCAL(NT1AMX*NT1AMX*NT1AMX,-1.0D0,WORK(KL3AM),1) 236 237 ELSE IF (LISTL(1:3).EQ.'L1 ' .OR. LISTL(1:3).EQ.'LE ' .OR. 238 & LISTL(1:3).EQ.'M1 ' .OR. LISTL(1:3).EQ.'N2 ' .OR. 239 & LISTL(1:3).EQ.'E0 ' ) THEN 240 241 CALL CCSDT_TBAR31_NODDY(WORK(KL3AM),FREQL,LISTL,IDLSTL, 242 & WORK(KLAMP0),WORK(KLAMH0), 243 & FOCK0,WORK(KFOCKD),WORK(KSCR1), 244 & WORK(KXIAJB),WORK(KINT1T0),WORK(KINT2T0), 245 & WORK(KEND3),LWRK3) 246 247 CALL DSCAL(NT1AMX*NT1AMX*NT1AMX,-1.0D0,WORK(KL3AM),1) 248 249 ELSE 250 251 CALL QUIT('CCSDT_ETA_NODDY> LISTL NOT AVAILABLE:'//LISTL) 252 253 END IF 254 255*---------------------------------------------------------------------* 256* Compute contribution from <L_3|[[A,T^0_3],\tau_nu_1|HF>: 257*---------------------------------------------------------------------* 258 CALL DZERO(WORK(KOME1),NT1AMX) 259 260 CALL CCSDT_E1AM(WORK(KOME1),WORK(KL3AM),WORK(KT3AM),FOCKA) 261 262 DO I = 1,NT1AMX 263 OMEGA1(I) = OMEGA1(I) + WORK(KOME1+I-1) 264 END DO 265 266 IF (LOCDBG) THEN 267 WRITE(LUPRI,*) 'CCSDT_ETA_NODDY> Contribution to Eta1:' 268 CALL CC_PRP(WORK(KOME1),WORK,1,1,0) 269 END IF 270 271*---------------------------------------------------------------------* 272* Compute contribution from <L_3|[[A,T^0_2],\tau_nu_2]|HF> 273*---------------------------------------------------------------------* 274 CALL DZERO(WORK(KOME2),NT1AMX*NT1AMX) 275 276 CALL CCSDT_E2AM(WORK(KOME2),WORK(KL3AM),WORK(KT2AM),FOCKA) 277 278 DO I = 1,NT1AMX 279 DO J = 1,I 280 IJ = NT1AMX*(I-1) + J 281 NIJ = INDEX(I,J) 282 OMEGA2(NIJ) = OMEGA2(NIJ) + WORK(KOME2+IJ-1) 283 END DO 284 END DO 285 286*---------------------------------------------------------------------* 287* finite difference test: 288*---------------------------------------------------------------------* 289 IF (FD_TEST) THEN 290 KFDETA1 = KEND3 291 KFDETA2 = KFDETA1 + NT1AMX 292 KFDETA3 = KFDETA2 + NT1AMX*NT1AMX 293 KEND4 = KFDETA3 + NT1AMX*NT1AMX*NT1AMX 294 LWRK4 = LWORK - KEND4 295 IF (LWRK4 .LT. 0) THEN 296 CALL QUIT('Insufficient space in CCSDT_ETA_NODDY') 297 ENDIF 298 299 CALL CCSDT_ETA_FD(WORK(KT1AMP0),WORK(KT2AM),WORK(KT3AM), 300 & WORK(KL3AM),WORK(KL2AM), 301 & WORK(KFDETA1),WORK(KFDETA2),WORK(KFDETA3), 302 & FOCKA,FOCKA_AO, 303 & WORK(KSCR1),WORK(KEND4),LWRK4) 304 305 WRITE(LUPRI,*) 'CCSDT_ETA_NODDY> my ETA1:' 306 CALL OUTPUT(WORK(KOME1),1,NVIRT,1,NRHFT, 307 & NVIRT,NRHFT,1,LUPRI) 308 WRITE(LUPRI,*) 'CCSDT_ETA_NODDY> finite difference Eta1:' 309 CALL OUTPUT(WORK(KFDETA1),1,NVIRT,1,NRHFT, 310 & NVIRT,NRHFT,1,LUPRI) 311 CALL DAXPY(NT1AMX,-1.0D0,WORK(KOME1),1,WORK(KFDETA1),1) 312 WRITE(LUPRI,*) 'CCSDT_ETA_NODDY> norm of difference:', 313 & DDOT(NT1AMX,WORK(KFDETA1),1,WORK(KFDETA1),1) 314 315 WRITE(LUPRI,*) 'CCSDT_ETA_NODDY> my ETA2:' 316 CALL OUTPUT(WORK(KOME2),1,NT1AMX,1,NT1AMX, 317 & NT1AMX,NT1AMX,1,LUPRI) 318 WRITE(LUPRI,*) 'CCSDT_ETA_NODDY> finite difference Eta2:' 319 CALL OUTPUT(WORK(KFDETA2),1,NT1AMX,1,NT1AMX, 320 & NT1AMX,NT1AMX,1,LUPRI) 321 CALL DAXPY(NT1AMX*NT1AMX,-1.0D0,WORK(KOME2),1,WORK(KFDETA2),1) 322 WRITE(LUPRI,*) 'CCSDT_ETA_NODDY> norm of difference:', 323 & DDOT(NT1AMX*NT1AMX,WORK(KFDETA2),1,WORK(KFDETA2),1) 324 325 WRITE(LUPRI,*) 'CCSDT_ETA_NODDY> finite difference Eta3:' 326 CALL OUTPUT(WORK(KFDETA3),1,NT1AMX*NT1AMX,1,NT1AMX, 327 & NT1AMX*NT1AMX,NT1AMX,1,LUPRI) 328 329 END IF 330 331*---------------------------------------------------------------------* 332* Compute triples result vector 333* <L_2|[A,\tau_nu_3]|HF> + <L_3|[A,\tau_nu_3]|HF> , 334* (CCSDT_E3AM accounts for the wrong sign of L_3) 335*---------------------------------------------------------------------* 336 ! overwrite T3 vector 337 KE3AM = KT3AM 338 339 CALL DZERO(WORK(KE3AM),NT1AMX*NT1AMX*NT1AMX) 340 341 L2INCL = .TRUE. 342 CALL CCSDT_E3AM(WORK(KE3AM),WORK(KL2AM),WORK(KL3AM),FOCKA,L2INCL) 343 344 IF (FD_TEST) THEN 345 WRITE(LUPRI,*) 'CCSDT_ETA_NODDY> my Eta3:' 346 CALL OUTPUT(WORK(KE3AM),1,NT1AMX*NT1AMX,1,NT1AMX, 347 & NT1AMX*NT1AMX,NT1AMX,1,LUPRI) 348 CALL DAXPY(NT1AMX*NT1AMX*NT1AMX,-1.0D0,WORK(KE3AM),1, 349 & WORK(KFDETA3),1) 350 WRITE(LUPRI,*) 'CCSDT_ETA_NODDY> norm of difference:', 351 & DDOT(NT1AMX*NT1AMX*NT1AMX,WORK(KFDETA3),1,WORK(KFDETA3),1) 352 END IF 353*---------------------------------------------------------------------* 354* Now we split: 355* for IOPTRES < 5 we compute the effective Eta vector 356* for IOPTRES = 5 we compute the contractions Eta^A T^B 357*---------------------------------------------------------------------* 358 IF (IOPTRES.GE.1 .AND. IOPTRES.LE.4) THEN 359 360 CALL DCOPY(NT1AMX,OMEGA1,1,OMEGA1EFF,1) 361 CALL DCOPY(NT2AMX,OMEGA2,1,OMEGA2EFF,1) 362 363 CALL CC_LHPART_NODDY(OMEGA1EFF,OMEGA2EFF,WORK(KE3AM),-FREQE, 364 & WORK(KFOCKD),WORK(KFIELD), 365 & WORK(K0IOOOO),WORK(K0IOVVO), 366 & WORK(K0IOOVV),WORK(K0IVVVV), 367 & WORK(KT2AM),WORK(KINT1S0),WORK(KINT2S0), 368 & WORK(KEND3),LWRK3) 369 370 371 ELSE IF (IOPTRES.EQ.5) THEN 372 373 SIGN = +1.0D0 374 CALL CCDOTRSP_NODDY(WORK(KOME1),WORK(KOME2),WORK(KE3AM),SIGN, 375 & ITRAN,LISTDP,IDOTS,DOTPROD,MXVEC, 376 & WORK(KLAMP0),WORK(KLAMH0), 377 & FOCK0,WORK(KFOCKD), 378 & WORK(KXIAJB), WORK(KYIAJB), 379 & WORK(KINT1T0),WORK(KINT2T0), 380 & WORK(KINT1S0),WORK(KINT2S0), 381 & 'CCSDT_ETA_NODDY',LOCDBG,LOCDBG,.FALSE., 382 & WORK(KEND3),LWRK3) 383 384 ELSE 385 CALL QUIT('Illegal value for IOPTRES IN CCSDT_ETA_NODDY') 386 END IF 387 388 CALL QEXIT('CCSDT_ETA_NODDY') 389 390 RETURN 391 END 392*---------------------------------------------------------------------* 393* END OF SUBROUTINE CCSDT_ETA_NODDY * 394*---------------------------------------------------------------------* 395*=====================================================================* 396 SUBROUTINE CCSDT_E3AM(E3AM,L2AM,L3AM,FOCKA,L2INCL) 397*---------------------------------------------------------------------* 398* 399* Purpose: compute triples exictation amplitudes of eta vector 400* 401* Eta_nu_3 = <L_2|[A,tau_nu_3]|HF> 402* + <L_3|[A,tau_nu_3]|HF> 403* 404* if L2INCL=.false. the first contribution (from L2) is skipped 405* 406* Written by Christof Haettig, April 2002 407*=====================================================================* 408#if defined (IMPLICIT_NONE) 409 IMPLICIT NONE 410#else 411# include "implicit.h" 412#endif 413#include "priunit.h" 414#include "ccsdinp.h" 415#include "maxorb.h" 416#include "ccsdsym.h" 417#include "ccfield.h" 418#include "ccorb.h" 419 420 LOGICAL L2INCL 421 422#if defined (SYS_CRAY) 423 REAL AIBJCK 424 REAL E3AM(NT1AMX,NT1AMX,NT1AMX) 425 REAL L2AM(NT1AMX,NT1AMX) 426 REAL L3AM(NT1AMX,NT1AMX,NT1AMX) 427 REAL FOCKA(NORBT,NORBT) 428 REAL HALF 429#else 430 DOUBLE PRECISION AIBJCK 431 DOUBLE PRECISION E3AM(NT1AMX,NT1AMX,NT1AMX) 432 DOUBLE PRECISION L2AM(NT1AMX,NT1AMX) 433 DOUBLE PRECISION L3AM(NT1AMX,NT1AMX,NT1AMX) 434 DOUBLE PRECISION FOCKA(NORBT,NORBT) 435 DOUBLE PRECISION HALF, THIRD 436#endif 437 PARAMETER ( HALF = 0.5D0, THIRD = 1.0D0/3.0D0 ) 438 439 INTEGER NK, NC, NCK, NJ, NB, NBJ, NBK, NI, NA, NAI, NL, ND, 440 & NDJ, NBL, NCI, NAK, NBI, NAJ 441 442 DO NK = 1,NRHFT 443 DO NC = 1,NVIRT 444 NCK = NVIRT*(NK-1) + NC 445 446 DO NJ = 1,NRHFT 447 DO NB = 1,NVIRT 448 NBJ = NVIRT*(NJ-1) + NB 449 NBK = NVIRT*(NK-1) + NB 450 451 DO NI = 1,NRHFT 452 DO NA = 1,NVIRT 453 NAI = NVIRT*(NI-1) + NA 454 455 AIBJCK = 0.0D0 456 457 IF (L2INCL) THEN 458 AIBJCK = AIBJCK + 459 & ( + FOCKA(NK,NRHFT+NC) * L2AM(NAI,NBJ) 460 & - FOCKA(NJ,NRHFT+NC) * L2AM(NAI,NBK) ) 461 END IF 462 463 DO ND = 1, NVIRT 464 NDJ = NVIRT*(NJ-1) + ND 465 AIBJCK = AIBJCK 466 & + HALF*L3AM(NAI,NDJ,NCK) * FOCKA(NRHFT+ND,NRHFT+NB) 467 END DO 468 469 DO NL = 1, NRHFT 470 NBL = NVIRT*(NL-1) + NB 471 AIBJCK = AIBJCK - HALF*L3AM(NAI,NBL,NCK) * FOCKA(NJ,NL) 472 END DO 473 474 E3AM(NAI,NBJ,NCK) = E3AM(NAI,NBJ,NCK) + AIBJCK 475 E3AM(NAI,NCK,NBJ) = E3AM(NAI,NCK,NBJ) + AIBJCK 476 E3AM(NBJ,NAI,NCK) = E3AM(NBJ,NAI,NCK) + AIBJCK 477 E3AM(NCK,NAI,NBJ) = E3AM(NCK,NAI,NBJ) + AIBJCK 478 E3AM(NBJ,NCK,NAI) = E3AM(NBJ,NCK,NAI) + AIBJCK 479 E3AM(NCK,NBJ,NAI) = E3AM(NCK,NBJ,NAI) + AIBJCK 480 481 END DO 482 END DO 483 END DO 484 END DO 485 END DO 486 END DO 487C 488C------------------------------------------------------ 489C Get rid of amplitudes that are not allowed. 490C------------------------------------------------------ 491 492 CALL CCSDT_CLEAN_T3(E3AM,NT1AMX,NVIRT,NRHFT) 493 494 RETURN 495 END 496*---------------------------------------------------------------------* 497* END OF SUBROUTINE CCSDT_E3AM * 498*---------------------------------------------------------------------* 499*=====================================================================* 500 SUBROUTINE CCSDT_E2AM(E2AM,L3AM,T2AM,FOCKA) 501*---------------------------------------------------------------------* 502* 503* Purpose: compute triples correction to doubles Eta vector 504* 505* Eta_nu_2(L3) = <L_3|[[A,T^0_2],tau_nu_2]|HF> 506* 507* Written by Christof Haettig, April 2002 508*=====================================================================* 509#if defined (IMPLICIT_NONE) 510 IMPLICIT NONE 511#else 512# include "implicit.h" 513#endif 514#include "priunit.h" 515#include "ccsdinp.h" 516#include "maxorb.h" 517#include "ccsdsym.h" 518#include "ccorb.h" 519 520#if defined (SYS_CRAY) 521 REAL E2AM(NT1AMX,NT1AMX) 522 REAL L3AM(NT1AMX,NT1AMX,NT1AMX) 523 REAL T2AM(NT1AMX,NT1AMX) 524 REAL FOCKA(NORBT,NORBT) 525 REAL CONTRIB, HALF 526#else 527 DOUBLE PRECISION E2AM(NT1AMX,NT1AMX) 528 DOUBLE PRECISION L3AM(NT1AMX,NT1AMX,NT1AMX) 529 DOUBLE PRECISION T2AM(NT1AMX,NT1AMX) 530 DOUBLE PRECISION FOCKA(NORBT,NORBT) 531 DOUBLE PRECISION CONTRIB, HALF 532#endif 533 PARAMETER (HALF = 0.5D0) 534 535 INTEGER NAI, NCK, NJ, NL, NB, ND, NBJ, NBL, NDL, NDJ 536 537 DO NAI = 1, NT1AMX 538 DO NCK = 1, NT1AMX 539 DO NJ = 1, NRHFT 540 DO NL = 1, NRHFT 541 DO NB = 1, NVIRT 542 DO ND = 1, NVIRT 543 NBJ = NVIRT*(NJ-1) + NB 544 NBL = NVIRT*(NL-1) + NB 545 NDL = NVIRT*(NL-1) + ND 546 NDJ = NVIRT*(NJ-1) + ND 547 548 CONTRIB = 549 & ( L3AM(NAI,NCK,NBL)*T2AM(NCK,NDL)*FOCKA(NJ,NRHFT+ND) + 550 & L3AM(NAI,NCK,NDJ)*T2AM(NCK,NDL)*FOCKA(NL,NRHFT+NB) ) 551 552 E2AM(NAI,NBJ) = E2AM(NAI,NBJ) - CONTRIB 553 E2AM(NBJ,NAI) = E2AM(NBJ,NAI) - CONTRIB 554 555 END DO 556 END DO 557 END DO 558 END DO 559 END DO 560 END DO 561 562 563 RETURN 564 END 565*---------------------------------------------------------------------* 566* END OF SUBROUTINE CCSDT_E2AM * 567*---------------------------------------------------------------------* 568*=====================================================================* 569 SUBROUTINE CCSDT_E2AXC(XMMAT,T2CAM,DIA,FOCKA) 570*---------------------------------------------------------------------* 571* 572* Purpose: noddy code for M matrix intermediate used in the `real` 573* code to compute the contraction Eta x R 574* 575* Written by Christof Haettig, May 2002 576*=====================================================================* 577#if defined (IMPLICIT_NONE) 578 IMPLICIT NONE 579#else 580# include "implicit.h" 581#endif 582#include "priunit.h" 583#include "ccsdinp.h" 584#include "maxorb.h" 585#include "ccsdsym.h" 586#include "ccorb.h" 587 588#if defined (SYS_CRAY) 589 REAL XMMAT(NRHFT,NRHFT,NT1AMX) 590 REAL T2CAM(NT1AMX,NT1AMX) 591 REAL FOCKA(NORBT,NORBT) 592 REAL DIA(NT1AMX), CONTRIB 593#else 594 DOUBLE PRECISION XMMAT(NRHFT,NRHFT,NT1AMX) 595 DOUBLE PRECISION T2CAM(NT1AMX,NT1AMX) 596 DOUBLE PRECISION FOCKA(NORBT,NORBT) 597 DOUBLE PRECISION DIA(NT1AMX), CONTRIB 598#endif 599 600 INTEGER NFN, NM, NI, NA, NAI, NAM 601 602 CALL DZERO(DIA,NT1AMX) 603 604 WRITE(LUPRI,*) 'M complete matrix:' 605 CALL OUTPUT(XMMAT,1,NRHFT*NRHFT,1,NT1AMX,NRHFT*NRHFT,NT1AMX, 606 & 1,LUPRI) 607 608 WRITE(LUPRI,*) 'response amplitudes:' 609 CALL OUTPUT(T2CAM,1,NT1AMX,1,NT1AMX,NT1AMX,NT1AMX, 610 & 1,LUPRI) 611 612 DO NFN = 1, NT1AMX 613 614 DO NM = 1, NRHFT 615 DO NI = 1, NRHFT 616 DO NA = 1, NVIRT 617 NAI = NVIRT*(NI-1) + NA 618 NAM = NVIRT*(NM-1) + NA 619 DIA(NAI) = DIA(NAI) + 620 & XMMAT(NI,NM,NFN) * T2CAM(NAM,NFN) 621 END DO 622 END DO 623 END DO 624 625 write(lupri,*) 'NFN:',NFN 626 write(lupri,*) 'amplitudes:' 627 call output(t2cam(1,nfn),1,nvirt,1,nrhft,nvirt,nrhft,1,lupri) 628 write(lupri,*) 'm matrix:' 629 call output(xmmat(1,1,nfn),1,nrhft,1,nrhft,nrhft,nrhft,1,lupri) 630 write(lupri,*) 'density:' 631 call output(dia,1,nvirt,1,nrhft,nvirt,nrhft,1,lupri) 632 633 END DO 634 635 WRITE(LUPRI,*) 'CCSDT_E2AXC> Density:' 636 CALL OUTPUT(DIA,1,NVIRT,1,NRHFT,NVIRT,NRHFT,1,LUPRI) 637 638 CONTRIB = 0.0D0 639 DO NI = 1, NRHFT 640 DO NA = 1, NVIRT 641 NAI = NVIRT*(NI-1) + NA 642 CONTRIB = CONTRIB + DIA(NAI) * FOCKA(NI,NRHFT+NA) 643 END DO 644 END DO 645 646 WRITE(LUPRI,*) 'CCSDT_E2AXC> Contribution:',CONTRIB 647 648 RETURN 649 END 650*---------------------------------------------------------------------* 651* END OF SUBROUTINE CCSDT_E2AXC * 652*---------------------------------------------------------------------* 653*=====================================================================* 654 SUBROUTINE CCSDT_E1AM(E1AM,L3AM,T3AM,FOCKA) 655*---------------------------------------------------------------------* 656* 657* Purpose: compute triples correction to singles Eta vector 658* 659* Eta_nu_1(L3,T3) = <L_3|[[A,T^0_3],tau_nu_1]|HF> 660* 661* Written by Christof Haettig, April 2002 662*=====================================================================* 663#if defined (IMPLICIT_NONE) 664 IMPLICIT NONE 665#else 666# include "implicit.h" 667#endif 668#include "priunit.h" 669#include "ccsdinp.h" 670#include "maxorb.h" 671#include "ccsdsym.h" 672#include "ccorb.h" 673 674 LOGICAL LOCDBG 675 PARAMETER (LOCDBG = .FALSE.) 676 677#if defined (SYS_CRAY) 678 REAL E1AM(NT1AMX), DDOT 679 REAL L3AM(NT1AMX,NT1AMX,NT1AMX) 680 REAL T3AM(NT1AMX,NT1AMX,NT1AMX) 681 REAL FOCKA(NORBT,NORBT) 682#else 683 DOUBLE PRECISION E1AM(NT1AMX), DDOT 684 DOUBLE PRECISION L3AM(NT1AMX,NT1AMX,NT1AMX) 685 DOUBLE PRECISION T3AM(NT1AMX,NT1AMX,NT1AMX) 686 DOUBLE PRECISION FOCKA(NORBT,NORBT) 687#endif 688 689 INTEGER NAI, NBJ, NK, NL, NCK, NC, ND, NCL, NDL, NDK 690 691 IF (LOCDBG) THEN 692 WRITE(LUPRI,*) 'CCSDT_E1AM> Norm^(E1AM) on input:', 693 & DDOT(NT1AMX,E1AM,1,E1AM,1) 694 WRITE(LUPRI,*) 'CCSDT_E1AM> Norm^(L3AM) on input:', 695 & DDOT(NT1AMX*NT1AMX*NT1AMX,L3AM,1,L3AM,1) 696 WRITE(LUPRI,*) 'CCSDT_E1AM> Norm^(T3AM) on input:', 697 & DDOT(NT1AMX*NT1AMX*NT1AMX,T3AM,1,T3AM,1) 698 WRITE(LUPRI,*) 'CCSDT_E1AM> Norm^(FOCKA) on input:', 699 & DDOT(NORBT*NORBT,FOCKA,1,FOCKA,1) 700 END IF 701 702 DO NAI = 1, NT1AMX 703 DO NBJ = 1, NT1AMX 704 DO NK = 1, NRHFT 705 DO NL = 1, NRHFT 706 DO NC = 1, NVIRT 707 DO ND = 1, NVIRT 708 NCL = NVIRT*(NL-1) + NC 709 NCK = NVIRT*(NK-1) + NC 710 NDL = NVIRT*(NL-1) + ND 711 NDK = NVIRT*(NK-1) + ND 712 713 E1AM(NDL) = E1AM(NDL) - 0.5D0 * 714 & ( L3AM(NAI,NBJ,NDK) * FOCKA(NL,NRHFT+NC) + 715 & L3AM(NAI,NBJ,NCL) * FOCKA(NK,NRHFT+ND) ) * 716 & T3AM(NAI,NBJ,NCK) 717 718 END DO 719 END DO 720 END DO 721 END DO 722 END DO 723 END DO 724 725 IF (LOCDBG) THEN 726 WRITE(LUPRI,*) 'CCSDT_E1AM> E1AM on output:' 727 CALL CC_PRP(E1AM,E1AM,1,1,0) 728 END IF 729 730 RETURN 731 END 732*---------------------------------------------------------------------* 733* END OF SUBROUTINE CCSDT_E1AM * 734*---------------------------------------------------------------------* 735*=====================================================================* 736 SUBROUTINE CCSDT_ETA_FD(T1AM,T2AM,T3AM,L3AM,L2AM,ETA1,ETA2,ETA3, 737 & FOCKA,FOCKA_AO,SCR1,WORK,LWORK) 738*---------------------------------------------------------------------* 739* 740* Purpose: Construct Eta vector by finite difference on Xi vector 741* 742* Written by Christof Haettig, April 2002 743*=====================================================================* 744#if defined (IMPLICIT_NONE) 745 IMPLICIT NONE 746#else 747# include "implicit.h" 748#endif 749#include "priunit.h" 750#include "ccsdinp.h" 751#include "maxorb.h" 752#include "ccsdsym.h" 753#include "ccorb.h" 754 755#if defined (SYS_CRAY) 756 REAL DELTA, ZERO, HALF, SIXTH 757#else 758 DOUBLE PRECISION DELTA, ZERO, HALF, SIXTH 759#endif 760 PARAMETER(DELTA=1.0D-6, ZERO=0.0D0, HALF=0.5D0, SIXTH=1.0D0/6.0D0) 761 762 INTEGER LWORK 763 764#if defined (SYS_CRAY) 765 REAL WORK(LWORK), SCR1(NT1AMX), 766 & T1AM(NT1AMX), T2AM(NT1AMX,NT1AMX), T3AM(NT1AMX,NT1AMX,NT1AMX), 767 & L2AM(NT1AMX,NT1AMX), L3AM(NT1AMX,NT1AMX,NT1AMX), 768 & ETA1(NT1AMX), ETA2(NT1AMX,NT1AMX), ETA3(NT1AMX,NT1AMX,NT1AMX), 769 & FOCKA(NORBT,NORBT), FOCKA_AO(NORBT,NORBT), 770 & TAIBJ, TAIBJCK, DDOT, EAIBJ 771#else 772 DOUBLE PRECISION WORK(LWORK), SCR1(NT1AMX), 773 & T1AM(NT1AMX), T2AM(NT1AMX,NT1AMX), T3AM(NT1AMX,NT1AMX,NT1AMX), 774 & L2AM(NT1AMX,NT1AMX), L3AM(NT1AMX,NT1AMX,NT1AMX), 775 & ETA1(NT1AMX), ETA2(NT1AMX,NT1AMX), ETA3(NT1AMX,NT1AMX,NT1AMX), 776 & FOCKA(NORBT,NORBT), FOCKA_AO(NORBT,NORBT), 777 & TAIBJ, TAIBJCK, DDOT, EAIBJCK, EAIBJ 778#endif 779 780 INTEGER KT2AM, KXI2, KXI3, NAI, NBJ, NCK, KEND1, LWRK1, KT1AM, 781 & KFOCKA, KLAMP0, KLAMH0, KFOCKD, NA, NI, NB, NBI, NC, 782 & NCI, NJ, NAJ, NK, NAK, KL3AM 783 784 785 CALL QUIT('CCSDT_ETA_FD needs to be addapted for intermediates.') 786*---------------------------------------------------------------------* 787* Allocations: 788*---------------------------------------------------------------------* 789 KL3AM = 1 790 KFOCKD= KL3AM + NT1AMX*NT1AMX*NT1AMX 791 KFOCKA= KFOCKD + NORBT 792 KLAMP0= KFOCKA + NORBT*NORBT 793 KLAMH0= KLAMP0 + NLAMDT 794 KT1AM = KLAMH0 + NLAMDT 795 KT2AM = KT1AM + NT1AMX 796 KXI2 = KT2AM + NT1AMX*NT1AMX 797 KXI3 = KXI2 + NT1AMX*NT1AMX 798 KEND1 = KXI3 + NT1AMX*NT1AMX*NT1AMX 799 LWRK1 = LWORK - KEND1 800 IF (LWRK1 .LT. 0) THEN 801 CALL QUIT('Insufficient space in CCSDT_ETA_FD') 802 ENDIF 803 804 DO I = 1, NRHFT 805 WORK(KFOCKD-1+I) = 0.0d0 806 ENDDO 807 DO I = 1, NVIRT 808 WORK(KFOCKD-1+NRHFT+I) = 1.0d0/3.0d0 809 ENDDO 810 811 ! turn sign of T3, because Xi vector routines work with -T3 812 CALL DSCAL(NT1AMX*NT1AMX*NT1AMX,-1.0D0,T3AM,1) 813 814C WRITE(LUPRI,*) 'in finite difference norm^2(L3AM):', 815C & DDOT(NT1AMX*NT1AMX*NT1AMX,L3AM,1,L3AM,1) 816C WRITE(LUPRI,*) 'in finite difference norm^2(T3AM):', 817C & DDOT(NT1AMX*NT1AMX*NT1AMX,T3AM,1,T3AM,1) 818C WRITE(LUPRI,*) 'in finite difference norm^2(FOCKA):', 819C & DDOT(NORBT*NORBT,FOCKA,1,FOCKA,1) 820 821C write(lupri,*) 'in finite difference focka:' 822C call output(focka,1,norbt,1,norbt,norbt,norbt,1,lupri) 823C write(lupri,*) 'in finite difference T3AM:' 824C call output(t3am,1,NT1AMX*NT1AMX,1,NT1AMX, 825C & NT1AMX*NT1AMX,NT1AMX,1,lupri) 826*---------------------------------------------------------------------* 827* take the derivative of Xksi_3 w.r.t. T1 828* use, that Xksi_3 is linear in T1 829*---------------------------------------------------------------------* 830 CALL DCOPY(NT1AMX,T1AM,1,WORK(KT1AM),1) 831 CALL DZERO(WORK(KT2AM),NT1AMX*NT1AMX) 832 CALL DZERO(ETA1,NT1AMX) 833 834 DO NAI = 1, NT1AMX 835 836 CALL DZERO(WORK(KT1AM),NT1AMX) 837 WORK(KT1AM-1+NAI) = T1AM(NAI) - DELTA*HALF 838 839 CALL LAMMAT(WORK(KLAMP0),WORK(KLAMH0),WORK(KT1AM), 840 & WORK(KEND1),LWRK1) 841 CALL DCOPY(NORBT*NORBT,FOCKA_AO,1,WORK(KFOCKA),1) 842 CALL CC_FCKMO(WORK(KFOCKA),WORK(KLAMP0),WORK(KLAMH0), 843 & WORK(KEND1),LWRK1,1,1,1) 844 845 CALL DZERO(WORK(KXI3),NT1AMX*NT1AMX*NT1AMX) 846 CALL CCSDT_XKSI3(WORK(KXI3),WORK(KFOCKA),T3AM,T2AM) 847 848 ETA1(NAI) = ETA1(NAI) - 849 & SIXTH*DDOT(NT1AMX*NT1AMX*NT1AMX,WORK(KXI3),1,L3AM,1) 850 851 WORK(KT1AM-1+NAI) = T1AM(NAI) + DELTA*HALF 852 853 CALL LAMMAT(WORK(KLAMP0),WORK(KLAMH0),WORK(KT1AM), 854 & WORK(KEND1),LWRK1) 855 CALL DCOPY(NORBT*NORBT,FOCKA_AO,1,WORK(KFOCKA),1) 856 CALL CC_FCKMO(WORK(KFOCKA),WORK(KLAMP0),WORK(KLAMH0), 857 & WORK(KEND1),LWRK1,1,1,1) 858 859 CALL DZERO(WORK(KXI3),NT1AMX*NT1AMX*NT1AMX) 860 CALL CCSDT_XKSI3(WORK(KXI3),WORK(KFOCKA),T3AM,T2AM) 861 862 ETA1(NAI) = ETA1(NAI) + 863 & SIXTH*DDOT(NT1AMX*NT1AMX*NT1AMX,WORK(KXI3),1,L3AM,1) 864 865 WORK(KT1AM-1+NAI) = T1AM(NAI) 866 END DO 867 868 869*---------------------------------------------------------------------* 870* take the derivative of Xksi2 & Xksi_3 w.r.t. T3 871*---------------------------------------------------------------------* 872 CALL DZERO(WORK(KT2AM),NT1AMX*NT1AMX) 873 CALL DZERO(ETA3,NT1AMX*NT1AMX*NT1AMX) 874 875C CALL DZERO(WORK(KL3AM),NT1AMX*NT1AMX*NT1AMX) 876C WORK(KL3AM-1+13) = L3AM(13,1,1) 877 CALL DCOPY(NT1AMX*NT1AMX*NT1AMX,L3AM,1,WORK(KL3AM),1) 878 879 if (.true.) then 880 881 DO NAI = 1, NT1AMX 882 DO NBJ = 1, NT1AMX 883 DO NCK = 1, NT1AMX 884 TAIBJCK = T3AM(NAI,NBJ,NCK) 885 T3AM(NAI,NBJ,NCK) = TAIBJCK - DELTA*HALF 886 887 CALL DZERO(WORK(KXI2),NT1AMX*NT1AMX) 888 CALL CCSDT_XKSI2(WORK(KXI2),FOCKA,T3AM) 889 890 CALL DZERO(WORK(KXI3),NT1AMX*NT1AMX*NT1AMX) 891 CALL CCSDT_XKSI3(WORK(KXI3),FOCKA,T3AM,WORK(KT2AM)) 892 893 EAIBJCK = + 894 & HALF*DDOT(NT1AMX*NT1AMX,WORK(KXI2),1,L2AM,1) + 895 & SIXTH*DDOT(NT1AMX*NT1AMX*NT1AMX,WORK(KXI3),1,L3AM,1) 896 897 T3AM(NAI,NBJ,NCK) = TAIBJCK + DELTA*HALF 898 899 CALL DZERO(WORK(KXI2),NT1AMX*NT1AMX) 900 CALL CCSDT_XKSI2(WORK(KXI2),FOCKA,T3AM) 901 902 CALL DZERO(WORK(KXI3),NT1AMX*NT1AMX*NT1AMX) 903 CALL CCSDT_XKSI3(WORK(KXI3),FOCKA,T3AM,WORK(KT2AM)) 904 905 EAIBJCK = EAIBJCK - 906 & HALF*DDOT(NT1AMX*NT1AMX,WORK(KXI2),1,L2AM,1) - 907 & SIXTH*DDOT(NT1AMX*NT1AMX*NT1AMX,WORK(KXI3),1,L3AM,1) 908 909 ETA3(NAI,NBJ,NCK) = ETA3(NAI,NBJ,NCK) + EAIBJCK 910 ETA3(NAI,NCK,NBJ) = ETA3(NAI,NCK,NBJ) + EAIBJCK 911 ETA3(NBJ,NAI,NCK) = ETA3(NBJ,NAI,NCK) + EAIBJCK 912 ETA3(NCK,NAI,NBJ) = ETA3(NCK,NAI,NBJ) + EAIBJCK 913 ETA3(NBJ,NCK,NAI) = ETA3(NBJ,NCK,NAI) + EAIBJCK 914 ETA3(NCK,NBJ,NAI) = ETA3(NCK,NBJ,NAI) + EAIBJCK 915 916 T3AM(NAI,NBJ,NCK) = TAIBJCK 917 END DO 918 END DO 919 END DO 920 end if 921 922*---------------------------------------------------------------------* 923* take the derivative of Xksi_3 w.r.t. T2 924*---------------------------------------------------------------------* 925 CALL DCOPY(NT1AMX*NT1AMX,T2AM,1,WORK(KT2AM),1) 926 CALL DZERO(ETA2,NT1AMX*NT1AMX) 927 928 if (.true.) then 929 930 DO NAI = 1, NT1AMX 931 DO NBJ = 1, NT1AMX 932 TAIBJ = T2AM(NAI,NBJ) 933 WORK(KT2AM-1+(NBJ-1)*NT1AMX+NAI) = TAIBJ - DELTA*HALF 934 935 CALL DZERO(WORK(KXI3),NT1AMX*NT1AMX*NT1AMX) 936 CALL CCSDT_XKSI3(WORK(KXI3),FOCKA,T3AM,WORK(KT2AM)) 937 938 EAIBJ = - 939 & SIXTH*DDOT(NT1AMX*NT1AMX*NT1AMX,WORK(KXI3),1,L3AM,1) 940 941 942 WORK(KT2AM-1+(NBJ-1)*NT1AMX+NAI) = TAIBJ + DELTA*HALF 943 944 CALL DZERO(WORK(KXI3),NT1AMX*NT1AMX*NT1AMX) 945 CALL CCSDT_XKSI3(WORK(KXI3),FOCKA,T3AM,WORK(KT2AM)) 946 947 EAIBJ = EAIBJ + 948 & SIXTH*DDOT(NT1AMX*NT1AMX*NT1AMX,WORK(KXI3),1,L3AM,1) 949 950 ETA2(NAI,NBJ) = ETA2(NAI,NBJ) + EAIBJ 951 ETA2(NBJ,NAI) = ETA2(NBJ,NAI) + EAIBJ 952 953 WORK(KT2AM-1+(NBJ-1)*NT1AMX+NAI) = TAIBJ 954 END DO 955 END DO 956 957 end if 958 959*---------------------------------------------------------------------* 960* scale result with 1/DELTA and remove forbidden triples elements 961*---------------------------------------------------------------------* 962 CALL DSCAL(NT1AMX, 1.0D0/DELTA,ETA1,1) 963 CALL DSCAL(NT1AMX*NT1AMX, 1.0D0/DELTA,ETA2,1) 964 CALL DSCAL(NT1AMX*NT1AMX*NT1AMX,1.0D0/DELTA,ETA3,1) 965 966 CALL CCSDT_CLEAN_T3(ETA3,NT1AMX,NVIRT,NRHFT) 967 968 ! restor correct sign of T3 969 CALL DSCAL(NT1AMX*NT1AMX*NT1AMX,-1.0D0,T3AM,1) 970 971 RETURN 972 END 973*---------------------------------------------------------------------* 974* END OF SUBROUTINE CCSDT_ETA_FD * 975*---------------------------------------------------------------------* 976