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*---------------------------------------------------------------------* 20C /* Deck cc_lhtr_noddy */ 21*---------------------------------------------------------------------* 22 SUBROUTINE CC_LHTR_NODDY(FREQ,OMEGA1,OMEGA2,T1AM,T2AM,C1AM,C2AM, 23 * XLAMDP,XLAMDH,WORK,LWORK,IOPT) 24*=====================================================================* 25C 26C Purpose: 27C 28C Calculate the iterative triples corrections to the CCSD 29C equations. 30C 31C NB! The T2 amplitudes are assumed to be a full square. 32C 33C 34C NB! It is assumed that the vectors are allocated in the following 35C order: 36C T1AM(*), OMEGA1(*), OMEGA2(*), T2AM(*), SCR(*), WRK(*). 37C 38C IOPT = 0 No storage of T3, L3 39C IOPT = 1 Store T3 and L3 in 'T3AMPL' and 'L3AMPL' 40C 41C 42*---------------------------------------------------------------------* 43#include "implicit.h" 44#include "priunit.h" 45#include "dummy.h" 46#include "maxorb.h" 47#include "maxash.h" 48#include "aovec.h" 49#include "inforb.h" 50#include "infind.h" 51#include "blocks.h" 52#include "inftap.h" 53#include "ccsdinp.h" 54#include "ccsdsym.h" 55#include "ccsdio.h" 56#include "mxcent.h" 57#include "eritap.h" 58#include "ccfield.h" 59#include "ccnoddy.h" 60C 61 LOGICAL LOCDBG 62 PARAMETER(LOCDBG = .FALSE.) 63C 64 PARAMETER (ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, TWO = 2.0D0) 65C 66 DIMENSION OMEGA1(*),OMEGA2(*),T1AM(*),T2AM(*), C1AM(*), C2AM(*) 67 DIMENSION XLAMDP(*),XLAMDH(*),WORK(LWORK) 68C 69C INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 70C 71 IF (NSYM .NE. 1) THEN 72 CALL QUIT('No symmetry in this part yet') 73 ENDIF 74C 75C------------------------ 76C Dynamic Allocation. 77C------------------------ 78C 79 KSCR1 = 1 80 KFOCK = KSCR1 + NT1AMX 81 KFOCKD = KFOCK + N2BST(ISYMOP) 82 KEND1 = KFOCKD + NORBT 83 84 IF ((NONHF) .AND. (NFIELD .GT. 0)) THEN 85 KONEP = KEND1 86 KONEPAO= KONEP + N2BST(ISYMOP) 87 KEND1 = KONEPAO+N2BST(ISYMOP) 88 LWRK1 = LWORK - KEND1 89 ENDIF 90 91 KINT1S = KEND1 92 KINT2S = KINT1S + NT1AMX*NVIRT*NVIRT 93 KEND1 = KINT2S + NRHFT*NRHFT*NT1AMX 94 95 KIOVVO = KEND1 96 KIOOVV = KIOVVO + NRHFT*NVIRT*NVIRT*NRHFT 97 KIOOOO = KIOOVV + NRHFT*NRHFT*NVIRT*NVIRT 98 KIVVVV = KIOOOO + NRHFT*NRHFT*NRHFT*NRHFT 99 KEND1 = KIVVVV + NVIRT*NVIRT*NVIRT*NVIRT 100 101 KL3AM = KEND1 102 KEND1 = KL3AM + NT1AMX*NT1AMX*NT1AMX 103 104 ! what is above has to be kept until the end... 105 ! everything below might be overwritten in CC_LHPART_NODDY 106 KEND1A = KEND1 107 LWRK1A = LWORK - KEND1A 108 109 KINT1T = KEND1 110 KINT2T = KINT1T + NT1AMX*NVIRT*NVIRT 111 KEND1 = KINT2T + NRHFT*NRHFT*NT1AMX 112 113 KT3AM = KEND1 114 KEND1 = KT3AM + NT1AMX*NT1AMX*NT1AMX 115 116 KOME1 = KEND1 117 KOME2 = KOME1 + NT1AMX 118 KEND1 = KOME2 + NT1AMX*NT1AMX 119 120 KIAJB = KEND1 121 KYIAJB = KIAJB + NT1AMX*NT1AMX 122 KEND1 = KYIAJB + NT1AMX*NT1AMX 123 124 LWRK1 = LWORK - KEND1 125 IF (LWRK1 .LT. 0) THEN 126 CALL QUIT('Insufficient space in CC_LHTR_NODDY') 127 ENDIF 128C 129C-------------------------------- 130C Initialize integral arrays. 131C-------------------------------- 132C 133 CALL DZERO(WORK(KL3AM),NT1AMX*NT1AMX*NT1AMX) 134 135 CALL CCSDT_READ_NODDY(.TRUE.,WORK(KFOCKD),WORK(KFOCK), 136 & WORK(KONEP), WORK(KONEPAO), 137 & .TRUE.,WORK(KIAJB), WORK(KYIAJB), 138 & .TRUE.,WORK(KINT1S),WORK(KINT2S), 139 & .TRUE.,WORK(KINT1T),WORK(KINT2T), 140 & .TRUE.,WORK(KIOVVO),WORK(KIOOVV), 141 & WORK(KIOOOO),WORK(KIVVVV), 142 & NORBT,NLAMDT,NRHFT,NVIRT,NT1AMX) 143C 144C------------------------------------- 145C Get the T3 amplitudes. 146C------------------------------------- 147 LUTEMP = -1 148 CALL GPOPEN(LUTEMP,FILNODT30,'UNKNOWN',' ','UNFORMATTED', 149 & IDUMMY,.FALSE.) 150 READ(LUTEMP) (WORK(KT3AM+I-1), I=1,NT1AMX*NT1AMX*NT1AMX) 151 CALL GPCLOSE(LUTEMP,'KEEP') 152C 153 IF (IOPT .EQ. 1 .OR. (NONHF .AND. NFIELD.GT.0)) THEN 154 LUSIFC = -1 155 CALL GPOPEN(LUSIFC,'T3AMPL','UNKNOWN',' ','UNFORMATTED',IDUMMY, 156 & .FALSE.) 157 REWIND LUSIFC 158 WRITE (LUSIFC) (WORK(KT3AM-1+I), I=1,NT1AMX*NT1AMX*NT1AMX) 159 CALL GPCLOSE(LUSIFC,'KEEP') 160 ENDIF 161C 162C------------------------------------ 163C Calculate the CC3 corrections. 164C------------------------------------ 165C 166 CALL DZERO(WORK(KOME1),NT1AMX) 167C 168 CALL T3_LEFT1(WORK(KOME1),WORK(KYIAJB),WORK(KIAJB), 169 * C2AM,WORK(KT3AM)) 170C 171 CALL T3_LEFT2(WORK(KOME1),WORK(KYIAJB),WORK(KIAJB), 172 * C2AM,WORK(KT3AM)) 173C 174 CALL T3_LEFT3(WORK(KOME1),WORK(KIAJB),C2AM,WORK(KT3AM)) 175C 176 DO I = 1,NT1AMX 177 OMEGA1(I) = OMEGA1(I) + WORK(KOME1+I-1) 178 ENDDO 179C 180C--------------------------------------------------------------- 181C Calculate rhs vector for triples left amplitude equations: 182C--------------------------------------------------------------- 183C 184 CALL DZERO(WORK(KL3AM),NT1AMX*NT1AMX*NT1AMX) 185 186 CALL CCSDT_L3AM_R(WORK(KL3AM),FREQ,WORK(KINT1T),WORK(KINT2T), 187 * WORK(KIAJB),WORK(KFOCK),C1AM,C2AM, 188 * WORK(KSCR1),WORK(KFOCKD),.FALSE.) 189C 190C---------------------------------------------------------- 191C Solve triples equations and partition solution vector 192C into the singles and doubles results vectors: 193C---------------------------------------------------------- 194C 195 CALL CC_LHPART_NODDY(OMEGA1,OMEGA2,WORK(KL3AM),FREQ, 196 & WORK(KFOCKD),WORK(KONEP), 197 & WORK(KIOOOO),WORK(KIOVVO), 198 & WORK(KIOOVV),WORK(KIVVVV), 199 & T2AM,WORK(KINT1S),WORK(KINT2S), 200 & WORK(KEND1A),LWRK1A) 201 202 203 IF (IOPT .EQ. 1) THEN 204 LUSIFC = -1 205 CALL GPOPEN(LUSIFC,'L3AMPL','UNKNOWN',' ','UNFORMATTED',IDUMMY, 206 & .FALSE.) 207 REWIND LUSIFC 208 WRITE (LUSIFC) (WORK(KL3AM-1+I), I=1,NT1AMX*NT1AMX*NT1AMX) 209 CALL GPCLOSE(LUSIFC,'KEEP') 210 ENDIF 211 212 IF (LOCDBG) THEN 213 WRITE(LUPRI,*) 'Result vector after CC_LHTR_NODDY:' 214 CALL CC_PRP(OMEGA1,OMEGA2,1,1,1) 215 END IF 216 217 RETURN 218 END 219*---------------------------------------------------------------------* 220C /* Deck t3_left1 */ 221*---------------------------------------------------------------------* 222 SUBROUTINE T3_LEFT1(OMEGA1,YIAJB,XIAJB,C2AM,T3AM) 223*=====================================================================* 224C 225C 226C 227*---------------------------------------------------------------------* 228#include "implicit.h" 229#include "priunit.h" 230#include "inforb.h" 231#include "ccsdinp.h" 232#include "ccsdsym.h" 233 DIMENSION OMEGA1(NT1AMX), YIAJB(NT1AMX,NT1AMX) 234 DIMENSION XIAJB(NT1AMX,NT1AMX) 235 DIMENSION C2AM(NT1AMX,NT1AMX), T3AM(NT1AMX,NT1AMX,NT1AMX) 236C 237C INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 238C 239 DO A = 1, NVIRT 240 DO I = 1, NRHFT 241 NAI = NVIRT*(I-1) + A 242 DO D = 1, NVIRT 243 DO L = 1, NRHFT 244 NDI = NVIRT*(I-1) + D 245 NDL = NVIRT*(L-1) + D 246 NAL = NVIRT*(L-1) + A 247 DO E = 1, NVIRT 248 DO M = 1, NRHFT 249 NEM = NVIRT*(M-1) + E 250 DO F = 1, NVIRT 251 DO N = 1, NRHFT 252 NEN = NVIRT*(N-1) + E 253 NFN = NVIRT*(N-1) + F 254 NFM = NVIRT*(M-1) + F 255C 256 OMEGA1(NAI) = OMEGA1(NAI) 257 * + T3AM(NDL,NEM,NFN) 258 * * C2AM(NFM,NDI)*YIAJB(NEN,NAL) 259 * - T3AM(NDL,NEN,NFM) 260 * * C2AM(NFM,NDI)*XIAJB(NEN,NAL) 261 ENDDO 262 ENDDO 263 ENDDO 264 ENDDO 265 ENDDO 266 ENDDO 267 ENDDO 268 ENDDO 269C 270 RETURN 271 END 272*---------------------------------------------------------------------* 273C /* Deck t3_left2 */ 274*---------------------------------------------------------------------* 275 SUBROUTINE T3_LEFT2(OMEGA1,YIAJB,XIAJB,C2AM,T3AM) 276*=====================================================================* 277C 278C 279C 280*---------------------------------------------------------------------* 281#include "implicit.h" 282#include "priunit.h" 283#include "inforb.h" 284#include "ccsdinp.h" 285#include "ccsdsym.h" 286 DIMENSION OMEGA1(NT1AMX), YIAJB(NT1AMX,NT1AMX) 287 DIMENSION XIAJB(NT1AMX,NT1AMX) 288 DIMENSION C2AM(NT1AMX,NT1AMX), T3AM(NT1AMX,NT1AMX,NT1AMX) 289C 290C INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 291C 292 DO A = 1, NVIRT 293 DO I = 1, NRHFT 294 NAI = NVIRT*(I-1) + A 295 DO D = 1, NVIRT 296 DO L = 1, NRHFT 297 NDL = NVIRT*(L-1) + D 298 DO E = 1, NVIRT 299 DO M = 1, NRHFT 300 NAM = NVIRT*(M-1) + A 301 NEI = NVIRT*(I-1) + E 302 NEM = NVIRT*(M-1) + E 303 DO F = 1, NVIRT 304 DO N = 1, NRHFT 305 NDN = NVIRT*(N-1) + D 306 NFL = NVIRT*(L-1) + F 307 NFN = NVIRT*(N-1) + F 308C 309 OMEGA1(NAI) = OMEGA1(NAI) 310 * + T3AM(NDL,NEM,NFN) 311 * * C2AM(NAM,NDN)*YIAJB(NEI,NFL) 312 * - T3AM(NDN,NEM,NFL) 313 * * C2AM(NAM,NDN)*XIAJB(NEI,NFL) 314 ENDDO 315 ENDDO 316 ENDDO 317 ENDDO 318 ENDDO 319 ENDDO 320 ENDDO 321 ENDDO 322C 323 RETURN 324 END 325*---------------------------------------------------------------------* 326C /* Deck t3_left3 */ 327*---------------------------------------------------------------------* 328 SUBROUTINE T3_LEFT3(OMEGA1,XIAJB,C2AM,T3AM) 329*=====================================================================* 330C 331C Note : XIAJB is L and not g. 332C 333C 334*---------------------------------------------------------------------* 335#include "implicit.h" 336#include "priunit.h" 337#include "inforb.h" 338#include "ccsdinp.h" 339#include "ccsdsym.h" 340 DIMENSION OMEGA1(NT1AMX), XIAJB(NT1AMX,NT1AMX) 341 DIMENSION C2AM(NT1AMX,NT1AMX), T3AM(NT1AMX,NT1AMX,NT1AMX) 342C 343C INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 344C 345 DO A = 1, NVIRT 346 DO I = 1, NRHFT 347 NAI = NVIRT*(I-1) + A 348 DO D = 1, NVIRT 349 DO L = 1, NRHFT 350 NDL = NVIRT*(L-1) + D 351 DO E = 1, NVIRT 352 DO M = 1, NRHFT 353 NEM = NVIRT*(M-1) + E 354 DO F = 1, NVIRT 355 DO N = 1, NRHFT 356 NEN = NVIRT*(N-1) + E 357 NFM = NVIRT*(M-1) + F 358 NFN = NVIRT*(N-1) + F 359C 360 OMEGA1(NAI) = OMEGA1(NAI) 361 * + (T3AM(NDL,NEM,NFN)-T3AM(NDL,NEN,NFM)) 362 * * C2AM(NEM,NDL)*XIAJB(NFN,NAI) 363 ENDDO 364 ENDDO 365 ENDDO 366 ENDDO 367 ENDDO 368 ENDDO 369 ENDDO 370 ENDDO 371C 372 RETURN 373 END 374*---------------------------------------------------------------------* 375C /* Deck ccsdt_l03am */ 376*=====================================================================* 377 SUBROUTINE CCSDT_L03AM(L3AM,XINT1T,XINT2T,XIAJB,FOCK, 378 * L1AM,L2AM,SCR1,FOCKD,FCKFLD,L3AM2) 379*---------------------------------------------------------------------* 380* 381* Purpose: compute zero-order triples lagrangian multipliers 382* 383*---------------------------------------------------------------------* 384 IMPLICIT NONE 385#include "priunit.h" 386#include "ccsdinp.h" 387#include "ccsdsym.h" 388#include "ccfield.h" 389#include "maxorb.h" 390#include "ccorb.h" 391 392#if defined (SYS_CRAY) 393 REAL XINT1T(*), XINT2T(*), XIAJB(*), FOCK(*) 394 REAL L3AM(*), L2AM(*), L1AM(*), L3AM2(*) 395 REAL SCR1(*), FOCKD(*), FCKFLD(*) 396#else 397 DOUBLE PRECISION XINT1T(*), XINT2T(*), XIAJB(*), FOCK(*) 398 DOUBLE PRECISION L3AM(*), L2AM(*), L1AM(*), L3AM2(*) 399 DOUBLE PRECISION SCR1(*), FOCKD(*), FCKFLD(*) 400#endif 401 402 CALL DZERO(L3AM,NT1AMX*NT1AMX*NT1AMX) 403 404 CALL CCSDT_L3AM_R(L3AM,0.0D0,XINT1T,XINT2T,XIAJB,FOCK, 405 * L1AM,L2AM,SCR1,FOCKD,.FALSE.) 406 407 CALL CCSDT_3AM(L3AM,0.0D0,SCR1,FOCKD,NONHF,FCKFLD,.TRUE.,L3AM2) 408 409 RETURN 410 END 411*---------------------------------------------------------------------* 412* END OF SUBROUTINE CCSDT_L03AM * 413*---------------------------------------------------------------------* 414*---------------------------------------------------------------------* 415C /* Deck ccsdt_l3am_r */ 416*=====================================================================* 417 SUBROUTINE CCSDT_L3AM_R(T3BAR,FREQ,XINT1,XINT2,XIAJB,FOCK,T1AM, 418 * T2AM,SCR1,FOCKD,DIVIDE) 419*---------------------------------------------------------------------* 420* 421* Purpose: compute contributions to left triples amplitudes 422* 423*---------------------------------------------------------------------* 424 IMPLICIT NONE 425#include "priunit.h" 426#include "maxorb.h" 427#include "ccorb.h" 428#include "ccsdinp.h" 429#include "ccsdsym.h" 430#include "ccfield.h" 431C 432#if defined (SYS_CRAY) 433 REAL XINT1(NT1AMX,NVIRT,NVIRT) 434 REAL XINT2(NT1AMX,NRHFT,NRHFT) 435 REAL XIAJB(NT1AMX,NT1AMX), FOCK(NORBT,NORBT) 436 REAL T3BAR(NT1AMX,NT1AMX,NT1AMX), SCR1(NT1AMX) 437 REAL T1AM(NT1AMX), T2AM(NT1AMX,NT1AMX), FOCKD(*) 438 REAL AIBJCK, FREQ, TWO 439#else 440 DOUBLE PRECISION XINT1(NT1AMX,NVIRT,NVIRT) 441 DOUBLE PRECISION XINT2(NT1AMX,NRHFT,NRHFT) 442 DOUBLE PRECISION XIAJB(NT1AMX,NT1AMX), FOCK(NORBT,NORBT) 443 DOUBLE PRECISION T3BAR(NT1AMX,NT1AMX,NT1AMX), SCR1(NT1AMX) 444 DOUBLE PRECISION T1AM(NT1AMX), T2AM(NT1AMX,NT1AMX), FOCKD(*) 445 DOUBLE PRECISION AIBJCK, FREQ, TWO 446#endif 447 448 PARAMETER (TWO = 2.0D0) 449 450 LOGICAL DIVIDE 451 452 INTEGER NI, NA, NAI, NK, NC, NCK, NJ, NCJ, NBJ, NBK, NBI, NCI, 453 & NAJ, NAK, ND, NDI, NDJ, NDK, NL, NAL, NBL, NB 454C 455 ! we can just divide by orbital energy denominators if we 456 ! have non-HF external fields 457 IF (DIVIDE .AND. NONHF) 458 & CALL QUIT('NONHF Problem in CCSDT_L3AM_R!') 459C 460 DO NI = 1,NRHFT 461 DO NA = 1,NVIRT 462 NAI = NVIRT*(NI-1) + NA 463 SCR1(NAI) = FOCKD(NRHFT+NA) - FOCKD(NI) 464 END DO 465 END DO 466C 467 DO 100 NK = 1,NRHFT 468 DO 105 NC = 1,NVIRT 469C 470 NCK = NVIRT*(NK-1) + NC 471C 472 DO 110 NJ = 1,NRHFT 473 NCJ = NVIRT*(NJ-1) + NC 474 DO 120 NB = 1,NVIRT 475C 476 NBJ = NVIRT*(NJ-1) + NB 477 NBK = NVIRT*(NK-1) + NB 478C 479 DO 130 NI = 1,NRHFT 480 NBI = NVIRT*(NI-1) + NB 481 NCI = NVIRT*(NI-1) + NC 482 DO 135 NA = 1,NVIRT 483C 484 NAI = NVIRT*(NI-1) + NA 485 NAJ = NVIRT*(NJ-1) + NA 486 NAK = NVIRT*(NK-1) + NA 487C 488 AIBJCK = 0.0D0 489C 490C T1* = TWO T1 => Factor of two 491 AIBJCK = AIBJCK + XIAJB(NBJ,NCK)*T1AM(NAI) 492C 493C T1* = TWO T1 => Factor of two 494 AIBJCK = AIBJCK - XIAJB(NBJ,NCI)*T1AM(NAK) 495C 496 AIBJCK = AIBJCK + FOCK(NK,NRHFT+NC)* 497 * T2AM(NAI,NBJ) 498C 499 AIBJCK = AIBJCK - FOCK(NJ,NRHFT+NC)* 500 * T2AM(NAI,NBK) 501C 502 DO 140 ND = 1,NVIRT 503C 504 NDI = NVIRT*(NI-1) + ND 505 NDJ = NVIRT*(NJ-1) + ND 506 NDK = NVIRT*(NK-1) + ND 507C 508 AIBJCK = AIBJCK + 509 * (TWO*XINT1(NCK,ND,NB)-XINT1(NBK,ND,NC))* 510 * T2AM(NDJ,NAI) 511C 512 AIBJCK = AIBJCK - 513 * XINT1(NBI,ND,NC)* 514 * T2AM(NDK,NAJ) 515C 516 140 CONTINUE 517C 518 DO 150 NL = 1,NRHFT 519C 520 NAL = NVIRT*(NL-1) + NA 521 NBL = NVIRT*(NL-1) + NB 522C 523 AIBJCK = AIBJCK - 524 * (TWO*XINT2(NCK,NJ,NL)-XINT2(NCJ,NK,NL))* 525 * T2AM(NBL,NAI) 526C 527 AIBJCK = AIBJCK + 528 * XINT2(NCJ,NI,NL)* 529 * T2AM(NBK,NAL) 530C 531 150 CONTINUE 532C 533 IF (DIVIDE) THEN 534 AIBJCK = AIBJCK/(SCR1(NAI)+SCR1(NBJ)+SCR1(NCK)-FREQ) 535 END IF 536C 537 T3BAR(NAI,NBJ,NCK) = T3BAR(NAI,NBJ,NCK) + AIBJCK 538 T3BAR(NAI,NCK,NBJ) = T3BAR(NAI,NCK,NBJ) + AIBJCK 539 T3BAR(NBJ,NAI,NCK) = T3BAR(NBJ,NAI,NCK) + AIBJCK 540 T3BAR(NCK,NAI,NBJ) = T3BAR(NCK,NAI,NBJ) + AIBJCK 541 T3BAR(NBJ,NCK,NAI) = T3BAR(NBJ,NCK,NAI) + AIBJCK 542 T3BAR(NCK,NBJ,NAI) = T3BAR(NCK,NBJ,NAI) + AIBJCK 543C 544 135 CONTINUE 545 130 CONTINUE 546 120 CONTINUE 547 110 CONTINUE 548 105 CONTINUE 549 100 CONTINUE 550C 551C------------------------------------------------------ 552C Get rid of amplitudes that are not allowed. 553C------------------------------------------------------ 554C 555 CALL CCSDT_CLEAN_T3(T3BAR,NT1AMX,NVIRT,NRHFT) 556C 557 RETURN 558 END 559*---------------------------------------------------------------------* 560* END OF SUBROUTINE CCSDT_L3AM_R * 561*---------------------------------------------------------------------* 562C /* Deck cc_lhpart_noddy */ 563*=====================================================================* 564 SUBROUTINE CC_LHPART_NODDY(OMEGA1,OMEGA2,L3AM,FREQ, 565 & FOCKD,FIELD, 566 & XOOOO,XOVVO,XOOVV,XVVVV, 567 & T2AM,XINT1S,XINT2S, 568 & WORK,LWORK) 569*---------------------------------------------------------------------* 570* 571* Purpose: solve 'left' triples equations and partition the 572* triples solution vector into an effective lhs vector 573* 574* Note: in case of non-HF external fields zero-order 575* cluster amplitudes must be available on file 'T3AMPL' 576* 577*---------------------------------------------------------------------* 578 IMPLICIT NONE 579#include "dummy.h" 580#include "ccsdinp.h" 581#include "ccsdsym.h" 582#include "ccfield.h" 583#include "maxorb.h" 584#include "ccorb.h" 585 586#if defined (SYS_CRAY) 587 REAL OMEGA1(*), OMEGA2(*), L3AM(*), FREQ, FOCKD(*) 588 REAL FIELD(*), T2AM(*), XINT1S(*), XINT2S(*), WORK(*) 589 REAL XOOOO(*), XOVVO(*), XOOVV(*), XVVVV(*) 590#else 591 DOUBLE PRECISION OMEGA1(*), OMEGA2(*), L3AM(*), FREQ, FOCKD(*) 592 DOUBLE PRECISION FIELD(*), T2AM(*), XINT1S(*), XINT2S(*), WORK(*) 593 DOUBLE PRECISION XOOOO(*), XOVVO(*), XOOVV(*), XVVVV(*) 594#endif 595 596 LOGICAL TRANSPOSE 597 INTEGER LWORK, KSCR1, KEND1, KL3SCR, LWRK1, KOME1, KOME2, KT03AM, 598 & LUTEMP, IJ, NIJ, INDEX 599 600 INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 601 602 CALL QENTER('CCLPARTN') 603 604*---------------------------------------------------------------------* 605* Solve triples equations: 606*---------------------------------------------------------------------* 607 KSCR1 = 1 608 KEND1 = KSCR1 + NT1AMX 609 IF (NONHF) THEN 610 KL3SCR = KEND1 611 KEND1 = KL3SCR + NT1AMX*NT1AMX*NT1AMX 612 END IF 613 614 LWRK1 = LWORK - KEND1 615 IF (LWRK1 .LT. 0) THEN 616 CALL QUIT('Insufficient space in CC_LHPART_NODDY (1)') 617 ENDIF 618 619 TRANSPOSE = .TRUE. 620 CALL CCSDT_3AM(L3AM,FREQ,WORK(KSCR1),FOCKD, 621 * NONHF,FIELD,TRANSPOSE,WORK(KL3SCR)) 622 623 CALL DSCAL(NT1AMX*NT1AMX*NT1AMX,-1.0d0,L3AM,1) 624 625*---------------------------------------------------------------------* 626* Add contribution to singles doubles result vectors: 627*---------------------------------------------------------------------* 628 KOME1 = 1 629 KOME2 = KOME1 + NT1AMX 630 KEND1 = KOME2 + NT1AMX*NT1AMX 631 632 IF (NONHF) THEN 633 KT03AM = KEND1 634 KEND1 = KT03AM + NT1AMX*NT1AMX*NT1AMX 635 END IF 636 637 LWRK1 = LWORK - KEND1 638 IF (LWRK1 .LT. 0) THEN 639 CALL QUIT('Insufficient space in CC_LHPART_NODDY (2)') 640 ENDIF 641 642 CALL DZERO(WORK(KOME1),NT1AMX) 643 CALL DZERO(WORK(KOME2),NT1AMX*NT1AMX) 644 645 CALL CC3_L3_OMEGA1_NODDY(WORK(KOME1),L3AM, 646 * XOOOO,XOVVO,XOOVV,XVVVV,T2AM) 647 648 CALL CC3_L3_OMEGA2_NODDY(WORK(KOME2),L3AM,XINT1S,XINT2S) 649 650c ------------------------------------------ 651c contributions from non-HF external fields: 652c ------------------------------------------ 653 IF (NONHF) THEN 654 LUTEMP = -1 655 CALL GPOPEN(LUTEMP,'T3AMPL','UNKNOWN',' ','UNFORMATTED', 656 & IDUMMY,.FALSE.) 657 REWIND LUTEMP 658 READ (LUTEMP) (WORK(KT03AM-1+I),I=1,NT1AMX*NT1AMX*NT1AMX) 659 CALL GPCLOSE(LUTEMP,'KEEP') 660 661 CALL CCSDT_E1AM(WORK(KOME1),L3AM,WORK(KT03AM),FIELD) 662 663 CALL CCSDT_E2AM(WORK(KOME2),L3AM,T2AM,FIELD) 664 ENDIF 665 666 667 DO I = 1,NT1AMX 668 OMEGA1(I) = OMEGA1(I) + WORK(KOME1+I-1) 669 END DO 670 671 DO I = 1,NT1AMX 672 DO J = 1,I 673 IJ = NT1AMX*(I-1) + J 674 NIJ = INDEX(I,J) 675 OMEGA2(NIJ) = OMEGA2(NIJ) + WORK(KOME2+IJ-1) 676 ENDDO 677 ENDDO 678 679 CALL QEXIT('CCLPARTN') 680 RETURN 681 END 682*---------------------------------------------------------------------* 683* END OF SUBROUTINE CC_LHPART_NODDY * 684*---------------------------------------------------------------------* 685C /* Deck cc3__l3_omega */ 686 SUBROUTINE CC3_L3_OMEGA2_NODDY(OMEGA2,L3AM,XINT1,XINT2) 687C 688C Kasper Hald, Spring 2002. 689C 690C Calculate the doubles of the L3 part of CC3 left hand side doubles 691C 692#include "implicit.h" 693#include "priunit.h" 694#include "inforb.h" 695#include "ccsdinp.h" 696#include "ccsdsym.h" 697C 698 INTEGER INDEX, NAI, NBJ, NBK, NAK, NCK, NCJ, NCI, NDJ, NDK, NDI 699 INTEGER NBL, NCL, NAL, NAJ, NBI, AI 700C 701#if defined (SYS_CRAY) 702 REAL OMEGA2(NT1AMX,NT1AMX) 703 REAL L3AM(NT1AMX,NT1AMX,NT1AMX) 704 REAL XINT1(NT1AMX,NVIRT,NVIRT) 705 REAL XINT2(NT1AMX,NRHFT,NRHFT) 706 REAL ONE, TWO, XAIBJ 707#else 708 DOUBLE PRECISION OMEGA2(NT1AMX,NT1AMX) 709 DOUBLE PRECISION L3AM(NT1AMX,NT1AMX,NT1AMX) 710 DOUBLE PRECISION XINT1(NT1AMX,NVIRT,NVIRT) 711 DOUBLE PRECISION XINT2(NT1AMX,NRHFT,NRHFT) 712 DOUBLE PRECISION ONE, TWO, XAIBJ 713#endif 714 PARAMETER (ONE = 1.0D0, TWO = 2.0D0) 715C 716C 717C INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 718C 719C=================================================== 720C Calculate the L3 contributions to omega2 721C=================================================== 722C 723 CALL DZERO(OMEGA2,NT1AMX*NT1AMX) 724C 725 DO NI = 1,NRHFT 726 DO NA = 1,NVIRT 727 NAI = NVIRT*(NI-1) + NA 728C 729 DO NJ = 1,NRHFT 730 DO NB = 1,NVIRT 731 NBJ = NVIRT*(NJ-1) + NB 732C 733 DO NK = 1,NRHFT 734 NBK = NVIRT*(NK-1) + NB 735 NAK = NVIRT*(NK-1) + NA 736 DO NC = 1,NVIRT 737C 738 NCK = NVIRT*(NK-1) + NC 739 NCJ = NVIRT*(NJ-1) + NC 740 NCI = NVIRT*(NI-1) + NC 741C 742 DO ND = 1,NVIRT 743C 744 NDJ = NVIRT*(NJ-1) + ND 745 NDK = NVIRT*(NK-1) + ND 746 NDI = NVIRT*(NI-1) + ND 747C 748 OMEGA2(NBJ,NAI) = OMEGA2(NBJ,NAI) + 749 * L3AM(NBJ,NCI,NDK)*XINT1(NDK,NC,NA) 750C 751C 752 ENDDO 753C 754 DO NL = 1,NRHFT 755C 756 NBL = NVIRT*(NL-1) + NB 757 NCL = NVIRT*(NL-1) + NC 758 NAL = NVIRT*(NL-1) + NA 759C 760 OMEGA2(NBJ,NAI) = OMEGA2(NBJ,NAI) 761 * - L3AM(NBJ,NAK,NCL)*XINT2(NCL,NI,NK) 762C 763 ENDDO 764C 765 ENDDO 766 ENDDO 767C 768 ENDDO 769 ENDDO 770C 771 ENDDO 772 ENDDO 773C 774C---------------------------------------------------- 775C Make P^{ab}_{ij} for the T3BAR contributions. 776C---------------------------------------------------- 777C 778 DO NAI = 1,NT1AMX 779 DO NBJ = 1,NAI 780C 781 XAIBJ = OMEGA2(NAI,NBJ) + OMEGA2(NBJ,NAI) 782 OMEGA2(NAI,NBJ) = XAIBJ 783 OMEGA2(NBJ,NAI) = XAIBJ 784C 785 ENDDO 786 ENDDO 787C 788 RETURN 789 END 790C /* Deck cc3__l3_omega1_noddy */ 791 SUBROUTINE CC3_L3_OMEGA1_NODDY(OMEGA1,L3AM, 792 * XOOOO,XOVVO,XOOVV,XVVVV,T2AM) 793C 794C Kasper Hald, Spring 2002. 795C 796C Calculate the doubles of the L3 part of CC3 left hand side singles 797C 798#include "implicit.h" 799#include "priunit.h" 800#include "inforb.h" 801#include "ccsdinp.h" 802#include "ccsdsym.h" 803 804 LOGICAL LOCDBG 805 PARAMETER (LOCDBG = .FALSE.) 806C 807 INTEGER INDEX, NAI, NBJ, NBK, NAK, NCK, NCJ, NCI, NDJ, NDK, NDI 808 INTEGER NBL, NCL, NAL, NAJ, NBI, AI 809C 810#if defined (SYS_CRAY) 811 REAL OMEGA1(NT1AMX) 812 REAL L3AM(NT1AMX,NT1AMX,NT1AMX) 813 REAL XOOOO(NRHFT,NRHFT,NRHFT,NRHFT) 814 REAL XOVVO(NRHFT,NVIRT,NVIRT,NRHFT) 815 REAL XOOVV(NRHFT,NRHFT,NVIRT,NVIRT) 816 REAL XVVVV(NVIRT,NVIRT,NVIRT,NVIRT) 817 REAL T2AM(NT1AMX,NT1AMX) 818 REAL ONE, TWO, XAIBJ 819#else 820 DOUBLE PRECISION OMEGA1(NT1AMX) 821 DOUBLE PRECISION L3AM(NT1AMX,NT1AMX,NT1AMX) 822 DOUBLE PRECISION XOOOO(NRHFT,NRHFT,NRHFT,NRHFT) 823 DOUBLE PRECISION XOVVO(NRHFT,NVIRT,NVIRT,NRHFT) 824 DOUBLE PRECISION XOOVV(NRHFT,NRHFT,NVIRT,NVIRT) 825 DOUBLE PRECISION XVVVV(NVIRT,NVIRT,NVIRT,NVIRT) 826 DOUBLE PRECISION T2AM(NT1AMX,NT1AMX) 827 DOUBLE PRECISION ONE, TWO, XAIBJ 828#endif 829 PARAMETER (ONE = 1.0D0, TWO = 2.0D0) 830C 831C 832C INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 833C 834C=================================================== 835C Calculate the L3BAR contributions to omega1 836C=================================================== 837C 838C------------------------------------------- 839C 4 occupied integrals in integrals 840C------------------------------------------- 841C 842 if (.true.) then 843 DO NA = 1, NVIRT 844 DO NB = 1, NVIRT 845 DO NC = 1, NVIRT 846 DO NI = 1, NRHFT 847 NAI = NVIRT*(NI-1) + NA 848 DO NJ = 1, NRHFT 849 NBJ = NVIRT*(NJ-1) + NB 850 DO NK = 1, NRHFT 851 NCK = NVIRT*(NK-1) + NC 852 DO NL = 1, NRHFT 853 NAL = NVIRT*(NL-1) + NA 854 DO NM = 1, NRHFT 855 NCM = NVIRT*(NM-1) + NC 856 OMEGA1(NAI) = OMEGA1(NAI) 857 * + L3AM(NBJ,NAL,NCM)*T2AM(NBJ,NCK)*XOOOO(NI,NL,NK,NM) 858 ENDDO 859 ENDDO 860 ENDDO 861 ENDDO 862 ENDDO 863 ENDDO 864 ENDDO 865 ENDDO 866 endif 867C 868C------------------------------------------- 869C 4 virtual integrals in integrals 870C------------------------------------------- 871C 872 if (.true.) then 873 DO NA = 1, NVIRT 874 DO NB = 1, NVIRT 875 DO NC = 1, NVIRT 876 DO ND = 1, NVIRT 877 DO NE = 1, NVIRT 878 DO NI = 1, NRHFT 879 NAI = NVIRT*(NI-1) + NA 880 NDI = NVIRT*(NI-1) + ND 881 DO NJ = 1, NRHFT 882 NBJ = NVIRT*(NJ-1) + NB 883 DO NK = 1, NRHFT 884 NEK = NVIRT*(NK-1) + NE 885 NCK = NVIRT*(NK-1) + NC 886C 887 OMEGA1(NAI) = OMEGA1(NAI) 888 * + L3AM(NBJ,NDI,NEK)*T2AM(NBJ,NCK)*XVVVV(ND,NA,NE,NC) 889 ENDDO 890 ENDDO 891 ENDDO 892 ENDDO 893 ENDDO 894 ENDDO 895 ENDDO 896 ENDDO 897 endif 898C 899C-------------------------------------------------------- 900C 2 terms with g_{ovvo} and 2 terms with g_{oovv} 901C-------------------------------------------------------- 902C 903 if (.true.) then 904 DO NA = 1, NVIRT 905 DO NB = 1, NVIRT 906 DO NC = 1, NVIRT 907 DO ND = 1, NVIRT 908 DO NI = 1, NRHFT 909 NAI = NVIRT*(NI-1) + NA 910 NCI = NVIRT*(NI-1) + NC 911 NDI = NVIRT*(NI-1) + ND 912 DO NJ = 1, NRHFT 913 NBJ = NVIRT*(NJ-1) + NB 914 DO NK = 1, NRHFT 915 NAK = NVIRT*(NK-1) + NA 916 NCK = NVIRT*(NK-1) + NC 917 NDK = NVIRT*(NK-1) + ND 918 DO NL = 1, NRHFT 919 NAL = NVIRT*(NL-1) + NA 920 NCL = NVIRT*(NL-1) + NC 921 NDL = NVIRT*(NL-1) + ND 922C 923C - L3^{daf}_{lmn} * t^{de}_{lm} * g_{iefn} 924 OMEGA1(NAI) = OMEGA1(NAI) 925 * - L3AM(NBJ,NAK,NDL)*T2AM(NBJ,NCK)*XOVVO(NI,NC,ND,NL) 926C 927C - L3^{daf}_{lnm} * t^{de}_{lm} * g_{infe} 928 OMEGA1(NAI) = OMEGA1(NAI) 929 * - L3AM(NBJ,NAL,NDK)*T2AM(NBJ,NCK)*XOOVV(NI,NL,ND,NC) 930C 931C - L3^{def}_{lin} * t^{de}_{lm} * g_{mafn} 932 OMEGA1(NAI) = OMEGA1(NAI) 933 * - L3AM(NBJ,NCI,NDL)*T2AM(NBJ,NCK)*XOVVO(NK,NA,ND,NL) 934C 935C - L3^{def}_{lni} * t^{de}_{lm} * g_{mnfa} 936 OMEGA1(NAI) = OMEGA1(NAI) 937 * - L3AM(NBJ,NCL,NDI)*T2AM(NBJ,NCK)*XOOVV(NK,NL,ND,NA) 938C 939 ENDDO 940 ENDDO 941 ENDDO 942 ENDDO 943 ENDDO 944 ENDDO 945 ENDDO 946 ENDDO 947 endif 948C 949C===================================== 950C Write out all results 951C===================================== 952C 953 if (locdbg) then 954 write(lupri,*) ' ' 955 do nai = 1, nt1amx 956 if (abs(omega1(nai)) .gt. 1.0D-12) then 957 write(lupri,*) 'Omega1_noddy(',nai,') = ', 958 * omega1(nai) 959 endif 960 enddo 961 write(lupri,*) ' ' 962 end if 963C 964 RETURN 965 END 966C /* Deck ccfop_tran1 */ 967 SUBROUTINE CCFOP_TRAN1(XINT1,XINT2,XINT3,XINT4, 968 * XLAMDP,XLAMDH,AOINT,IDEL) 969C 970C 971C 972#include "implicit.h" 973#include "priunit.h" 974#include "inforb.h" 975#include "ccsdinp.h" 976#include "ccsdsym.h" 977 DIMENSION XINT1(NRHFT,NVIRT,NVIRT,NRHFT) 978 DIMENSION XINT2(NRHFT,NRHFT,NVIRT,NVIRT) 979 DIMENSION XINT3(NRHFT,NRHFT,NRHFT,NRHFT) 980 DIMENSION XINT4(NVIRT,NVIRT,NVIRT,NVIRT) 981 DIMENSION XLAMDP(NBAST,NORBT), XLAMDH(NBAST,NORBT) 982 DIMENSION AOINT(NNBAST,NBAST) 983C 984 LOGICAL LDEBUG 985C 986 INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 987C 988 LDEBUG = .TRUE. 989C 990C---------------------------------------- 991C Calculate integrals : 992C---------------------------------------- 993C 994 DO 100 G = 1,NBAST 995 DO 110 IB = 1,NBAST 996 DO 120 A = 1,NBAST 997 NAB = INDEX(A,IB) 998C 999 if (aoint(nab,g) .eq. 0.0d0) goto 120 1000C 1001 DO NC = 1,NVIRT 1002 DO ND = 1,NVIRT 1003 DO NK = 1,NRHFT 1004 DO NL = 1,NRHFT 1005C 1006 XINT1(NK,NC,ND,NL) = XINT1(NK,NC,ND,NL) 1007 * + AOINT(NAB,G)*XLAMDP(A,NK)*XLAMDH(IB,NRHFT+NC) 1008 * *XLAMDP(G,NRHFT+ND)*XLAMDH(IDEL,NL) 1009C 1010 ENDDO 1011 ENDDO 1012 ENDDO 1013 ENDDO 1014C 1015 DO NC = 1,NVIRT 1016 DO ND = 1,NVIRT 1017 DO NK = 1,NRHFT 1018 DO NL = 1,NRHFT 1019C 1020 XINT2(NK,NL,NC,ND) = XINT2(NK,NL,NC,ND) 1021 * + AOINT(NAB,G)*XLAMDP(A,NK)*XLAMDH(IB,NL) 1022 * *XLAMDP(G,NRHFT+NC)*XLAMDH(IDEL,NRHFT+ND) 1023C 1024 ENDDO 1025 ENDDO 1026 ENDDO 1027 ENDDO 1028C 1029 DO NK = 1,NRHFT 1030 DO NL = 1,NRHFT 1031 DO NM = 1,NRHFT 1032 DO NN = 1,NRHFT 1033C 1034 XINT3(NK,NL,NM,NN) = XINT3(NK,NL,NM,NN) 1035 * + AOINT(NAB,G)*XLAMDP(A,NK)*XLAMDH(IB,NL) 1036 * *XLAMDP(G,NM)*XLAMDH(IDEL,NN) 1037C 1038 ENDDO 1039 ENDDO 1040 ENDDO 1041 ENDDO 1042C 1043 DO NC = 1,NVIRT 1044 DO ND = 1,NVIRT 1045 DO NE = 1,NVIRT 1046 DO NF = 1,NVIRT 1047C 1048 XINT4(NC,ND,NE,NF) = XINT4(NC,ND,NE,NF) 1049 * + AOINT(NAB,G)*XLAMDP(A,NRHFT+NC)*XLAMDH(IB,NRHFT+ND) 1050 * *XLAMDP(G,NRHFT+NE)*XLAMDH(IDEL,NRHFT+NF) 1051C 1052 ENDDO 1053 ENDDO 1054 ENDDO 1055 ENDDO 1056C 1057 120 CONTINUE 1058 110 CONTINUE 1059 100 CONTINUE 1060C 1061 RETURN 1062 END 1063C /* Deck cc_lhtr_noddy2 */ 1064 SUBROUTINE CC_LHTR_NODDY2(OMEGA1,OMEGA2,OMEGA3,T1AM,T2AM,T3AM, 1065 * C1AM,C2AM,C3AM,XLAMDP,XLAMDH,FOCK, 1066 * WORK,LWORK) 1067C 1068C Written by Henrik Koch 19-Sep-1994 1069C 1070C Version 1.0 1071C 1072C Purpose: 1073C 1074C Calculate the iterative triples corrections to the CCSD 1075C equations. 1076C 1077C NB! The T2 amplitudes are assumed to be a full square. 1078C 1079C 1080C NB! It is assumed that the vectors are allocated in the following 1081C order: 1082C T1AM(*), OMEGA1(*), OMEGA2(*), T2AM(*), SCR(*), WRK(*). 1083C 1084C IOPT = 0 No storage of T3, L3 1085C IOPT = 1 Store T3 and L3 in 'T3AMPL' and 'L3AMPL' 1086C 1087C 1088#include "implicit.h" 1089#include "priunit.h" 1090#include "dummy.h" 1091#include "iratdef.h" 1092#include "maxorb.h" 1093#include "maxash.h" 1094#include "aovec.h" 1095 PARAMETER (ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, TWO = 2.0D0) 1096C 1097 DIMENSION INDEXA(MXCORB_CC) 1098 DIMENSION OMEGA1(*),OMEGA2(*), OMEGA3(*) 1099 DIMENSION T1AM(*),T2AM(*), T3AM(*) 1100 DIMENSION C1AM(*), C2AM(*), C3AM(*) 1101 DIMENSION XLAMDP(*),XLAMDH(*), FOCK(*), WORK(LWORK) 1102 LOGICAL L2INCL 1103C 1104#include "ccfield.h" 1105#include "inforb.h" 1106CCN #include "infind.h" not compatible with R12! 1107#include "ccisao.h" 1108#include "r12int.h" 1109#include "blocks.h" 1110#include "inftap.h" 1111#include "ccsdinp.h" 1112#include "ccsdsym.h" 1113#include "ccsdio.h" 1114COMMENT COMMENT 1115COMMENT COMMENT 1116#include "mxcent.h" 1117#include "eritap.h" 1118COMMENT COMMENT 1119COMMENT COMMENT 1120C 1121C INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 1122C 1123 IF (NSYM .NE. 1) THEN 1124 CALL QUIT('CC_LHTR_NODDY2: No symmetry in this part yet') 1125 ENDIF 1126C 1127C------------------------ 1128C Dynamic Allocation. 1129C------------------------ 1130C 1131 KCMO = 1 1132 KFOCKD = KCMO + NLAMDT 1133 KINT1 = KFOCKD + NORBT 1134 KINT2 = KINT1 + NT1AMX*NVIRT*NVIRT 1135 KINT1T = KINT2 + NRHFT*NRHFT*NT1AMX 1136 KINT2T = KINT1T + NT1AMX*NVIRT*NVIRT 1137 KINT1S = KINT2T + NRHFT*NRHFT*NT1AMX 1138 KINT2S = KINT1S + NT1AMX*NVIRT*NVIRT 1139 KINT3S = KINT2S + NRHFT*NRHFT*NT1AMX 1140 KINT4S = KINT3S + NT1AMX*NVIRT*NVIRT 1141 KOME1 = KINT4S + NRHFT*NRHFT*NT1AMX 1142 KOME2 = KOME1 + NT1AMX 1143 KOME3 = KOME2 + NT1AMX*NT1AMX 1144 KIAJB = KOME3 + NT1AMX*NT1AMX*NT1AMX 1145 KYIAJB = KIAJB + NT1AMX*NT1AMX 1146 KEND1 = KYIAJB + NT1AMX*NT1AMX 1147 LWRK1 = LWORK - KEND1 1148C 1149C New Integrals. 1150 KIOVVO = KEND1 1151 KIOOVV = KIOVVO + NRHFT*NVIRT*NVIRT*NRHFT 1152 KIOOOO = KIOOVV + NRHFT*NRHFT*NVIRT*NVIRT 1153 KIVVVV = KIOOOO + NRHFT*NRHFT*NRHFT*NRHFT 1154 KEND1 = KIVVVV + NVIRT*NVIRT*NVIRT*NVIRT 1155 LWRK1 = LWORK - KEND1 1156C 1157 IF ((NONHF) .AND. (NFIELD .GT. 0)) THEN 1158 KONEP = KEND1 1159 KT3SCR = KONEP + NORBT*NORBT 1160 KEND1 = KT3SCR + NT1AMX*NT1AMX*NT1AMX 1161 LWRK1 = LWORK - KEND1 1162 ENDIF 1163C 1164C 1165 IF (LWRK1 .LT. 0) THEN 1166 CALL QUIT('Insufficient space in CCSD_TRIPLE') 1167 ENDIF 1168C 1169C-------------------------------- 1170C Initialize integral arrays. 1171C-------------------------------- 1172C 1173 CALL DZERO(WORK(KINT1),NT1AMX*NVIRT*NVIRT) 1174 CALL DZERO(WORK(KINT2),NT1AMX*NRHFT*NRHFT) 1175 CALL DZERO(WORK(KINT1T),NT1AMX*NVIRT*NVIRT) 1176 CALL DZERO(WORK(KINT2T),NT1AMX*NRHFT*NRHFT) 1177 CALL DZERO(WORK(KINT1S),NT1AMX*NVIRT*NVIRT) 1178 CALL DZERO(WORK(KINT2S),NT1AMX*NRHFT*NRHFT) 1179 CALL DZERO(WORK(KINT3S),NT1AMX*NVIRT*NVIRT) 1180 CALL DZERO(WORK(KINT4S),NT1AMX*NRHFT*NRHFT) 1181 CALL DZERO(WORK(KOME1),NT1AMX) 1182 CALL DZERO(WORK(KOME2),NT1AMX*NT1AMX) 1183 CALL DZERO(WORK(KOME3),NT1AMX*NT1AMX*NT1AMX) 1184 CALL DZERO(WORK(KIAJB),NT1AMX*NT1AMX) 1185 CALL DZERO(WORK(KYIAJB),NT1AMX*NT1AMX) 1186C 1187 CALL DZERO(WORK(KIOVVO),NRHFT*NVIRT*NVIRT*NRHFT) 1188 CALL DZERO(WORK(KIOOVV),NRHFT*NVIRT*NVIRT*NRHFT) 1189 CALL DZERO(WORK(KIOOOO),NRHFT*NRHFT*NRHFT*NRHFT) 1190 CALL DZERO(WORK(KIVVVV),NVIRT*NVIRT*NVIRT*NVIRT) 1191C 1192C======================================================= 1193C Get the one electron integrals from the fields 1194C======================================================= 1195C 1196 IF ((NONHF) .AND. NFIELD .GT. 0) THEN 1197 CALL DZERO(WORK(KONEP),N2BST(ISYMOP)) 1198 DO I = 1, NFIELD 1199 FF = EFIELD(I) 1200 CALL CC_ONEP(WORK(KONEP),WORK(KEND1),LWRK1,FF,1,LFIELD(I)) 1201 ENDDO 1202 CALL CC_FCKMO(WORK(KONEP),XLAMDP,XLAMDH, 1203 * WORK(KEND1),LWRK1,1,1,1) 1204 ENDIF 1205C 1206C==================================================== 1207C Start the loop over distributions of integrals. 1208C==================================================== 1209C 1210 IF (DIRECT) THEN 1211 CALL QUIT('Direct not implemented') 1212#ifdef NOT_IMPLEMENTED 1213 NTOSYM = 1 1214 DTIME = SECOND() 1215 CALL HERDI1(WORK(KEND1),LWRK1,IPRINT) 1216 DTIME = SECOND() - DTIME 1217 TIMHER1 = TIMHER1 + DTIME 1218#endif 1219 ELSE 1220 NTOSYM = NSYM 1221 ENDIF 1222C 1223 DO 100 ISYMD1 = 1,NTOSYM 1224C 1225 IF (DIRECT) THEN 1226 NTOT = MAXSHL 1227 ELSE 1228 NTOT = NBAS(ISYMD1) 1229 ENDIF 1230C 1231 DO 110 ILLL = 1,NTOT 1232C 1233C----------------------------------------------------------------- 1234C If direct calculate the integrals and transposed t2am. 1235C----------------------------------------------------------------- 1236C 1237 IF (DIRECT) THEN 1238 CALL QUIT('Direct not implemented') 1239#ifdef NOT_IMPLEMENTED 1240 DTIME = SECOND() 1241 CALL HERDI2(WORK(KEND1),LWRK1,INDEXA,ILLL,NUMDIS,IPRINT) 1242 DTIME = SECOND() - DTIME 1243 TIMHER2 = TIMHER2 + DTIME 1244COMMENT COMMENT 1245 KRECNR = KEND1 1246 KEND1 = KRECNR + (NBUFX(0) - 1)/IRAT + 1 1247 LWRK1 = LWORK - KEND1 1248 IF (LWRK1 .LT. 0) THEN 1249 CALL QUIT('Insufficient core in CC_FOP3') 1250 END IF 1251COMMENT COMMENT 1252#endif 1253 ELSE 1254 NUMDIS = 1 1255 ENDIF 1256C 1257C----------------------------------------------------- 1258C Loop over number of distributions in disk. 1259C----------------------------------------------------- 1260C 1261 DO 120 IDEL2 = 1,NUMDIS 1262C 1263 IF (DIRECT) THEN 1264#ifdef NOT_IMPLEMENTED 1265 IDEL = INDEXA(IDEL2) 1266 IF (NOAUXB) THEN 1267 IDUM = 1 1268 CALL IJKAUX(IDEL,IDUM,IDUM,IDUM) 1269 END IF 1270 ISYMD = ISAO(IDEL) 1271#endif 1272 ELSE 1273 IDEL = IBAS(ISYMD1) + ILLL 1274 ISYMD = ISYMD1 1275 ENDIF 1276C 1277 ISYDIS = MULD2H(ISYMD,ISYMOP) 1278C 1279C------------------------------------------ 1280C Work space allocation no. 2. 1281C------------------------------------------ 1282C 1283 KXINT = KEND1 1284 KEND2 = KXINT + NDISAO(ISYDIS) 1285 LWRK2 = LWORK - KEND2 1286C 1287 IF (LWRK2 .LT. 0) THEN 1288 WRITE(LUPRI,*) 'Need : ',KEND2,'Available : ',LWORK 1289 CALL QUIT('Insufficient space in CCSD_TRIPLE') 1290 ENDIF 1291C 1292C 1293C----------------------------------------- 1294C Read in batch of integrals. 1295C----------------------------------------- 1296C 1297c DTIME = SECOND() 1298 call quit('CC_LHTR_NODDY2: KRECNR not defined, FIXME') 1299 krecnr = 1 ! hjaaj Aug 07: to silence ftnchek 1300 CALL CCRDAO(WORK(KXINT),IDEL,IDEL2,WORK(KEND2),LWRK2, 1301 * WORK(KRECNR),DIRECT) 1302c DTIME = SECOND() - DTIME 1303c TIMRDAO = TIMRDAO + DTIME 1304C 1305C---------------------------------------------------- 1306C Calculate integrals needed in CCSD[T]. 1307C---------------------------------------------------- 1308C 1309 CALL CCSDT_TRAN1(WORK(KINT1T),WORK(KINT2T),XLAMDP, 1310 * XLAMDH,WORK(KXINT),IDEL) 1311C 1312 CALL CC3_TRAN2(WORK(KIAJB),WORK(KYIAJB),XLAMDP, 1313 * XLAMDH,WORK(KXINT),IDEL) 1314C 1315 CALL CCSDT_TRAN3(WORK(KINT1S),WORK(KINT2S),XLAMDP, 1316 * XLAMDH,WORK(KXINT),IDEL) 1317C 1318 CALL CCFOP_TRAN1(WORK(KIOVVO),WORK(KIOOVV), 1319 * WORK(KIOOOO),WORK(KIVVVV), 1320 * XLAMDP,XLAMDH,WORK(KXINT),IDEL) 1321 120 CONTINUE 1322 110 CONTINUE 1323 100 CONTINUE 1324C 1325C------------------------------------ 1326C Calculate the corrections from T3. 1327C------------------------------------ 1328C 1329 CALL T3_LEFT1(WORK(KOME1),WORK(KYIAJB),WORK(KIAJB), 1330 * C2AM,T3AM) 1331C 1332 CALL T3_LEFT2(WORK(KOME1),WORK(KYIAJB),WORK(KIAJB), 1333 * C2AM,T3AM) 1334C 1335 CALL T3_LEFT3(WORK(KOME1),WORK(KIAJB),C2AM,T3AM) 1336C 1337 DO I = 1,NT1AMX 1338 OMEGA1(I) = OMEGA1(I) + WORK(KOME1+I-1) 1339 ENDDO 1340C 1341C---------------------------------------------------------- 1342C Calculate the contributions to Omega2 from C3 1343C The result is written out inside the routine. 1344C---------------------------------------------------------- 1345C 1346 CALL DZERO(WORK(KOME1),NT1AMX) 1347 CALL DZERO(WORK(KOME2),NT1AMX*NT1AMX) 1348 CALL DZERO(WORK(KOME3),NT1AMX*NT1AMX*NT1AMX) 1349 1350 CALL CC3_L3_OMEGA1_NODDY(WORK(KOME1),C3AM, 1351 * WORK(KIOOOO),WORK(KIOVVO), 1352 * WORK(KIOOVV),WORK(KIVVVV),T2AM) 1353 1354 CALL CC3_L3_OMEGA2_NODDY(WORK(KOME2),C3AM, 1355 * WORK(KINT1S),WORK(KINT2S)) 1356 1357 IF (NONHF .AND. NFIELD.GT.0) THEN 1358 CALL CCSDT_E1AM(WORK(KOME1),C3AM,T3AM,WORK(KONEP)) 1359 1360 CALL CCSDT_E2AM(WORK(KOME2),C3AM,T2AM,WORK(KONEP)) 1361 ENDIF 1362 1363 DO 300 I = 1,NT1AMX 1364 OMEGA1(I) = OMEGA1(I) + WORK(KOME1+I-1) 1365 300 CONTINUE 1366 1367 CALL DAXPY(NT1AMX*NT1AMX,1.0D0,WORK(KOME2),1,OMEGA2,1) 1368C 1369C--------------------------------------------------------- 1370C Calculate the contributions to omega3 from C1, C2 1371C--------------------------------------------------------- 1372 DO I = 1, NRHFT 1373 WORK(KFOCKD-1+I) = 0.0d0 1374 ENDDO 1375 DO I = 1, NVIRT 1376 WORK(KFOCKD-1+NRHFT+I) = 1.0d0/3.0d0 1377 ENDDO 1378C 1379C for non-HF external fields the above trick for the 1380C orbital energies does solve all problems... 1381C --> change to CCSDT_L3AM_R 1382 IF (NONHF) CALL QUIT('fix me!') 1383C 1384C hjaaj-aug-2007: changed undefined KSCR1 to KEND1, hope it is OK. 1385 CALL CCSDT_L03AM(WORK(KOME3),WORK(KINT1T),WORK(KINT2T), 1386 * WORK(KIAJB),FOCK,C1AM,C2AM, 1387 * WORK(KEND1),WORK(KFOCKD), 1388 * WORK(KONEP),WORK(KT3SCR)) 1389 1390 IF (NONHF .AND. NFIELD.GT.0) THEN 1391 L2INCL = .TRUE. 1392 CALL CCSDT_E3AM(WORK(KOME3),C2AM,C3AM,WORK(KONEP),L2INCL) 1393 END IF 1394 1395 CALL DAXPY(NT1AMX*NT1AMX*NT1AMX,1.0D0,WORK(KOME3),1,OMEGA3,1) 1396C 1397C---------------------------------------- 1398C Print norm of result vectors. 1399C---------------------------------------- 1400C 1401 XNORM = DDOT(NT1AMX,OMEGA1,1,OMEGA1,1) 1402 write(lupri,*) 'Norm Omega1 at the end of noddy :',xnorm 1403 XNORM = DDOT(NT1AMX**2,OMEGA2,1,OMEGA2,1) 1404 write(lupri,*) 'Norm Omega2 at the end of noddy :',xnorm 1405 XNORM = DDOT(NT1AMX**3,OMEGA3,1,OMEGA3,1) 1406 write(lupri,*) 'Norm Omega3 at the end of noddy :',xnorm 1407C 1408 RETURN 1409 END 1410