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 19C /* Deck ccsd_triple */ 20 SUBROUTINE CC_FOPTRIPLES(OMEGA1,DENAB,DENIJ,T1AM,T2AM, 21 * XLAMDP,XLAMDH,WORK,LWORK) 22C 23C Written by K. Hald, Fall 2001 based on 24C CCSD_TRIPLE() written by Henrik Koch 19-Sep-1994 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 39#include "implicit.h" 40#include "priunit.h" 41#include "dummy.h" 42#include "maxorb.h" 43#include "maxash.h" 44#include "aovec.h" 45 PARAMETER (ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, TWO = 2.0D0) 46C 47 DIMENSION INDEXA(MXCORB_CC) 48 DIMENSION OMEGA1(*),DENAB(*),DENIJ(*),T1AM(*),T2AM(*) 49 DIMENSION XLAMDP(*),XLAMDH(*),WORK(LWORK) 50C 51COMMENT COMMENT COMMENT 52C#include "inforb.h" 53#include "ccorb.h" 54COMMENT COMMENT COMMENT 55CCN #include "infind.h" not compatible with R12! 56#include "ccisao.h" 57#include "r12int.h" 58#include "blocks.h" 59#include "inftap.h" 60#include "ccsdinp.h" 61#include "ccsdsym.h" 62#include "ccsdio.h" 63COMMENT COMMENT 64COMMENT COMMENT 65#include "mxcent.h" 66#include "eritap.h" 67#include "ccfield.h" 68#include "ccfop.h" 69COMMENT COMMENT 70COMMENT COMMENT 71C 72 LOGICAL ITERATIVE 73C 74 INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 75C 76 IF (NSYM .NE. 1) THEN 77 CALL QUIT('No symmetry in this part yet') 78 ENDIF 79C 80COMMENT COMMENT 81COMMENT COMMENT 82 write(lupri,*) 'Entered FOPTRIPLES' 83 write(lupri,*) 'NLAMDT = ',nlamdt 84 write(lupri,*) 'NLAMDS = ',nlamds 85 write(lupri,*) 'NORBT = ',norbt 86 write(lupri,*) 'NORBTS = ',norbts 87 call flshfo(lupri) 88C------------------------ 89C Dynamic Allocation. 90C------------------------ 91C 92 KCMO = 1 93 KOME1 = KCMO + NLAMDS 94COMMENT COMMENT 95COMMENT COMMENT CHANGED NLAMDT TO NLAMDS 96COMMENT COMMENT 97 KOME2 = KOME1 + NT1AMX 98 KOME11 = KOME2 + NT1AMX*NT1AMX 99 KOME22 = KOME11 + NT1AMX 100 KSCR1 = KOME22 + NT1AMX*NT1AMX 101 KFOCK = KSCR1 + NT1AMX 102 KFOCKD = KFOCK + N2BST(ISYMOP) 103 KINT1S = KFOCKD + NORBTS 104COMMENT COMMENT 105COMMENT COMMENT CHANGED NORBT TO NORBTS 106COMMENT COMMENT 107 KINT2S = KINT1S + NT1AMX*NVIRT*NVIRT 108 KINT3S = KINT2S + NRHFT*NRHFT*NT1AMX 109 KINT4S = KINT3S + NT1AMX*NVIRT*NVIRT 110 KINT1T = KINT4S + NRHFT*NRHFT*NT1AMX 111 KINT2T = KINT1T + NT1AMX*NVIRT*NVIRT 112 KIAJB = KINT2T + NRHFT*NRHFT*NT1AMX 113 KYIAJB = KIAJB + NT1AMX*NT1AMX 114 KT3AM = KYIAJB + NT1AMX*NT1AMX 115 KT3BAR = KT3AM + NT1AMX*NT1AMX*NT1AMX 116 KVIRT = KT3BAR + NT1AMX*NT1AMX*NT1AMX 117 KOCC = KVIRT + NVIRT*NVIRT*NVIRT*NRHFT 118 KKAPAA = KOCC + NRHFT*NRHFT*NRHFT*NVIRT 119 KKAPII = KKAPAA + NVIRT 120 KEND1 = KKAPII + NRHFT 121 LWRK1 = LWORK - KEND1 122C 123 IF (LWRK1 .LT. 0) THEN 124 CALL QUIT('Insufficient space in CCSD_TRIPLE') 125 ENDIF 126C 127C-------------------------------- 128C Initialize integral arrays. 129C-------------------------------- 130C 131 CALL DZERO(WORK(KINT1T),NT1AMX*NVIRT*NVIRT) 132 CALL DZERO(WORK(KINT2T),NT1AMX*NRHFT*NRHFT) 133 CALL DZERO(WORK(KINT1S),NT1AMX*NVIRT*NVIRT) 134 CALL DZERO(WORK(KINT2S),NT1AMX*NRHFT*NRHFT) 135 CALL DZERO(WORK(KINT3S),NT1AMX*NVIRT*NVIRT) 136 CALL DZERO(WORK(KINT4S),NT1AMX*NRHFT*NRHFT) 137 CALL DZERO(WORK(KIAJB),NT1AMX*NT1AMX) 138 CALL DZERO(WORK(KYIAJB),NT1AMX*NT1AMX) 139 CALL DZERO(WORK(KT3AM),NT1AMX*NT1AMX*NT1AMX) 140 CALL DZERO(WORK(KT3BAR),NT1AMX*NT1AMX*NT1AMX) 141 CALL DZERO(WORK(KOME1),NT1AMX) 142 CALL DZERO(WORK(KOME2),NT1AMX*NT1AMX) 143 CALL DZERO(WORK(KOME11),NT1AMX) 144 CALL DZERO(WORK(KOME22),NT1AMX*NT1AMX) 145 CALL DZERO(WORK(KVIRT),NVIRT*NVIRT*NVIRT*NRHFT) 146 CALL DZERO(WORK(KOCC),NRHFT*NRHFT*NRHFT*NVIRT) 147C 148C---------------------------------------------- 149C Read MO-coefficients from interface file. 150C---------------------------------------------- 151C 152 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY, 153 & .FALSE.) 154 REWIND LUSIFC 155C 156 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 157C 158 READ (LUSIFC) 159 READ (LUSIFC) 160 READ (LUSIFC) (WORK(KCMO+I-1), I=1,NLAMDS) 161COMMENT COMMENT COMMENT 162COMMENT COMMENT COMMENT CHANGED NLAMDT TO NLAMDS 163COMMENT COMMENT COMMENT 164C 165C 166 CALL CMO_REORDER(WORK(KCMO),WORK(KEND1),LWRK1) 167C 168C------------------------------------- 169C Read canonical orbital energies. 170C------------------------------------- 171C 172 REWIND LUSIFC 173C 174 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 175 READ (LUSIFC) 176 READ (LUSIFC) (WORK(KFOCKD+I-1), I=1,NORBTS) 177COMMENT COMMENT COMMENT COMMENT 178COMMENT COMMENT COMMENT COMMENT CHANGED NORBT TO NORBTS 179COMMENT COMMENT COMMENT COMMENT 180C 181 CALL GPCLOSE(LUSIFC,'KEEP') 182C 183 IF (FROIMP .OR. FROEXP) THEN 184 CALL CCSD_DELFRO(WORK(KFOCKD),WORK(KEND1),LWRK1) 185 ENDIF 186C 187 CALL FOCK_REORDER(WORK(KFOCKD),WORK(KEND1),LWRK1) 188C 189C-------------------------------------------------- 190C Read in AO fock matrix and transform to MO 191C-------------------------------------------------- 192C 193 LUFCK = -1 194C The T1 transformed density is used to calculate this AO fock matrix. 195C CALL GPOPEN(LUFCK,'CC_FCKH','OLD',' ','UNFORMATTED', 196C * IDUMMY,.FALSE.) 197C The CMO transformed density is used to calculate this AO fock matrix. 198 CALL GPOPEN(LUFCK,'CC_FCKREF','OLD',' ','UNFORMATTED', 199 * IDUMMY,.FALSE.) 200C 201 REWIND(LUFCK) 202 READ(LUFCK)(WORK(KFOCK + I-1),I = 1,N2BST(ISYMOP)) 203 CALL GPCLOSE(LUFCK,'KEEP' ) 204C 205 IF (IPRINT .GT. 50) THEN 206 CALL AROUND( 'FOP3 : Usual Fock AO matrix' ) 207 CALL CC_PRFCKAO(WORK(KFOCK),ISYMOP) 208 ENDIF 209C 210 CALL CC_FCKMO(WORK(KFOCK),WORK(KCMO),WORK(KCMO), 211 * WORK(KEND1),LWRK1,1,1,1) 212C 213COMMENT 214C IF (IPRINT .GT. 50) THEN 215COMMENT 216 CALL AROUND('OUTPUT OF FOCK') 217 CALL OUTPUT(WORK(KFOCK),1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI) 218COMMENT 219C ENDIF 220COMMENT 221C 222C---------------------------------------------- 223C Calculate the one electron integrals 224C from the field. Add this contri 225C to the fock matrix as well. 226C---------------------------------------------- 227C 228 IF ((NONHF) .AND. NFIELD .GT. 0) THEN 229C 230 KONEP = KEND1 231 KEND1 = KONEP + N2BST(ISYMOP) 232 LWRK1 = LWORK - KEND1 233C 234 CALL DZERO(WORK(KONEP),N2BST(ISYMOP)) 235C 236 DO I = 1, NFIELD 237C 238 KONEP2 = KEND1 239 KEND2 = KONEP2 + N2BST(ISYMOP) 240 LWRK2 = LWORK - KEND2 241C 242 IF (LWRK2 .LT. N2BST(ISYMOP)) THEN 243 CALL QUIT('Exceeded workmemory in CC_FOP3 (Field)') 244 ENDIF 245C 246 CALL DZERO(WORK(KONEP2),N2BST(ISYMOP)) 247 FF = EFIELD(I) 248 CALL CC_ONEP(WORK(KONEP2),WORK(KEND2),LWRK2,FF,1,LFIELD(I)) 249C 250 CALL DAXPY(N2BST(ISYMOP),1.0D0,WORK(KONEP2),1,WORK(KONEP),1) 251 ENDDO 252C 253 IF (IPRINT .GT. 50) THEN 254 CALL AROUND( 'FOP3 : Usual AO matrix from field(s)' ) 255 CALL CC_PRFCKAO(WORK(KONEP),ISYMOP) 256 ENDIF 257C 258 CALL CC_FCKMO(WORK(KONEP),WORK(KCMO),WORK(KCMO), 259 * WORK(KEND1),LWRK1,1,1,1) 260C 261COMMENT 262C IF (IPRINT .GT. 50) THEN 263COMMENT 264 CALL AROUND('OUTPUT OF STRENGTH ONEP') 265 CALL OUTPUT(WORK(KONEP),1,NORBT,1,NORBT,NORBT,NORBT, 266 * 1,LUPRI) 267COMMENT 268C ENDIF 269COMMENT 270C 271 CALL DAXPY(N2BST(ISYMOP),1.0D0,WORK(KONEP),1,WORK(KFOCK),1) 272COMMENT 273C IF (IPRINT .GT. 50) THEN 274COMMENT 275 CALL AROUND('OUTPUT OF FOCK+STRENGTH ONEP') 276 CALL OUTPUT(WORK(KFOCK),1,NORBT,1,NORBT,NORBT,NORBT, 277 * 1,LUPRI) 278COMMENT 279C ENDIF 280COMMENT 281 ENDIF 282C 283C--------------------------------------------------------- 284C Calculate the one electron integrals 285C from the field (no strength). 286C--------------------------------------------------------- 287C 288 IF (((NONHF) .AND. (NFIELD .GT. 0)) .OR. (DIPMOM)) THEN 289C 290 KONEP2 = KEND1 291 KEND1 = KONEP2 + N2BST(ISYMOP) 292 LWRK1 = LWORK - KEND1 293C 294 CALL DZERO(WORK(KONEP2),N2BST(ISYMOP)) 295C 296 IF (DIPMOM) THEN 297 FF = 1.0D0 298 ISYONE = -1 299 CALL CC_ONEP(WORK(KONEP2),WORK(KEND1),LWRK1,FF, 300 * ISYONE,'ZDIPLEN ') 301 ELSE 302 DO I = 1, NFIELD 303C 304 KONEP3 = KEND1 305 KEND2 = KONEP3 + N2BST(ISYMOP) 306 LWRK2 = LWORK - KEND2 307C 308 IF (LWRK2 .LT. N2BST(ISYMOP)) THEN 309 CALL QUIT('Exceeded workmemory in CC_FOP3 (Field)') 310 ENDIF 311C 312 CALL DZERO(WORK(KONEP3),N2BST(ISYMOP)) 313 FF = 1.0D0 314 ISYONE = 1 315 CALL CC_ONEP(WORK(KONEP3),WORK(KEND2),LWRK2,FF, 316 * ISYONE,LFIELD(I)) 317C 318 CALL DAXPY(N2BST(ISYMOP),1.0D0,WORK(KONEP3),1, 319 * WORK(KONEP2),1) 320 ENDDO 321 ENDIF 322C 323 IF (IPRINT .GT. 50) THEN 324 CALL AROUND( 'FOP3 : Usual AO matrix from field(s)' ) 325 CALL CC_PRFCKAO(WORK(KONEP2),ISYMOP) 326 ENDIF 327C 328 CALL CC_FCKMO(WORK(KONEP2),WORK(KCMO),WORK(KCMO), 329 * WORK(KEND1),LWRK1,1,1,1) 330C 331COMMENT 332C IF (IPRINT .GT. 50) THEN 333COMMENT 334 CALL AROUND('OUTPUT OF ONEP WITHOUT STREGTH') 335 CALL OUTPUT(WORK(KONEP2),1,NORBT,1,NORBT,NORBT,NORBT, 336 * 1,LUPRI) 337COMMENT 338C ENDIF 339COMMENT 340C 341 ENDIF 342C 343C==================================================== 344C Start the loop over distributions of integrals. 345C==================================================== 346C 347 IF (DIRECT) THEN 348 NTOSYM = 1 349 DTIME = SECOND() 350c CALL HERDI1(WORK(KEND1),LWRK1,IPRINT) 351 CALL QUIT('Direct not implemented') 352 DTIME = SECOND() - DTIME 353 TIMHER1 = TIMHER1 + DTIME 354 ELSE 355 NTOSYM = NSYM 356 ENDIF 357C 358 DO 100 ISYMD1 = 1,NTOSYM 359C 360 IF (DIRECT) THEN 361 NTOT = MAXSHL 362 ELSE 363 NTOT = NBAS(ISYMD1) 364 ENDIF 365C 366 DO 110 ILLL = 1,NTOT 367C 368C----------------------------------------------------------------- 369C If direct calculate the integrals and transposed t2am. 370C----------------------------------------------------------------- 371C 372 IF (DIRECT) THEN 373 DTIME = SECOND() 374c CALL HERDI2(WORK(KEND1),LWRK1,INDEXA,ILLL,NUMDIS,IPRINT) 375 CALL QUIT('Direct not implemented') 376 DTIME = SECOND() - DTIME 377 TIMHER2 = TIMHER2 + DTIME 378COMMENT COMMENT 379 KRECNR = KEND1 380 KEND1 = KRECNR + (NBUFX(0) - 1)/IRAT + 1 381 LWRK1 = LWORK - KEND1 382 IF (LWRK1 .LT. 0) THEN 383 CALL QUIT('Insufficient core in CC_FOP3') 384 END IF 385COMMENT COMMENT 386 ELSE 387 NUMDIS = 1 388 ENDIF 389C 390C----------------------------------------------------- 391C Loop over number of distributions in disk. 392C----------------------------------------------------- 393C 394 DO 120 IDEL2 = 1,NUMDIS 395C 396 IF (DIRECT) THEN 397 IDEL = INDEXA(IDEL2) 398 IF (NOAUXB) THEN 399 IDUM = 1 400 CALL IJKAUX(IDEL,IDUM,IDUM,IDUM) 401 END IF 402 ISYMD = ISAO(IDEL) 403 ELSE 404 IDEL = IBAS(ISYMD1) + ILLL 405 ISYMD = ISYMD1 406 ENDIF 407C 408 ISYDIS = MULD2H(ISYMD,ISYMOP) 409C 410C------------------------------------------ 411C Work space allocation no. 2. 412C------------------------------------------ 413C 414 KXINT = KEND1 415 KEND2 = KXINT + NDISAO(ISYDIS) 416 LWRK2 = LWORK - KEND2 417C 418 IF (LWRK2 .LT. 0) THEN 419 WRITE(LUPRI,*) 'Need : ',KEND2,'Available : ',LWORK 420 CALL QUIT('Insufficient space in CCSD_TRIPLE') 421 ENDIF 422C 423C 424C----------------------------------------- 425C Read in batch of integrals. 426C----------------------------------------- 427C 428 DTIME = SECOND() 429 CALL CCRDAO(WORK(KXINT),IDEL,IDEL2,WORK(KEND2),LWRK2, 430 * WORK(KRECNR),DIRECT) 431 DTIME = SECOND() - DTIME 432 TIMRDAO = TIMRDAO + DTIME 433C 434C---------------------------------------------------- 435C Calculate integrals needed in CCSD[T]. 436C---------------------------------------------------- 437C 438C CALL CCSDT_TRAN1(WORK(KINT1),WORK(KINT2),WORK(KCMO), 439C * WORK(KCMO),WORK(KXINT),IDEL) 440C 441COMMENT COMMENT 442COMMENT COMMENT 443 CALL CCSDT_TRAN1(WORK(KINT1T),WORK(KINT2T),WORK(KCMO), 444 * WORK(KCMO),WORK(KXINT),IDEL) 445C CALL CCSDT_TRAN1(WORK(KINT1T),WORK(KINT2T),XLAMDP, 446C * XLAMDH,WORK(KXINT),IDEL) 447C 448 CALL CCSDT_TRAN2(WORK(KIAJB),WORK(KYIAJB),WORK(KCMO), 449 * WORK(KXINT),IDEL) 450C 451 CALL CCSDT_TRAN3(WORK(KINT1S),WORK(KINT2S),XLAMDP, 452 * XLAMDH,WORK(KXINT),IDEL) 453COMMENT COMMENT 454COMMENT COMMENT 455 CALL CCSDT_TRAN3PT(WORK(KINT3S),WORK(KINT4S),WORK(KCMO), 456 * WORK(KXINT),IDEL) 457C 458 120 CONTINUE 459 110 CONTINUE 460 100 CONTINUE 461C 462COMMENT COMMENT COMMENT 463COMMENT COMMENT COMMENT 464 call print_integrals(work(kint3s),work(kint4s)) 465COMMENT COMMENT COMMENT 466COMMENT COMMENT COMMENT 467C------------------------------------------------ 468C Calculate the triple cluster amplitudes. 469C------------------------------------------------ 470C 471 rho1n = ddot(nt1amx,t1am,1,t1am,1) 472 write(lupri,*) 'Norm of T1AM before T3AM : ',rho1n 473 rho1n = ddot(nt1amx*nt1amx,t2am,1,t2am,1) 474 write(lupri,*) 'Norm of T2AM before T3AM : ',rho1n 475 rho1n = ddot(nt1amx*nt1amx*nt1amx,work(kt3am),1,work(kt3am),1) 476 write(lupri,*) 'Norm of T3AM before T3AM : ',rho1n 477 rho1n = ddot(nt1amx*nvirt*nvirt,work(kint3s),1,work(kint3s),1) 478 write(lupri,*) 'Norm of KINT3S before T3AM : ',rho1n 479 rho1n = ddot(nt1amx*nrhft*nrhft,work(kint4s),1,work(kint4s),1) 480 write(lupri,*) 'Norm of KINT4S before T3AM : ',rho1n 481 rho1n = ddot(nt1amx,work(kscr1),1,work(kscr1),1) 482 write(lupri,*) 'Norm of KSCR1 before call : ',rho1n 483 write(lupri,*) 'Before CCSDT_T3AM' 484 call flshfo(lupri) 485C 486 IF ((NFIELD .GT. 0) .AND. (NONHF)) THEN 487 ITERATIVE = .TRUE. 488 CALL CCSDT_T3AM(WORK(KT3AM),WORK(KINT3S),WORK(KINT4S),T2AM, 489 * WORK(KSCR1),WORK(KFOCKD),WORK(KT3BAR),WORK(KONEP), 490 * ITERATIVE) 491 ELSE 492 ITERATIVE = .FALSE. 493 CALL CCSDT_T3AM(WORK(KT3AM),WORK(KINT3S),WORK(KINT4S),T2AM, 494 * WORK(KSCR1),WORK(KFOCKD),DUMMY,DUMMY,ITERATIVE) 495 ENDIF 496C 497 call temppr(work(kt3am),1) 498C 499 rho1n = ddot(nt1amx*nt1amx*nt1amx,work(kt3am),1,work(kt3am),1) 500 write(lupri,*) 'Norm of T3AM after T3AM : ',rho1n 501C 502C----------------------------------------------------- 503C For reference calculate the CCSD(T) energy. 504C----------------------------------------------------- 505C 506 CALL AROUND('START OF CCSD(T) ENERGY') 507C 508COMMENT COMMENT 509COMMENT COMMENT 510 rho1n = ddot(nt1amx*nt1amx,work(kiajb),1,work(kiajb),1) 511 write(lupri,*) 'Norm of IAJB integrals = ',rho1n 512 rho1n = ddot(nt1amx*nvirt*nvirt,work(kint1t),1,work(kint1t),1) 513 write(lupri,*) 'Norm of INT1T integrals = ',rho1n 514 rho1n = ddot(nt1amx*nrhft*nrhft,work(kint2t),1,work(kint2t),1) 515 write(lupri,*) 'Norm of INT2T integrals = ',rho1n 516 rho1n = ddot(n2bst(isymop),work(kfock),1,work(kfock),1) 517 write(lupri,*) 'Norm of FOCK iintermediate = ',rho1n 518 rho1n = ddot(nt1amx*nt1amx*nt1amx,work(kt3am),1,work(kt3am),1) 519 write(lupri,*) 'Norm of T3AM before calc. = ',rho1n 520 xnorm = ddot(nt1amx,work(kome1),1,work(kome1),1) 521 write(lupri,*) 'Omega1 norm before calc. = ',xnorm 522 xnorm = ddot(nt1amx*nt1amx,work(kome2),1,work(kome2),1) 523 write(lupri,*) 'Omega2 norm before calc. = ',xnorm 524COMMENT COMMENT 525COMMENT COMMENT 526 CALL CCSDT_OMEGA1OLD(WORK(KOME1),WORK(KIAJB),WORK(KT3AM)) 527C 528 CALL CCSDT_OMEGA2OLD(WORK(KOME2),WORK(KINT1T),WORK(KINT2T), 529 * WORK(KFOCK),WORK(KT3AM)) 530C 531COMMENT COMMENT 532COMMENT COMMENT 533 xnorm = ddot(nt1amx,work(kome1),1,work(kome1),1) 534 write(lupri,*) 'Omega1 norm after calc. = ',xnorm 535 xnorm = ddot(nt1amx*nt1amx,work(kome2),1,work(kome2),1) 536 write(lupri,*) 'Omega2 norm after calc. = ',xnorm 537 xnorm = ddot(nt1amx,t1am,1,t1am,1) 538 write(lupri,*) 'T1 norm = ',xnorm 539 xnorm = ddot(nt1amx*nt1amx,t2am,1,t2am,1) 540 write(lupri,*) 'T2 norm = ',xnorm 541COMMENT COMMENT 542COMMENT COMMENT 543 ENERT1 = ddot(nt1amx,t1am,1,work(kome1),1) 544 write(LUPRI,*) 'The E5(ST) contribution :',2.0D0*ENERT1 545C 546 CALL OUTPUT(WORK(KOME2),1,NT1AMX,1,NT1AMX,NT1AMX,NT1AMX,1,LUPRI) 547 CALL ECCSDT(T2AM,WORK(KOME2)) 548C 549 CALL AROUND('END OF CCSD(T) ENERGY') 550C 551C--------------------------------------------- 552C Calculate the triple bar amplitudes. 553C--------------------------------------------- 554C 555COMMENT COMMENT COMMENT 556COMMENT COMMENT COMMENT 557 rho1n = ddot(nt1amx*nt1amx*nt1amx,work(kt3bar),1,work(kt3bar),1) 558 write(lupri,*) 'Norm of T3BAR before call : ',rho1n 559 rho1n = ddot(nt1amx*nvirt*nvirt,work(kint3s),1,work(kint3s),1) 560 write(lupri,*) 'Norm of KINT3S before call : ',rho1n 561 rho1n = ddot(nt1amx*nrhft*nrhft,work(kint4s),1,work(kint4s),1) 562 write(lupri,*) 'Norm of KINT4S before call : ',rho1n 563 rho1n = ddot(nt1amx*nt1amx,work(kiajb),1,work(kiajb),1) 564 write(lupri,*) 'Norm of KIAJB before call : ',rho1n 565 rho1n = ddot(n2bst(isymop),work(kfock),1,work(kfock),1) 566 write(lupri,*) 'Norm of KFOCK before call : ',rho1n 567 rho1n = ddot(nt1amx,t1am,1,t1am,1) 568 write(lupri,*) 'Norm of T1AM before call : ',rho1n 569 rho1n = ddot(nt1amx*nt1amx,t2am,1,t2am,1) 570 write(lupri,*) 'Norm of T2AM before call : ',rho1n 571 rho1n = ddot(norbt,work(kfockd),1,work(kfockd),1) 572 write(lupri,*) 'Norm of KFOCKD before call : ',rho1n 573 rho1n = ddot(nt1amx,work(kscr1),1,work(kscr1),1) 574 write(lupri,*) 'Norm of KSCR1 before call : ',rho1n 575 call flshfo(lupri) 576C 577C CALL AROUND( 'In CC_FOP3: Triples Fock MO matrix' ) 578C CALL CC_PRFCKMO(WORK(KFOCK),ISYMOP) 579C do a = 1, nvirt 580C do i = 1, nrhft 581C nai = nvirt*(I-1) + a 582C if (abs(t1am(nai)) .gt. 1.0D-11) then 583C write(lupri,*) 'T1AM^{',a,'}_{',i,'} = ',t1am(nai) 584C endif 585C enddo 586C enddo 587COMMENT COMMENT COMMENT 588COMMENT COMMENT COMMENT 589C 590C 591 CALL CCSDT_T3BAR(WORK(KT3BAR),WORK(KINT3S),WORK(KINT4S), 592 * WORK(KIAJB),WORK(KFOCK),T1AM,T2AM, 593 * WORK(KSCR1),WORK(KFOCKD)) 594C 595 call temppr(work(kt3bar),2) 596C 597 rho1n = ddot(nt1amx*nt1amx*nt1amx,work(kt3am),1,work(kt3am),1) 598 write(lupri,*) 'Norm of T3AM after T3BAR = ',rho1n 599 rho1n = ddot(nt1amx*nt1amx*nt1amx,work(kt3bar),1,work(kt3bar),1) 600 write(lupri,*) 'Norm of T3BAR after T3BAR = ',rho1n 601C 602C----------------------------------------------------------- 603C Calculate the CCSD(T) corrections to the density. 604C----------------------------------------------------------- 605C 606 CALL DZERO(WORK(KOME1),NT1AMX) 607C ia block 608 CALL CCSDT_OMEGA1(WORK(KOME1),T2AM,WORK(KT3AM)) 609C 610COMMENT COMMENT COMMENT 611COMMENT COMMENT COMMENT 612 call around('1e- density in NODDY') 613 call cc_prp(work(kome1),DUMMY,1,1,0) 614COMMENT COMMENT COMMENT 615COMMENT COMMENT COMMENT 616 DO I = 1,NT1AMX 617 OMEGA1(I) = OMEGA1(I) + WORK(KOME1+I-1) 618 ENDDO 619C------------------------------------------------------ 620C Calculate the term to the unrelaxed gradient. 621C------------------------------------------------------ 622C 623COMMENT COMMENT 624 write(lupri,*) 'Field before contraction (fop3)' 625 call output(work(konep2),1,norbt,1,norbt,norbt,norbt,1,lupri) 626COMMENT COMMENT 627 CALL CCSDT_UNRE(T2AM,WORK(KT3AM),WORK(KT3BAR),WORK(KOME1), 628 * WORK(KONEP2)) 629C 630C ij and ab block 631C CALL CCSDT_NONREL_TERM(DENAB,DENIJ,WORK(KT3AM),WORK(KT3BAR)) 632C----------------------------------------------------------- 633C For debugging purposes calculate the right hand side 634C of the equations for the kappa amplitudes. 635C----------------------------------------------------------- 636C 637 CALL DZERO(WORK(KOME2),NT1AMX*NT1AMX) 638C 639 rho1n = ddot(nt1amx*nt1amx*nt1amx,work(kt3am),1,work(kt3am),1) 640 write(lupri,*) 'Norm of T3AM before kappa : ',rho1n 641 rho1n = ddot(nt1amx*nt1amx*nt1amx,work(kt3bar),1,work(kt3bar),1) 642 write(lupri,*) 'Norm of T3BAR before kappa : ',rho1n 643 rho1n = ddot(nt1amx,t1am,1,t1am,1) 644 write(lupri,*) 'Norm of T1AM before kappa : ',rho1n 645 rho1n = ddot(nt1amx*nt1amx,t2am,1,t2am,1) 646 write(lupri,*) 'Norm of T2AM before kappa : ',rho1n 647 rho1n = ddot(nvirt*nvirt*nvirt*nrhft,work(kvirt),1,work(kvirt),1) 648 write(lupri,*) 'Norm of VIRT before call : ',rho1n 649 rho1n = ddot(nvirt*nrhft*nrhft*nrhft,work(kocc),1,work(kocc),1) 650 write(lupri,*) 'Norm of OCC before kappa : ',rho1n 651 rho1n = ddot(nt1amx*nt1amx,work(kome2),1,work(kome2),1) 652 write(lupri,*) 'Norm of OMEGA2 before kappa: ',rho1n 653C 654 CALL CCSDT_KAPPA(WORK(KT3AM),WORK(KT3BAR),T1AM,T2AM, 655 * WORK(KVIRT),WORK(KOCC),WORK(KOME2), 656 * WORK(KOME22)) 657C 658C----------------------------------------------------------- 659C For debugging purposes calculate the right hand side 660C of the equations for the t1-bar and t2-bar 661C amplitudes. 662C----------------------------------------------------------- 663C 664C CALL CCSDT_OMEGA1OLD(WORK(KOME11),WORK(KIAJB),WORK(KT3AM)) 665C 666COMMENT COMMENT 667COMMENT COMMENT 668C do i = 1, nt1amx 669C if (abs(work(kome11-1+i)) .gt. 1.0D-9) then 670C write(lupri,*) 'Omega1_eta_debug(',i,') = ',work(kome11-1+i) 671C endif 672C enddo 673COMMENT COMMENT 674COMMENT COMMENT 675C 676 CALL DZERO(WORK(KOME2),NT1AMX*NT1AMX) 677 CALL DZERO(WORK(KOME22),NT1AMX*NT1AMX) 678 CALL CCETA_OMEGA2(WORK(KOME22),WORK(KOME2),WORK(KT3AM), 679 * WORK(KT3BAR),WORK(KINT1T),WORK(KINT2T), 680 * WORK(KFOCK)) 681C 682COMMENT COMMENT COMMENT COMMENT COMMENT 683COMMENT COMMENT COMMENT COMMENT COMMENT 684C call around('OMEGA2_ETA in NODDY :') 685C call output(work(kome22),1,nt1amx,1,nt1amx,nt1amx,nt1amx,1,lupri) 686COMMENT COMMENT COMMENT COMMENT COMMENT 687COMMENT COMMENT COMMENT COMMENT COMMENT 688C----------------------------------------------------------- 689C For debugging purposes calculate the diagonal 690C of the kappa amplitudes 691C----------------------------------------------------------- 692C 693 CALL DZERO(WORK(KKAPAA),NVIRT) 694 CALL DZERO(WORK(KKAPII),NRHFT) 695 CALL DEBUG_KAPPADIAG(WORK(KKAPAA),WORK(KKAPII), 696 * WORK(KT3AM),WORK(KT3BAR)) 697C 698 DO I = 1, NVIRT 699C IF (ABS(WORK(KKAPAA+I-1)) .GT. 1.0D-9) THEN 700 WRITE(LUPRI,*) 'Kappa_{aa}(',i,') = ',WORK(KKAPAA+I-1) 701C ENDIF 702 ENDDO 703 DO I = 1, NRHFT 704C IF (ABS(WORK(KKAPII+I-1)) .GT. 1.0D-9) THEN 705 WRITE(LUPRI,*) 'Kappa_{ii}(',i,') = ',WORK(KKAPII+I-1) 706C ENDIF 707 ENDDO 708C 709 RETURN 710 END 711C /* Deck ccsdt_tran1 */ 712 SUBROUTINE CCSDT_TRAN1(XINT1,XINT2,XLAMDP,XLAMDH,AOINT,IDEL) 713C 714C 715C 716#include "implicit.h" 717 DIMENSION XINT1(NT1AMX,NVIRT,NVIRT), XINT2(NT1AMX,NRHFT,NRHFT) 718 DIMENSION XLAMDP(NBAST,NORBT), XLAMDH(NBAST,NORBT) 719 DIMENSION AOINT(NNBAST,NBAST) 720#include "priunit.h" 721COMMENT COMMENT COMMENT 722C#include "inforb.h" 723#include "ccorb.h" 724COMMENT COMMENT COMMENT 725#include "ccsdinp.h" 726#include "ccsdsym.h" 727C 728 INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 729COMMENT COMMENT 730COMMENT COMMENT 731 write(lupri,*) 'NRHFT, NVIRT, NT1AMX : ',nrhft,nvirt,nt1amx 732 write(lupri,*) 'NORBT, NBAST : ',norbt,nbast 733COMMENT COMMENT 734COMMENT COMMENT 735C 736 DO 100 G = 1,NBAST 737 DO 110 IB = 1,NBAST 738 DO 120 A = 1,NBAST 739 NAB = INDEX(A,IB) 740C 741 if (aoint(nab,g) .eq. 0.0d0) goto 120 742 DO 200 D = 1,NVIRT 743 DO 210 B = 1,NVIRT 744 DO 220 K = 1,NRHFT 745 DO 230 C = 1,NVIRT 746C 747 NCK = NVIRT*(K-1) + C 748C 749 XINT1(NCK,B,D) = XINT1(NCK,B,D) 750 * + AOINT(NAB,G)*XLAMDH(A,NRHFT+C)*XLAMDP(IB,K) 751 * *XLAMDP(G,NRHFT+B)*XLAMDH(IDEL,NRHFT+D) 752C 753 230 CONTINUE 754 220 CONTINUE 755 210 CONTINUE 756 200 CONTINUE 757C 758 DO 300 J = 1,NRHFT 759 DO 310 L = 1,NRHFT 760 DO 320 K = 1,NRHFT 761 DO 330 C = 1,NVIRT 762C 763 NCK = NVIRT*(K-1) + C 764C 765 XINT2(NCK,L,J) = XINT2(NCK,L,J) 766 * + AOINT(NAB,G)*XLAMDH(A,NRHFT+C)*XLAMDP(IB,K) 767 * *XLAMDP(G,L)*XLAMDH(IDEL,J) 768C 769 330 CONTINUE 770 320 CONTINUE 771 310 CONTINUE 772 300 CONTINUE 773C 774 120 CONTINUE 775 110 CONTINUE 776 100 CONTINUE 777C 778 RETURN 779 END 780C /* Deck ccsdt_tran2 */ 781 SUBROUTINE CCSDT_TRAN2(XIAJB,YIAJB,CMO,AOINT,IDEL) 782C 783C 784C 785#include "implicit.h" 786 PARAMETER (TWO = 2.0D0) 787 DIMENSION XIAJB(NT1AMX,NT1AMX), AOINT(NNBAST,NBAST) 788 DIMENSION YIAJB(NT1AMX,NT1AMX) 789 DIMENSION CMO(NBAST,NORBT) 790#include "priunit.h" 791COMMENT COMMENT COMMENT 792C#include "inforb.h" 793#include "ccorb.h" 794COMMENT COMMENT COMMENT 795#include "ccsdinp.h" 796#include "ccsdsym.h" 797C 798 INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 799C 800 DO 100 G = 1,NBAST 801 DO 110 B = 1,NBAST 802 DO 120 A = 1,NBAST 803 NAB = INDEX(A,B) 804C 805 if (aoint(nab,g) .eq. 0.0d0) goto 120 806 DO 200 L = 1,NRHFT 807 DO 210 D = 1,NVIRT 808 NDL = NVIRT*(L-1) + D 809 DO 220 K = 1,NRHFT 810 DO 230 C = 1,NVIRT 811C 812 NCK = NVIRT*(K-1) + C 813C 814 XIAJB(NCK,NDL) = XIAJB(NCK,NDL) 815 * + AOINT(NAB,G)* 816 * (TWO*CMO(A,NRHFT+C)*CMO(B,K)*CMO(G,NRHFT+D)*CMO(IDEL,L) 817 * - CMO(A,NRHFT+C)*CMO(B,L)*CMO(G,NRHFT+D)*CMO(IDEL,K)) 818C 819 YIAJB(NCK,NDL) = YIAJB(NCK,NDL) 820 * + AOINT(NAB,G)* 821 * CMO(A,NRHFT+C)*CMO(B,K)*CMO(G,NRHFT+D)*CMO(IDEL,L) 822C 823 230 CONTINUE 824 220 CONTINUE 825 210 CONTINUE 826 200 CONTINUE 827C 828 120 CONTINUE 829 110 CONTINUE 830 100 CONTINUE 831C 832 RETURN 833 END 834C /* Deck ccsdt_t3am */ 835 SUBROUTINE CCSDT_T3AM(T3AM,XINT1,XINT2,T2AM,SCR1,FOCKD, 836 * T3AM2,FCKFIELD,ITERATIVE) 837C 838C 839C 840#include "implicit.h" 841 DIMENSION XINT1(NT1AMX,NVIRT,NVIRT), XINT2(NT1AMX,NRHFT,NRHFT) 842 DIMENSION T3AM(NT1AMX,NT1AMX,NT1AMX),SCR1(NT1AMX),FOCKD(*) 843 DIMENSION T3AM2(NT1AMX,NT1AMX,NT1AMX) 844 DIMENSION T2AM(NT1AMX,NT1AMX) 845 DIMENSION FCKFIELD(NORBT,NORBT) 846#include "priunit.h" 847COMMENT COMMENT COMMENT 848C#include "inforb.h" 849#include "ccorb.h" 850COMMENT COMMENT COMMENT 851#include "ccsdinp.h" 852#include "ccsdsym.h" 853#include "ccfield.h" 854C 855 PARAMETER(ONE = 1.0D0, HALF = 0.5D0) 856C 857 LOGICAL LDBG1, LDBG2, FIRST, NONCONV, ITERATIVE 858C 859 INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 860C 861 DO 50 I = 1,NRHFT 862 DO 60 A = 1,NVIRT 863 NAI = NVIRT*(I-1) + A 864 SCR1(NAI) = FOCKD(NRHFT+A) - FOCKD(I) 865 60 CONTINUE 866 50 CONTINUE 867C 868 IF (ITERATIVE) THEN 869 CALL DZERO(T3AM2,NT1AMX*NT1AMX*NT1AMX) 870 ENDIF 871C 872 FIRST = .TRUE. 873 NONCONV = .TRUE. 874 MAXITE = 0 875C 876COMMENT COMMENT 877 write(lupri,*) 'ITERATIVE = ',iterative 878COMMENT COMMENT 879 DO WHILE ((FIRST) .OR. 880 * ((ITERATIVE) .AND. (NONCONV))) 881C 882 FIRST = .FALSE. 883C 884 DO NCK = 1,NT1AMX 885C 886 DO J = 1,NRHFT 887 DO B = 1,NVIRT 888C 889 NBJ = NVIRT*(J-1) + B 890C 891 DO NAI = 1,NT1AMX 892C 893 AIBJCK = 0.0D0 894 DO D = 1,NVIRT 895C 896 NDJ = NVIRT*(J-1) + D 897C 898 AIBJCK = AIBJCK + XINT1(NCK,B,D)*T2AM(NDJ,NAI) 899C 900 ENDDO 901C 902 DO L = 1,NRHFT 903C 904 NBL = NVIRT*(L-1) + B 905C 906 AIBJCK = AIBJCK - XINT2(NCK,L,J)*T2AM(NBL,NAI) 907COMMENT COMMENT 908COMMENT COMMENT 909C if (abs(XINT2(NCK,L,J)*T2AM(NBL,NAI)) .gt. 1.0D-12) then 910C write(lupri,'(A,5I3)') 'NAI, B, J, NCK, L : ', 911C * NAI, B, J, NCK, L 912C write(lupri,*) 'INT(NCK.L.J) T2AM(NBL,NAI) : ', 913C * XINT2(NCK,L,J),T2AM(NBL,NAI) 914C endif 915COMMENT COMMENT 916COMMENT COMMENT 917C 918 ENDDO 919C 920C--------------------------------------------------------------------- 921C Add field dependent terms 922C N^8 scaling but who cares in a noddy!!! 923C Divide the first multiplication in two to get N^7 scaling 924C--------------------------------------------------------------------- 925C 926 IF (ITERATIVE) THEN 927C 928 DO D = 1, NVIRT 929 NDJ = NVIRT*(J-1) + D 930 DO L = 1, NRHFT 931 NBL = NVIRT*(L-1) + B 932COMMENT COMMENT 933 AIBJCK = AIBJCK 934 * + T2AM(NAI,NBL) 935 * *T2AM(NCK,NDJ) 936 * *FCKFIELD(L,NRHFT+D) 937COMMENT COMMENT 938 ENDDO 939 ENDDO 940C 941 DO D = 1, NVIRT 942 NDJ = NVIRT*(J-1) + D 943COMMENT COMMENT 944 AIBJCK = AIBJCK 945 * - HALF*T3AM2(NAI,NDJ,NCK) 946 * *FCKFIELD(NRHFT+B,NRHFT+D) 947COMMENT COMMENT 948 ENDDO 949C 950 DO L = 1, NRHFT 951 NBL = NVIRT*(L-1) + B 952COMMENT COMMENT 953 AIBJCK = AIBJCK 954 * + HALF*T3AM2(NAI,NBL,NCK) 955 * *FCKFIELD(L,J) 956COMMENT COMMENT 957 ENDDO 958C 959 ENDIF 960C 961C----------------------------------- 962C Final summation 963C----------------------------------- 964C 965 AIBJCK = AIBJCK/(SCR1(NAI) + SCR1(NBJ) + SCR1(NCK)) 966COMMENT COMMENT 967COMMENT COMMENT 968C if (abs(aibjck) .gt. 1.0D-12) then 969C write(lupri,'(A,e20.10)') 'Contri : ',AIBJCK 970C endif 971COMMENT COMMENT 972COMMENT COMMENT 973C 974 T3AM(NAI,NBJ,NCK) = T3AM(NAI,NBJ,NCK) + AIBJCK 975 T3AM(NAI,NCK,NBJ) = T3AM(NAI,NCK,NBJ) + AIBJCK 976 T3AM(NBJ,NAI,NCK) = T3AM(NBJ,NAI,NCK) + AIBJCK 977 T3AM(NCK,NAI,NBJ) = T3AM(NCK,NAI,NBJ) + AIBJCK 978 T3AM(NBJ,NCK,NAI) = T3AM(NBJ,NCK,NAI) + AIBJCK 979 T3AM(NCK,NBJ,NAI) = T3AM(NCK,NBJ,NAI) + AIBJCK 980C 981 ENDDO ! NAI 982 ENDDO ! B 983 ENDDO ! J 984 ENDDO ! NCK 985C 986C---------------------------------------- 987C Remove the forbidden elements. 988C---------------------------------------- 989C 990 do a = 1, nvirt 991 do i = 1, nrhft 992 nai = nvirt*(i-1) + a 993 do b = 1, nvirt 994 nbi = nvirt*(i-1) + b 995 do c = 1, nvirt 996 nci = nvirt*(i-1) + c 997C 998 t3am(nai,nbi,nci) = 0.0D0 999 t3am(nai,nci,nbi) = 0.0D0 1000 t3am(nbi,nai,nci) = 0.0D0 1001 t3am(nbi,nci,nai) = 0.0D0 1002 t3am(nci,nai,nbi) = 0.0D0 1003 t3am(nci,nbi,nai) = 0.0D0 1004 enddo 1005 enddo 1006C 1007 do j = 1, nrhft 1008 naj = nvirt*(j-1) + a 1009 do k = 1, nrhft 1010 nak = nvirt*(k-1) + a 1011C 1012 t3am(nai,naj,nak) = 0.0D0 1013 t3am(nai,nak,naj) = 0.0D0 1014 t3am(naj,nai,nak) = 0.0D0 1015 t3am(naj,nak,nai) = 0.0D0 1016 t3am(nak,nai,naj) = 0.0D0 1017 t3am(nak,naj,nai) = 0.0D0 1018 enddo 1019 enddo 1020C 1021 enddo 1022 enddo 1023C 1024C------------------------------------------------------------ 1025C If iterative field thing compare the new and the old 1026C T3 amplitudes and see if they differ. 1027C------------------------------------------------------------ 1028C 1029 IF (ITERATIVE) THEN 1030C 1031 CALL DAXPY(NT1AMX*NT1AMX*NT1AMX,-ONE,T3AM,1,T3AM2,1) 1032C 1033 XNORM1 = DDOT(NT1AMX*NT1AMX*NT1AMX,T3AM2,1,T3AM2,1) 1034C 1035 XNORM = SQRT(XNORM1) 1036COMMENT COMMENT 1037 write(lupri,*) 'Norm of difference in T3AM = ',XNORM 1038COMMENT COMMENT 1039C 1040 IF (XNORM .LT. 1.0D-10) THEN 1041 NONCONV = .FALSE. 1042 ELSE 1043 CALL DCOPY(NT1AMX*NT1AMX*NT1AMX,T3AM,1,T3AM2,1) 1044 CALL DZERO(T3AM,NT1AMX*NT1AMX*NT1AMX) 1045 MAXITE = MAXITE + 1 1046 IF (MAXITE .GT. 30) THEN 1047 CALL QUIT('MAX ITERATIONS EXCEEDED in CCSDT_T3AM') 1048 ENDIF 1049 ENDIF 1050 ENDIF 1051C 1052C-------------------------- 1053C End do while loop 1054C-------------------------- 1055C 1056 ENDDO 1057C 1058C 1059 RETURN 1060 END 1061C /* Deck ccsdt_omega1 */ 1062 SUBROUTINE CCSDT_OMEGA1(OMEGA1,T2AM,T3AM) 1063C 1064C 1065C 1066#include "implicit.h" 1067 PARAMETER (TWO = 2.0D0) 1068 DIMENSION OMEGA1(NT1AMX),T2AM(NT1AMX,NT1AMX) 1069 DIMENSION T3AM(NT1AMX,NT1AMX,NT1AMX) 1070#include "priunit.h" 1071COMMENT COMMENT COMMENT 1072C#include "inforb.h" 1073#include "ccorb.h" 1074COMMENT COMMENT COMMENT 1075#include "ccsdinp.h" 1076#include "ccsdsym.h" 1077C 1078 INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 1079C 1080 DO 100 I = 1,NRHFT 1081 DO 110 A = 1,NVIRT 1082 NAI = NVIRT*(I-1) + A 1083C 1084 DO 115 B = 1,NVIRT 1085 DO 120 J = 1,NRHFT 1086 NBJ = NVIRT*(J-1) + B 1087C 1088 DO 130 K = 1,NRHFT 1089 NAK = NVIRT*(K-1) + A 1090 NBK = NVIRT*(K-1) + B 1091 DO 140 C = 1,NVIRT 1092 NCK = NVIRT*(K-1) + C 1093 NCI = NVIRT*(I-1) + C 1094 NCJ = NVIRT*(J-1) + C 1095C 1096 XTEMP = TWO* 1097 * (TWO*T2AM(NCK,NBJ)-T2AM(NCJ,NBK))* 1098 * ( T3AM(NAI,NBJ,NCK) - T3AM(NAK,NBJ,NCI) ) 1099C 1100 OMEGA1(NAI) = OMEGA1(NAI) + XTEMP 1101C 1102COMMENT COMMENT 1103COMMENT COMMENT 1104C if (abs(xtemp) .gt. 1.0D-9) then 1105C write(lupri,*) 'Cont to NAI (',nai,') = ',xtemp 1106C if (i .EQ. 1 .AND. a .EQ. 3) then 1107C write(lupri,*) 'A, B, C, I, J, K = ',a,b,c,i,j,k 1108C write(lupri,*) 'NAI, NBJ, NCK = ',nai,nbj,nck 1109C write(lupri,*) 'NAK, NCI = ',nak,nci 1110C write(lupri,*) 'TWO*T2AM(CK,BJ) = ',TWO*T2AM(NCK,NBJ) 1111C write(lupri,*) 'T2AM(CJ,BK) = ',T2AM(NCJ,NBK) 1112C write(lupri,*) 'T3AM(AI,BJ,CK) = ',T3AM(NAI,NBJ,NCK) 1113C write(lupri,*) 'T3AM(AK,BJ,CI) = ', T3AM(NAK,NBJ,NCI) 1114C endif 1115C endif 1116COMMENT COMMENT 1117COMMENT COMMENT 1118C 1119C 1120 140 CONTINUE 1121 130 CONTINUE 1122C 1123 120 CONTINUE 1124 115 CONTINUE 1125C 1126 110 CONTINUE 1127 100 CONTINUE 1128C 1129 RETURN 1130 END 1131C /* Deck ccsdt_tran3 */ 1132 SUBROUTINE CCSDT_TRAN3(XINT1,XINT2,XLAMDP,XLAMDH,AOINT,IDEL) 1133C 1134C 1135C 1136#include "implicit.h" 1137 DIMENSION XINT1(NT1AMX,NVIRT,NVIRT), XINT2(NT1AMX,NRHFT,NRHFT) 1138 DIMENSION XLAMDP(NBAST,NORBT), XLAMDH(NBAST,NORBT) 1139 DIMENSION AOINT(NNBAST,NBAST) 1140#include "priunit.h" 1141COMMENT COMMENT COMMENT 1142C#include "inforb.h" 1143#include "ccorb.h" 1144COMMENT COMMENT COMMENT 1145#include "ccsdinp.h" 1146#include "ccsdsym.h" 1147C 1148 LOGICAL LDEBUG 1149C 1150 INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 1151C 1152 LDEBUG = .FALSE. 1153C 1154 DO 100 G = 1,NBAST 1155 DO 110 IB = 1,NBAST 1156 DO 120 A = 1,NBAST 1157 NAB = INDEX(A,IB) 1158C 1159 if (aoint(nab,g) .eq. 0.0d0) goto 120 1160 DO 200 D = 1,NVIRT 1161 DO 210 B = 1,NVIRT 1162 DO 220 K = 1,NRHFT 1163 DO 230 C = 1,NVIRT 1164C 1165 NCK = NVIRT*(K-1) + C 1166C 1167 XINT1(NCK,B,D) = XINT1(NCK,B,D) 1168 * + AOINT(NAB,G)*XLAMDP(A,NRHFT+C)*XLAMDH(IB,K) 1169 * *XLAMDP(G,NRHFT+B)*XLAMDH(IDEL,NRHFT+D) 1170C 1171 230 CONTINUE 1172 220 CONTINUE 1173 210 CONTINUE 1174 200 CONTINUE 1175C 1176 DO 300 J = 1,NRHFT 1177 DO 310 L = 1,NRHFT 1178 DO 320 K = 1,NRHFT 1179 DO 330 C = 1,NVIRT 1180C 1181 NCK = NVIRT*(K-1) + C 1182C 1183 XINT2(NCK,L,J) = XINT2(NCK,L,J) 1184 * + AOINT(NAB,G)*XLAMDP(A,NRHFT+C)*XLAMDH(IB,K) 1185 * *XLAMDP(G,L)*XLAMDH(IDEL,J) 1186C 1187 330 CONTINUE 1188 320 CONTINUE 1189 310 CONTINUE 1190 300 CONTINUE 1191C 1192 120 CONTINUE 1193 110 CONTINUE 1194 100 CONTINUE 1195C 1196 IF (LDEBUG) THEN 1197 do a = 1, nt1amx 1198 do b = 1, nrhtf 1199 do c = 1, nrhtf 1200 if (abs(xint1(a,b,c)) .gt. 1.0D-9) then 1201 write(lupri,*) 'Lam-int1(',a,',',b,',',c,') = ', 1202 * xint1(a,b,c) 1203 endif 1204 enddo 1205 enddo 1206 enddo 1207 do a = 1, nt1amx 1208 do b = 1, nvirt 1209 do c = 1, nvirt 1210 if (abs(xint2(a,b,c)) .gt. 1.0D-9) then 1211 write(lupri,*) 'Lam-int2(',a,',',b,',',c,') = ', 1212 * xint2(a,b,c) 1213 endif 1214 enddo 1215 enddo 1216 enddo 1217 ENDIF 1218C 1219 RETURN 1220 END 1221C /* Deck ccsdt_tran3pt */ 1222 SUBROUTINE CCSDT_TRAN3PT(XINT1,XINT2,CMO,AOINT,IDEL) 1223C 1224C 1225C 1226#include "implicit.h" 1227#include "priunit.h" 1228COMMENT COMMENT COMMENT 1229C#include "inforb.h" 1230#include "ccorb.h" 1231COMMENT COMMENT COMMENT 1232#include "ccsdinp.h" 1233#include "ccsdsym.h" 1234 DIMENSION XINT1(NT1AMX,NVIRT,NVIRT), XINT2(NT1AMX,NRHFT,NRHFT) 1235 DIMENSION CMO(NBAST,NORBT) 1236 DIMENSION AOINT(NNBAST,NBAST) 1237C 1238 INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 1239C 1240 DO 100 G = 1,NBAST 1241 DO 110 IB = 1,NBAST 1242 DO 120 A = 1,NBAST 1243 NAB = INDEX(A,IB) 1244C 1245 if (aoint(nab,g) .eq. 0.0d0) goto 120 1246 DO 200 D = 1,NVIRT 1247 DO 210 B = 1,NVIRT 1248 DO 220 K = 1,NRHFT 1249 DO 230 C = 1,NVIRT 1250C 1251 NCK = NVIRT*(K-1) + C 1252C 1253 XINT1(NCK,B,D) = XINT1(NCK,B,D) 1254 * + AOINT(NAB,G)*CMO(A,NRHFT+C)*CMO(IB,K) 1255 * *CMO(G,NRHFT+B)*CMO(IDEL,NRHFT+D) 1256C 1257 230 CONTINUE 1258 220 CONTINUE 1259 210 CONTINUE 1260 200 CONTINUE 1261C 1262 DO 300 J = 1,NRHFT 1263 DO 310 L = 1,NRHFT 1264 DO 320 K = 1,NRHFT 1265 DO 330 C = 1,NVIRT 1266C 1267 NCK = NVIRT*(K-1) + C 1268C 1269 XINT2(NCK,L,J) = XINT2(NCK,L,J) 1270 * + AOINT(NAB,G)*CMO(A,NRHFT+C)*CMO(IB,K) 1271 * *CMO(G,L)*CMO(IDEL,J) 1272C 1273 330 CONTINUE 1274 320 CONTINUE 1275 310 CONTINUE 1276 300 CONTINUE 1277C 1278 120 CONTINUE 1279 110 CONTINUE 1280 100 CONTINUE 1281C 1282 RETURN 1283 END 1284C /* Deck ccsdt_t3bar */ 1285 SUBROUTINE CCSDT_T3BAR(T3BAR,XINT1,XINT2,XIAJB,FOCK,T1AM, 1286 * T2AM,SCR1,FOCKD) 1287C 1288C 1289C 1290#include "implicit.h" 1291C 1292 PARAMETER (TWO = 2.0D0) 1293C 1294 DIMENSION XINT1(NT1AMX, NVIRT,NVIRT), XINT2(NT1AMX, NRHFT,NRHFT) 1295 DIMENSION XIAJB(NT1AMX, NT1AMX), FOCK(NORBT,NORBT) 1296 DIMENSION T3BAR(NT1AMX, NT1AMX,NT1AMX), SCR1(NT1AMX), FOCKD(*) 1297 DIMENSION T1AM(NT1AMX), T2AM(NT1AMX,NT1AMX) 1298 LOGICAL LDBG1, LDBG2 1299#include "priunit.h" 1300COMMENT COMMENT COMMENT 1301C#include "inforb.h" 1302#include "ccorb.h" 1303COMMENT COMMENT COMMENT 1304#include "ccsdinp.h" 1305#include "ccsdsym.h" 1306C 1307 INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 1308C 1309 DO 50 NI = 1,NRHFT 1310 DO 60 NA = 1,NVIRT 1311 NAI = NVIRT*(NI-1) + NA 1312 SCR1(NAI) = FOCKD(NRHFT+NA) - FOCKD(NI) 1313 60 CONTINUE 1314 50 CONTINUE 1315COMMENT COMMENT COMMENT COMMENT 1316COMMENT COMMENT COMMENT COMMENT 1317COMMENT COMMENT COMMENT COMMENT 1318COMMENT COMMENT COMMENT COMMENT 1319C do a = 1, nvirt 1320C do i = 1, nrhft 1321C nai = nvirt*(i-1) + a 1322C do b = 1, nvirt 1323C do j = 1, nvirt 1324C nbj = nvirt*(j-1) + b 1325CC 1326C if (abs(t2am(nai,nbj)) .gt. 1.0D-10) then 1327C write(lupri,*) 'T2AM^{',a,',',b,'}_{',i,',',j,'}) = ', 1328C * t2am(nai,nbj) 1329C endif 1330CC 1331C enddo 1332C enddo 1333C enddo 1334C enddo 1335C 1336C do ni = 1, nt1amx 1337C if (abs(scr1(ni)) .gt. 1.0D-9) then 1338C write(lupri,*) 'SCR1(',ni,') = ',SCR1(ni) 1339C endif 1340C enddo 1341C 1342C do na = 1, nvirt 1343C do ni = 1, nrhft 1344C do nb = 1, nvirt 1345C do nc = 1, nvirt 1346CC 1347C nai = NVIRT*(NI-1) + NA 1348C nbi = NVIRT*(NI-1) + NB 1349CC 1350C if (abs(TWO*XINT1(nai,nb,nc)-XINT1(nbi,na,nc)) 1351C * .gt. 1.0D-9) then 1352C write(lupri,*)'L_vir(',nai,',',nb,',',nc,') =', 1353C * TWO*XINT1(nai,nb,nc)-XINT1(nbi,na,nc) 1354C endif 1355C enddo 1356C enddo 1357C enddo 1358C enddo 1359C do na = 1, nvirt 1360C do ni = 1, nrhft 1361C do nb = 1, nvirt 1362C do nc = 1, nvirt 1363CC 1364C nai = NVIRT*(NI-1) + NA 1365C nbi = NVIRT*(NI-1) + NB 1366CC 1367C if (abs(XINT1(nai,nb,nc)) .gt. 1.0D-9) then 1368C write(lupri,*)'g_vir(',nai,',',nb,',',nc,') =', 1369C * XINT1(nai,nb,nc) 1370C endif 1371C enddo 1372C enddo 1373C enddo 1374C enddo 1375CC 1376CC 1377C xnorm = 0.0D0 1378C do na = 1, nvirt 1379C do ni = 1, nrhft 1380C do nj = 1, nrhft 1381C do nk = 1, nrhft 1382CC 1383C nai = NVIRT*(NI-1) + NA 1384C nak = NVIRT*(NK-1) + NA 1385CC 1386C xnorm = xnorm + (TWO*XINT2(nai,nj,nk)-XINT2(nak,nj,ni)) 1387C * *(TWO*XINT2(nai,nj,nk)-XINT2(nak,nj,ni)) 1388CC 1389C if (abs(TWO*XINT2(nai,nj,nk)-XINT2(nak,nj,ni)) 1390C * .gt. 1.0D-9) then 1391C write(lupri,*)'L_occ(',nai,',',nj,',',nk,') =', 1392C * TWO*XINT2(nai,nj,nk)-XINT2(nak,nj,ni) 1393C endif 1394C enddo 1395C enddo 1396C enddo 1397C enddo 1398CC 1399C write(lupri,*) 'Norm squared of the L_occ : ',xnorm 1400C do na = 1, nvirt 1401C do ni = 1, nrhft 1402C do nj = 1, nrhft 1403C do nk = 1, nrhft 1404CC 1405C nai = NVIRT*(NI-1) + NA 1406CC 1407C if (abs(XINT2(nai,nj,nk)) .gt. 1.0D-9) then 1408C write(lupri,*)'g_occ(',nai,',',nj,',',nk,') =', 1409C * XINT2(nai,nj,nk) 1410C endif 1411C enddo 1412C enddo 1413C enddo 1414C enddo 1415COMMENT COMMENT COMMENT COMMENT 1416COMMENT COMMENT COMMENT COMMENT 1417COMMENT COMMENT COMMENT COMMENT 1418COMMENT COMMENT COMMENT COMMENT 1419C 1420 DO 100 NK = 1,NRHFT 1421 DO 105 NC = 1,NVIRT 1422C 1423 NCK = NVIRT*(NK-1) + NC 1424C 1425 DO 110 NJ = 1,NRHFT 1426 NCJ = NVIRT*(NJ-1) + NC 1427 DO 120 NB = 1,NVIRT 1428C 1429 NBJ = NVIRT*(NJ-1) + NB 1430 NBK = NVIRT*(NK-1) + NB 1431C 1432 DO 130 NI = 1,NRHFT 1433 NBI = NVIRT*(NI-1) + NB 1434 NCI = NVIRT*(NI-1) + NC 1435 DO 135 NA = 1,NVIRT 1436C 1437 NAI = NVIRT*(NI-1) + NA 1438 NAJ = NVIRT*(NJ-1) + NA 1439 NAK = NVIRT*(NK-1) + NA 1440C 1441 AIBJCK = 0.0D0 1442C 1443 if (.true.) then 1444C T1* = TWO T1 => Factor of two 1445 AIBJCK = AIBJCK + TWO*XIAJB(NBJ,NCK)*T1AM(NAI) 1446COMMENT COMMENT 1447COMMENT COMMENT 1448C if (abs(XIAJB(NBJ,NCK)*T1AM(NAI)) .gt. 1.0D-10) then 1449C write(lupri,*) 'Contribution from 1.st term : ', 1450C * XIAJB(NBJ,NCK)*T1AM(NAI) 1451C write(lupri,*) 'With NAI, NBJ, NCK : ',NAI,NBJ,NCK 1452C endif 1453COMMENT COMMENT 1454COMMENT COMMENT 1455 endif 1456C 1457 if (.true.) then 1458C T1* = TWO T1 => Factor of two 1459 AIBJCK = AIBJCK - TWO*XIAJB(NBJ,NCI)*T1AM(NAK) 1460COMMENT COMMENT 1461COMMENT COMMENT 1462C if (abs(XIAJB(NBJ,NCK)*T1AM(NAI)) .gt. 1.0D-10) then 1463C write(lupri,*) 'Contribution from 2.nd term : ', 1464C * XIAJB(NBJ,NCK)*T1AM(NAI) 1465C write(lupri,*) 'With NAI, NBJ, NCK : ',NAI,NBJ,NCK 1466C endif 1467COMMENT COMMENT 1468COMMENT COMMENT 1469 endif 1470C 1471 if (.true.) then 1472 AIBJCK = AIBJCK + TWO*FOCK(NK,NRHFT+NC)* 1473 * (TWO*T2AM(NAI,NBJ)-T2AM(NAJ,NBI)) 1474 endif 1475C 1476 if (.true.) then 1477 AIBJCK = AIBJCK - TWO*FOCK(NJ,NRHFT+NC)* 1478 * (TWO*T2AM(NAI,NBK)-T2AM(NAK,NBI)) 1479 endif 1480C 1481 DO 140 ND = 1,NVIRT 1482C 1483 NDI = NVIRT*(NI-1) + ND 1484 NDJ = NVIRT*(NJ-1) + ND 1485 NDK = NVIRT*(NK-1) + ND 1486C 1487 if (.true.) then 1488 AIBJCK = AIBJCK + TWO* 1489 * (TWO*XINT1(NCK,NB,ND)-XINT1(NBK,NC,ND))* 1490 * (TWO*T2AM(NDJ,NAI)-T2AM(NDI,NAJ)) 1491C 1492COMMENT COMMENT 1493COMMENT COMMENT 1494C if (abs((TWO*XINT1(NCK,NB,ND)-XINT1(NBK,NC,ND))* 1495C * (TWO*T2AM(NDJ,NAI)-T2AM(NDI,NAJ))) 1496C * .gt. 1.0D-10) then 1497C write(lupri,*) 'Contribution from 5.th term :', 1498C * TWO*(TWO*XINT1(NCK,NB,ND)-XINT1(NBK,NC,ND))* 1499C * (TWO*T2AM(NDJ,NAI)-T2AM(NDI,NAJ)) 1500C write(lupri,*) '2*XINT1, XINT1, 2*T2AM, T2AM', 1501C * TWO*XINT1(NCK,NB,ND),XINT1(NBK,NC,ND), 1502C * TWO*T2AM(NDJ,NAI),T2AM(NDI,NAJ) 1503C endif 1504COMMENT COMMENT 1505COMMENT COMMENT 1506 endif 1507C 1508 if (.true.) then 1509 AIBJCK = AIBJCK - TWO* 1510 * XINT1(NBI,NC,ND)* 1511 * (TWO*T2AM(NDK,NAJ)-T2AM(NDJ,NAK)) 1512COMMENT COMMENT 1513COMMENT COMMENT 1514C if (abs(TWO*XINT1(NBI,NC,ND)* 1515C * (TWO*T2AM(NDK,NAJ)-T2AM(NDJ,NAK))) 1516C * .gt. 1.0D-10) then 1517C write(lupri,*) 'Contribution from 6.th term :', 1518C * TWO*XINT1(NBI,NC,ND)* 1519C * (TWO*T2AM(NDK,NAJ)-T2AM(NDJ,NAK)) 1520C endif 1521COMMENT COMMENT 1522COMMENT COMMENT 1523 endif 1524C 1525 140 CONTINUE 1526C 1527 DO 150 NL = 1,NRHFT 1528C 1529 NAL = NVIRT*(NL-1) + NA 1530 NBL = NVIRT*(NL-1) + NB 1531C 1532 if (.true.) then 1533 AIBJCK = AIBJCK - TWO* 1534 * (TWO*XINT2(NCK,NL,NJ)-XINT2(NCJ,NL,NK))* 1535 * (TWO*T2AM(NBL,NAI)-T2AM(NBI,NAL)) 1536COMMENT COMMENT 1537COMMENT COMMENT 1538C if (abs(TWO* 1539C * (TWO*XINT2(NCK,NL,NJ)-XINT2(NCJ,NL,NK))* 1540C * (TWO*T2AM(NBL,NAI)-T2AM(NBI,NAL))) .gt. 1.0D-10) then 1541CC 1542C write(lupri,*) 'Contribution from L_occ term :',TWO* 1543C * (TWO*XINT2(NCK,NL,NJ)-XINT2(NCJ,NL,NK))* 1544C * (TWO*T2AM(NBL,NAI)-T2AM(NBI,NAL)) 1545C write(lupri,*) 'T2TCME, L : ', 1546C * (TWO*T2AM(NBL,NAI)-T2AM(NBI,NAL)), 1547C * (TWO*XINT2(NCK,NL,NJ)-XINT2(NCJ,NL,NK)) 1548C endif 1549COMMENT COMMENT 1550COMMENT COMMENT 1551 endif 1552C 1553 if (.true.) then 1554 AIBJCK = AIBJCK + TWO* 1555 * XINT2(NCJ,NL,NI)* 1556 * (TWO*T2AM(NBK,NAL)-T2AM(NBL,NAK)) 1557 endif 1558C 1559 150 CONTINUE 1560C 1561C 1562 AIBJCK = AIBJCK/(SCR1(NAI) + SCR1(NBJ) + SCR1(NCK)) 1563C 1564COMMENT COMMENT 1565COMMENT COMMENT 1566C if (abs(aibjck).gt.(1.0D-10/(scr1(nai)+scr1(nbj)+scr1(nck)))) 1567C * then 1568C write(lupri,'(1X,A,I3,I3,I3,I3,I3,I3)') 1569C * 'Contribution for A,B,C,I,J,K = ', 1570C * NA,NB,NC,NI,NJ,NK 1571C write(lupri,'(A,3I3,e20.10)') 'NAI,NBJ,NCK,AIBJCK = ', 1572C * NAI,NBJ,NCK,AIBJCK 1573CC write(lupri,*) 'SCR1(NAI),SCR1(NBJ),SCR1(NCK), SUM = ', 1574CC * SCR1(NAI),SCR1(NBJ),SCR1(NCK), 1575CC * SCR1(NAI)+SCR1(NBJ)+SCR1(NCK) 1576C write(lupri,*) ' ' 1577C endif 1578COMMENT COMMENT 1579COMMENT COMMENT 1580 T3BAR(NAI,NBJ,NCK) = T3BAR(NAI,NBJ,NCK) + AIBJCK 1581 T3BAR(NAI,NCK,NBJ) = T3BAR(NAI,NCK,NBJ) + AIBJCK 1582 T3BAR(NBJ,NAI,NCK) = T3BAR(NBJ,NAI,NCK) + AIBJCK 1583 T3BAR(NCK,NAI,NBJ) = T3BAR(NCK,NAI,NBJ) + AIBJCK 1584 T3BAR(NBJ,NCK,NAI) = T3BAR(NBJ,NCK,NAI) + AIBJCK 1585 T3BAR(NCK,NBJ,NAI) = T3BAR(NCK,NBJ,NAI) + AIBJCK 1586C 1587 135 CONTINUE 1588 130 CONTINUE 1589 120 CONTINUE 1590 110 CONTINUE 1591 105 CONTINUE 1592 100 CONTINUE 1593C 1594C------------------------------------------------------ 1595C Get rid of amplitudes that are not allowed. 1596C------------------------------------------------------ 1597C 1598 do a = 1, nvirt 1599 do i = 1, nrhft 1600 nai = nvirt*(i-1) + a 1601 do b = 1, nvirt 1602 nbi = nvirt*(i-1) + b 1603 do c = 1, nvirt 1604 nci = nvirt*(i-1) + c 1605C 1606 t3bar(nai,nbi,nci) = 0.0D0 1607 t3bar(nai,nci,nbi) = 0.0D0 1608 t3bar(nbi,nai,nci) = 0.0D0 1609 t3bar(nbi,nci,nai) = 0.0D0 1610 t3bar(nci,nai,nbi) = 0.0D0 1611 t3bar(nci,nbi,nai) = 0.0D0 1612 enddo 1613 enddo 1614C 1615 do j = 1, nrhft 1616 naj = nvirt*(j-1) + a 1617 do k = 1, nrhft 1618 nak = nvirt*(k-1) + a 1619C 1620 t3bar(nai,naj,nak) = 0.0D0 1621 t3bar(nai,nak,naj) = 0.0D0 1622 t3bar(naj,nai,nak) = 0.0D0 1623 t3bar(naj,nak,nai) = 0.0D0 1624 t3bar(nak,nai,naj) = 0.0D0 1625 t3bar(nak,naj,nai) = 0.0D0 1626 enddo 1627 enddo 1628C 1629 enddo 1630 enddo 1631C 1632 RETURN 1633 END 1634C /* Deck temppr */ 1635 SUBROUTINE TEMPPR(T3AM,IOPT) 1636C 1637#include "implicit.h" 1638 DIMENSION T3AM(NT1AMX,NT1AMX,NT1AMX) 1639#include "priunit.h" 1640#include "ccorb.h" 1641#include "ccsdsym.h" 1642C 1643 CALL AROUND('ENTERED TEMPPR') 1644 WRITE(LUPRI,*) 'IOPT = ',IOPT 1645C 1646 DO NA = 1,NVIRT 1647 DO NB = 1,NVIRT 1648 DO NC = 1,NVIRT 1649 DO NI = 1,NRHFT 1650 NAI = NVIRT*(NI-1) + NA 1651 DO NJ = 1,NRHFT 1652 NBJ = NVIRT*(NJ-1) + NB 1653 DO NK = 1,NRHFT 1654 NCK = NVIRT*(NK-1) + NC 1655C 1656 if (abs(t3am(nai,nbj,nck)) .gt. 1.0D-11) then 1657 if (iopt .eq. 1) then 1658 write(lupri,1) 't3am (^{',na,',',nb,',',nc,'}_{',ni,',', 1659 * nj,',',nk,'}) -> (', 1660 * NAI,',',NBJ,',',NCK,') = ',t3am(nai,nbj,nck) 1661 else if (iopt .eq. 2) then 1662 write(lupri,1) 't3bar(^{',na,',',nb,',',nc,'}_{',ni,',', 1663 * nj,',',nk,'}) -> (', 1664 * NAI,',',NBJ,',',NCK,') = ',t3am(nai,nbj,nck) 1665 else 1666 call quit('Wrong IOPT in temppr') 1667 endif 1668 endif 1669C 1670 ENDDO 1671 ENDDO 1672 ENDDO 1673 ENDDO 1674 ENDDO 1675 ENDDO 1676C 1677 1 FORMAT(1x,A8,I3,A1,I3,A1,I3,A3,I3,A1,I3,A1,I3,A7,I3,A1,I3,A1,I3, 1678 * A4,E20.10) 1679C 1680 RETURN 1681 END 1682C /* Deck ccsdt_nonrelterm */ 1683 SUBROUTINE CCSDT_NONREL_TERM(DENAB,DENIJ,T3AM,T3BAR) 1684C 1685C 1686C 1687#include "implicit.h" 1688 PARAMETER (TWO = 2.0D0, THREE = 3.0D0, HALF = 0.5D0) 1689 DIMENSION DENAB(NVIRT,NVIRT), DENIJ(NRHFT,NRHFT) 1690 DIMENSION T3AM(NT1AMX,NT1AMX,NT1AMX) 1691 DIMENSION T3BAR(NT1AMX,NT1AMX,NT1AMX) 1692#include "priunit.h" 1693COMMENT COMMENT COMMENT 1694C#include "inforb.h" 1695#include "ccorb.h" 1696COMMENT COMMENT COMMENT 1697#include "ccsdinp.h" 1698#include "ccsdsym.h" 1699C 1700 INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 1701C 1702 DO 100 I = 1,NRHFT 1703 DO 110 A = 1,NVIRT 1704 NAI = NVIRT*(I-1) + A 1705C 1706 DO 115 B = 1,NVIRT 1707 DO 120 J = 1,NRHFT 1708 NBJ = NVIRT*(J-1) + B 1709C 1710 DO 130 K = 1,NRHFT 1711 NAK = NVIRT*(K-1) + A 1712 NBK = NVIRT*(K-1) + B 1713 DO 140 C = 1,NVIRT 1714 NCK = NVIRT*(K-1) + C 1715 NCI = NVIRT*(I-1) + C 1716 NCJ = NVIRT*(J-1) + C 1717C 1718 DO 150 L = 1, NRHFT 1719C 1720 NCL = NVIRT*(L-1) + C 1721C 1722 DENIJ(K,L) = DENIJ(K,L) 1723 * - HALF*T3BAR(NAI,NBJ,NCL)* 1724 * T3AM(NAI,NBJ,NCK) 1725C 1726 150 CONTINUE 1727C 1728 DO 160 D = 1, NVIRT 1729C 1730 NDK = NVIRT*(K-1) + D 1731C 1732 DENAB(C,D) = DENAB(C,D) 1733 * + HALF*T3BAR(NAI,NBJ,NCK)* 1734 * T3AM(NAI,NBJ,NDK) 1735C 1736 160 CONTINUE 1737C 1738 140 CONTINUE 1739 130 CONTINUE 1740C 1741 120 CONTINUE 1742 115 CONTINUE 1743C 1744 110 CONTINUE 1745 100 CONTINUE 1746C 1747 RETURN 1748 END 1749 SUBROUTINE ECCSDT(T2AM,OMEGA2) 1750C 1751C 1752C 1753#include "implicit.h" 1754 PARAMETER (TWO=2.0d0) 1755 DIMENSION T2AM(NT1AMX,NT1AMX),OMEGA2(NT1AMX,NT1AMX) 1756#include "priunit.h" 1757COMMENT COMMENT COMMENT 1758C#include "inforb.h" 1759#include "ccorb.h" 1760COMMENT COMMENT COMMENT 1761#include "ccsdinp.h" 1762#include "ccsdsym.h" 1763C 1764 INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 1765C 1766 ENERGY = 0.0d0 1767 DO 100 I = 1,NRHFT 1768 DO 110 A = 1,NVIRT 1769 NAI = NVIRT*(I-1) + A 1770 DO 120 J = 1,NRHFT 1771 NAJ = NVIRT*(J-1) + A 1772 DO 130 B = 1,NVIRT 1773 NBJ = NVIRT*(J-1) + B 1774 NBI = NVIRT*(I-1) + B 1775 DELTA = OMEGA2(NAI,NBJ) 1776 * *(TWO*T2AM(NAI,NBJ) - T2AM(NAJ,NBI)) 1777 ENERGY = ENERGY + DELTA 1778 130 CONTINUE 1779 120 CONTINUE 1780 110 CONTINUE 1781 100 CONTINUE 1782C 1783 WRITE(LUPRI,*) 'Energy contribution E4(TT) :',ENERGY 1784C 1785 RETURN 1786 END 1787C /* Deck ccsdt_omega2 */ 1788 SUBROUTINE CCSDT_OMEGA2OLD(OMEGA2,XINT1T,XINT2T,FOCK,T3AM) 1789C 1790C 1791C 1792#include "implicit.h" 1793 PARAMETER (TWO = 2.0D0) 1794 DIMENSION XINT1T(NT1AMX,NVIRT,NVIRT) 1795 DIMENSION XINT2T(NT1AMX,NRHFT,NRHFT) 1796 DIMENSION OMEGA2(NT1AMX,NT1AMX) 1797 DIMENSION FOCK(NORBT,NORBT) 1798 DIMENSION T3AM(NT1AMX,NT1AMX,NT1AMX) 1799#include "priunit.h" 1800COMMENT COMMENT COMMENT 1801C#include "inforb.h" 1802#include "ccorb.h" 1803COMMENT COMMENT COMMENT 1804#include "ccsdinp.h" 1805#include "ccsdsym.h" 1806C 1807 INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 1808C 1809 DO 100 I = 1,NRHFT 1810 DO 110 A = 1,NVIRT 1811 NAI = NVIRT*(I-1) + A 1812C 1813 DO 120 J = 1,NRHFT 1814 DO 130 B = 1,NVIRT 1815 NBJ = NVIRT*(J-1) + B 1816C 1817 DO 140 K = 1,NRHFT 1818 NBK = NVIRT*(K-1) + B 1819 NAK = NVIRT*(K-1) + A 1820 DO 150 C = 1,NVIRT 1821C 1822 NCK = NVIRT*(K-1) + C 1823 NCJ = NVIRT*(J-1) + C 1824 NCI = NVIRT*(I-1) + C 1825C 1826 OMEGA2(NBJ,NAI) = OMEGA2(NBJ,NAI) - 1827 * (T3AM(NAI,NBJ,NCK)-T3AM(NAK,NBJ,NCI))*FOCK(K,NRHFT+C) 1828C 1829 DO 160 D = 1,NVIRT 1830C 1831 NDJ = NVIRT*(J-1) + D 1832 NDK = NVIRT*(K-1) + D 1833 NDI = NVIRT*(I-1) + D 1834C 1835 OMEGA2(NBJ,NAI) = OMEGA2(NBJ,NAI) + 1836 * (T3AM(NBK,NCI,NDJ) - TWO*T3AM(NBJ,NCI,NDK) + T3AM(NBJ,NCK,NDI)) 1837 * *XINT1T(NDK,A,C) 1838C 1839 160 CONTINUE 1840C 1841 DO 170 L = 1,NRHFT 1842C 1843 NBL = NVIRT*(L-1) + B 1844 NCL = NVIRT*(L-1) + C 1845 NAL = NVIRT*(L-1) + A 1846C 1847 OMEGA2(NBJ,NAI) = OMEGA2(NBJ,NAI) - 1848 * (T3AM(NBL,NAK,NCJ) - TWO*T3AM(NBJ,NAK,NCL) + T3AM(NBJ,NAL,NCK)) 1849 * *XINT2T(NCL,K,I) 1850C 1851 170 CONTINUE 1852C 1853 150 CONTINUE 1854 140 CONTINUE 1855C 1856 130 CONTINUE 1857 120 CONTINUE 1858C 1859 110 CONTINUE 1860 100 CONTINUE 1861C 1862 DO 200 NAI = 1,NT1AMX 1863 DO 210 NBJ = 1,NAI 1864C 1865 XAIBJ = OMEGA2(NAI,NBJ) + OMEGA2(NBJ,NAI) 1866 OMEGA2(NAI,NBJ) = XAIBJ 1867 OMEGA2(NBJ,NAI) = XAIBJ 1868C 1869 210 CONTINUE 1870 200 CONTINUE 1871C 1872 RETURN 1873 END 1874C /* Deck ccsdt_omega1old */ 1875 SUBROUTINE CCSDT_OMEGA1OLD(OMEGA1,XIAJB,T3AM) 1876C 1877C 1878C 1879#include "implicit.h" 1880 PARAMETER (TWO = 2.0D0) 1881 DIMENSION OMEGA1(NT1AMX),XIAJB(NT1AMX,NT1AMX) 1882 DIMENSION T3AM(NT1AMX,NT1AMX,NT1AMX) 1883#include "priunit.h" 1884COMMENT COMMENT COMMENT 1885C#include "inforb.h" 1886#include "ccorb.h" 1887COMMENT COMMENT COMMENT 1888#include "ccsdinp.h" 1889#include "ccsdsym.h" 1890C 1891 INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 1892C 1893 DO 100 I = 1,NRHFT 1894 DO 110 A = 1,NVIRT 1895 NAI = NVIRT*(I-1) + A 1896C 1897 DO 120 NBJ = 1,NT1AMX 1898C 1899 DO 130 K = 1,NRHFT 1900 NAK = NVIRT*(K-1) + A 1901 DO 140 C = 1,NVIRT 1902 NCK = NVIRT*(K-1) + C 1903 NCI = NVIRT*(I-1) + C 1904 OMEGA1(NAI) = OMEGA1(NAI) - XIAJB(NCK,NBJ)* 1905 * ( T3AM(NAI,NBJ,NCK) - T3AM(NAK,NBJ,NCI) ) 1906 140 CONTINUE 1907 130 CONTINUE 1908C 1909 120 CONTINUE 1910C 1911 110 CONTINUE 1912 100 CONTINUE 1913C 1914 RETURN 1915 END 1916C /* Deck cceta_omega2 */ 1917 SUBROUTINE CCETA_OMEGA2(OMEGA2,OMEGA22,T3AM,T3BAR,XINT1T, 1918 * XINT2T,FOCK) 1919C 1920C Kasper Hald, Fall 2001. 1921C 1922C Calculate the doubles of the right hand side of the t1-bar and t2-bar 1923C equation. 1924C NOTE : OMEGA22 will be erased. 1925C 1926#include "implicit.h" 1927C 1928 INTEGER INDEX, NAI, NBJ, NBK, NAK, NCK, NCJ, NCI, NDJ, NDK, NDI 1929 INTEGER NBL, NCL, NAL, NAJ, NBI, AI 1930C 1931#if defined (SYS_CRAY) 1932 REAL OMEGA2(NT1AMX,NT1AMX) 1933 REAL OMEGA22(NT1AMX,NT1AMX) 1934 REAL T3AM(NT1AMX,NT1AMX,NT1AMX) 1935 REAL T3BAR(NT1AMX,NT1AMX,NT1AMX) 1936 REAL XINT1T(NT1AMX,NVIRT,NVIRT) 1937 REAL XINT2T(NT1AMX,NRHFT,NRHFT) 1938 REAL FOCK(NORBT,NORBT) 1939 REAL ONE, TWO, XAIBJ 1940#else 1941 DOUBLE PRECISION OMEGA2(NT1AMX,NT1AMX) 1942 DOUBLE PRECISION OMEGA22(NT1AMX,NT1AMX) 1943 DOUBLE PRECISION T3AM(NT1AMX,NT1AMX,NT1AMX) 1944 DOUBLE PRECISION T3BAR(NT1AMX,NT1AMX,NT1AMX) 1945 DOUBLE PRECISION XINT1T(NT1AMX,NVIRT,NVIRT) 1946 DOUBLE PRECISION XINT2T(NT1AMX,NRHFT,NRHFT) 1947 DOUBLE PRECISION FOCK(NORBT,NORBT) 1948 DOUBLE PRECISION ONE, TWO, XAIBJ 1949#endif 1950 PARAMETER (ONE = 1.0D0, TWO = 2.0D0) 1951C 1952#include "priunit.h" 1953C#include "inforb.h" 1954#include "ccorb.h" 1955#include "ccsdinp.h" 1956#include "ccsdsym.h" 1957C 1958 INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 1959C 1960C DO AI = 1, NT1AMX 1961C DO A = 1, NVIRT 1962C DO B = 1, NVIRT 1963C if (abs(xint1t(ai,a,b)) .gt. 1.0D-9) then 1964C write(lupri,*) 'int1t(',ai,',',a,',',b,') = ', 1965C * xint1t(ai,a,b) 1966C endif 1967C ENDDO 1968C ENDDO 1969C ENDDO 1970C DO AI = 1, NT1AMX 1971C DO I = 1, NRHFT 1972C DO J = 1, NRHFT 1973C if (abs(xint2t(ai,i,j)) .gt. 1.0D-9) then 1974C write(lupri,*) 'int2t(',ai,',',i,',',j,') = ', 1975C * xint2t(ai,i,j) 1976C endif 1977C ENDDO 1978C ENDDO 1979C ENDDO 1980C 1981C--------------------------------------------- 1982C Start by calculating the T3AM terms 1983C and sum up. 1984C--------------------------------------------- 1985C 1986 DO I = 1,NRHFT 1987 DO A = 1,NVIRT 1988 NAI = NVIRT*(I-1) + A 1989C 1990 DO J = 1,NRHFT 1991 DO B = 1,NVIRT 1992 NBJ = NVIRT*(J-1) + B 1993C 1994 DO K = 1,NRHFT 1995 NBK = NVIRT*(K-1) + B 1996 NAK = NVIRT*(K-1) + A 1997 DO C = 1,NVIRT 1998C 1999 NCK = NVIRT*(K-1) + C 2000 NCJ = NVIRT*(J-1) + C 2001 NCI = NVIRT*(I-1) + C 2002C 2003 if (.true.) then 2004 OMEGA2(NBJ,NAI) = OMEGA2(NBJ,NAI) + 2005 * ( 2006 * T3AM(NBJ,NAI,NCK) 2007 * -T3AM(NBJ,NAK,NCI) 2008 * )*FOCK(K,NRHFT+C) 2009C 2010 endif 2011C 2012 DO D = 1,NVIRT 2013C 2014 NDJ = NVIRT*(J-1) + D 2015 NDK = NVIRT*(K-1) + D 2016 NDI = NVIRT*(I-1) + D 2017C 2018 if (.false.) then 2019 OMEGA2(NBJ,NAI) = OMEGA2(NBJ,NAI) - 2020 * (T3AM(NBK,NCI,NDJ) - TWO*T3AM(NBJ,NCI,NDK) + T3AM(NBJ,NCK,NDI)) 2021 * *XINT1T(NDK,A,C) 2022C 2023 endif 2024C 2025 ENDDO 2026C 2027 DO L = 1,NRHFT 2028C 2029 NBL = NVIRT*(L-1) + B 2030 NCL = NVIRT*(L-1) + C 2031 NAL = NVIRT*(L-1) + A 2032C 2033 if (.false.) then 2034 OMEGA2(NBJ,NAI) = OMEGA2(NBJ,NAI) + 2035 * (T3AM(NBL,NAK,NCJ) - TWO*T3AM(NBJ,NAK,NCL) + T3AM(NBJ,NAL,NCK)) 2036 * *XINT2T(NCL,K,I) 2037C 2038 endif 2039C 2040 ENDDO 2041C 2042 ENDDO 2043 ENDDO 2044C 2045 ENDDO 2046 ENDDO 2047C 2048 ENDDO 2049 ENDDO 2050C 2051C------------------------------------------------ 2052C Make P^{ab}_{ij} and P^{2CME}_{aibj} for 2053C the T3AM contributions. 2054C------------------------------------------------ 2055C 2056 DO NAI = 1,NT1AMX 2057 DO NBJ = 1,NAI 2058C 2059 XAIBJ = OMEGA2(NAI,NBJ) + OMEGA2(NBJ,NAI) 2060 OMEGA2(NAI,NBJ) = XAIBJ 2061 OMEGA2(NBJ,NAI) = XAIBJ 2062C 2063 ENDDO 2064 ENDDO 2065C 2066 CALL DZERO(OMEGA22,NT1AMX*NT1AMX) 2067C 2068 DO A = 1,NVIRT 2069 DO I = 1,NRHFT 2070 DO B = 1,NVIRT 2071 DO J = 1,NRHFT 2072C 2073 NAI = NVIRT*(I-1) + A 2074 NAJ = NVIRT*(J-1) + A 2075 NBI = NVIRT*(I-1) + B 2076 NBJ = NVIRT*(J-1) + B 2077 OMEGA22(NAI,NBJ) = OMEGA22(NAI,NBJ) + TWO*OMEGA2(NAI,NBJ) 2078 OMEGA22(NAJ,NBI) = OMEGA22(NAJ,NBI) - OMEGA2(NAI,NBJ) 2079C 2080 ENDDO 2081 ENDDO 2082 ENDDO 2083 ENDDO 2084C 2085C--------------------------------------------- 2086C Now calculate the T3BAR terms 2087C--------------------------------------------- 2088C 2089 CALL DZERO(OMEGA2,NT1AMX*NT1AMX) 2090C 2091 DO I = 1,NRHFT 2092 DO A = 1,NVIRT 2093 NAI = NVIRT*(I-1) + A 2094C 2095 DO J = 1,NRHFT 2096 DO B = 1,NVIRT 2097 NBJ = NVIRT*(J-1) + B 2098C 2099 DO K = 1,NRHFT 2100 NBK = NVIRT*(K-1) + B 2101 NAK = NVIRT*(K-1) + A 2102 DO C = 1,NVIRT 2103C 2104 NCK = NVIRT*(K-1) + C 2105 NCJ = NVIRT*(J-1) + C 2106 NCI = NVIRT*(I-1) + C 2107C 2108 DO D = 1,NVIRT 2109C 2110 NDJ = NVIRT*(J-1) + D 2111 NDK = NVIRT*(K-1) + D 2112 NDI = NVIRT*(I-1) + D 2113C 2114 if (.true.) then 2115 OMEGA2(NBJ,NAI) = OMEGA2(NBJ,NAI) + 2116 * T3BAR(NBJ,NCI,NDK)*XINT1T(NDK,A,C) 2117C 2118 endif 2119C 2120 ENDDO 2121C 2122 DO L = 1,NRHFT 2123C 2124 NBL = NVIRT*(L-1) + B 2125 NCL = NVIRT*(L-1) + C 2126 NAL = NVIRT*(L-1) + A 2127C 2128 if (.true.) then 2129 OMEGA2(NBJ,NAI) = OMEGA2(NBJ,NAI) + 2130 * - T3BAR(NBJ,NAK,NCL)*XINT2T(NCL,K,I) 2131C 2132 endif 2133C 2134 ENDDO 2135C 2136 ENDDO 2137 ENDDO 2138C 2139 ENDDO 2140 ENDDO 2141C 2142 ENDDO 2143 ENDDO 2144C 2145C---------------------------------------------------- 2146C Make P^{ab}_{ij} for the T3BAR contributions. 2147C---------------------------------------------------- 2148C 2149 DO NAI = 1,NT1AMX 2150 DO NBJ = 1,NAI 2151C 2152 XAIBJ = OMEGA2(NAI,NBJ) + OMEGA2(NBJ,NAI) 2153 OMEGA2(NAI,NBJ) = XAIBJ 2154 OMEGA2(NBJ,NAI) = XAIBJ 2155C 2156 ENDDO 2157 ENDDO 2158C 2159C------------------------------- 2160C Sum to final result. 2161C------------------------------- 2162C 2163 CALL DAXPY(NT1AMX*NT1AMX,ONE,OMEGA22,1,OMEGA2,1) 2164C 2165 RETURN 2166 END 2167C /* Deck ccsdt_kappa */ 2168 SUBROUTINE CCSDT_KAPPA(T3AM,T3BAR,T1AM,T2AM,VIRT,OCC, 2169 * OMEGA2,OMEGA22) 2170C 2171C Kasper Hald, Fall 2001. 2172C 2173C Calculate the doubles of the right hand side of the Kappa 2174C equation. 2175C 2176#include "implicit.h" 2177C 2178 INTEGER INDEX, NAI, NBJ, NBK, NAK, NCK, NCJ, NCI, NDJ, NDK, NDI 2179 INTEGER NBL, NCL, NAL, NAJ, NBI, AI 2180C 2181#if defined (SYS_CRAY) 2182 REAL T3AM(NT1AMX,NT1AMX,NT1AMX) 2183 REAL T3BAR(NT1AMX,NT1AMX,NT1AMX) 2184 REAL T1AM(NT1AMX) 2185 REAL T2AM(NT1AMX,NT1AMX) 2186 REAL VIRT(NVIRT,NVIRT,NRHFT,NVIRT) 2187 REAL OCC(NRHFT,NVIRT,NRHFT,NRHFT) 2188 REAL OMEGA2(NT1AMX,NT1AMX) 2189 REAL OMEGA22(NT1AMX,NT1AMX) 2190 REAL ONE, TWO, XAIBJ 2191#else 2192 DOUBLE PRECISION T3AM(NT1AMX,NT1AMX,NT1AMX) 2193 DOUBLE PRECISION T3BAR(NT1AMX,NT1AMX,NT1AMX) 2194 DOUBLE PRECISION T1AM(NT1AMX) 2195 DOUBLE PRECISION T2AM(NT1AMX,NT1AMX) 2196 DOUBLE PRECISION VIRT(NVIRT,NVIRT,NRHFT,NVIRT) 2197 DOUBLE PRECISION OCC(NRHFT,NVIRT,NRHFT,NRHFT) 2198 DOUBLE PRECISION OMEGA2(NT1AMX,NT1AMX) 2199 DOUBLE PRECISION OMEGA22(NT1AMX,NT1AMX) 2200 DOUBLE PRECISION ONE, TWO, XAIBJ 2201#endif 2202C 2203 LOGICAL L3ABIC, L3IAJB, L3IAJK1, L3IAJK2, L3IAJK3 2204 LOGICAL LBARABCI, LBARAIJK, LDEBUG 2205 PARAMETER (ONE = 1.0D0, TWO = 2.0D0) 2206C 2207#include "priunit.h" 2208C#include "inforb.h" 2209#include "ccorb.h" 2210#include "ccsdinp.h" 2211#include "ccsdsym.h" 2212C 2213 INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 2214C 2215C------------------------------------------------------- 2216C Specify which terms that should be calculated. 2217C------------------------------------------------------- 2218C 2219 LDEBUG = .false. 2220 L3ABIC = .true. 2221 L3IAJB = .true. 2222 L3IAJK1 = .true. 2223 L3IAJK2 = .true. 2224 L3IAJK3 = .true. 2225 LBARABCI = .true. 2226 LBARAIJK = .true. 2227C 2228C--------------------------------------------- 2229C Start by calculating the T3AM terms 2230C and sum up. 2231C--------------------------------------------- 2232C 2233 DO I = 1,NRHFT 2234 DO A = 1,NVIRT 2235 NAI = NVIRT*(I-1) + A 2236C 2237 DO J = 1,NRHFT 2238 DO B = 1,NVIRT 2239 NBI = NVIRT*(I-1) + B 2240 NBJ = NVIRT*(J-1) + B 2241 NAJ = NVIRT*(J-1) + A 2242C 2243 DO K = 1,NRHFT 2244 NBK = NVIRT*(K-1) + B 2245 NAK = NVIRT*(K-1) + A 2246 DO C = 1,NVIRT 2247C 2248 NCK = NVIRT*(K-1) + C 2249 NCJ = NVIRT*(J-1) + C 2250 NCI = NVIRT*(I-1) + C 2251C 2252 if (L3IAJB) then 2253 OMEGA2(NBJ,NAI) = OMEGA2(NBJ,NAI) + 2254 * ( T3AM(NCK,NAI,NBJ) -T3AM(NCI,NAK,NBJ) )*T1AM(NCK) 2255C 2256COMMENT COMMENT 2257COMMENT COMMENT 2258C if (abs(T3AM(NCI,NAK,NBJ)*T1AM(NCK)) .gt. 1.0D-11) then 2259C write(lupri,'(A,6I3)') 'Cont. from A,B,C,I,J,K : ',a,b,c,i,j,k 2260C write(lupri,*) 'Term is equal to ', 2261C * -TWO*T3AM(NCI,NAK,NBJ)*T1AM(NCK) 2262C write(lupri,*) 'With T1 and T3 : ',T1AM(NCK),T3AM(NCI,NAK,NBJ) 2263C endif 2264COMMENT COMMENT 2265COMMENT COMMENT 2266 endif 2267C 2268 DO D = 1,NVIRT 2269C 2270 NDJ = NVIRT*(J-1) + D 2271 NDK = NVIRT*(K-1) + D 2272 NDI = NVIRT*(I-1) + D 2273C 2274 if (L3ABIC) then 2275 VIRT(A,B,I,C) = VIRT(A,B,I,C) - TWO* 2276 * ( 2277 * T3AM(NDK,NCJ,NBI) 2278 * -TWO*T3AM(NBJ,NCI,NDK) 2279 * + T3AM(NBJ,NCK,NDI) 2280 * ) 2281 * *(TWO*T2AM(NDK,NAJ)-T2AM(NDJ,NAK)) 2282C 2283COMMENT COMMENT 2284COMMENT COMMENT 2285C if (abs(TWO* 2286C * (T3AM(NDK,NCJ,NBI)-TWO*T3AM(NBJ,NCI,NDK)+T3AM(NBJ,NCK,NDI)) 2287C * (-TWO*T3AM(NBJ,NCI,NDK)) 2288C * *(TWO*T2AM(NDK,NAJ)-T2AM(NDJ,NAK)) 2289C * ) .gt. 1.0D-12) then 2290C write(lupri,*) 'ABIC contribution : ',TWO* 2291C * (T3AM(NDK,NCJ,NBI)-TWO*T3AM(NBJ,NCI,NDK)+T3AM(NBJ,NCK,NDI)) 2292C * *(TWO*T2AM(NDK,NAJ)-T2AM(NDJ,NAK)) 2293C write(lupri,*) 'T3_1 = ',T3AM(NDK,NCJ,NBI) 2294C write(lupri,*) 'T3_2 = ',-TWO*T3AM(NBJ,NCI,NDK) 2295C write(lupri,*) 'T3_3 = ',T3AM(NBJ,NCK,NDI) 2296C write(lupri,*) '2CME OF T2 = ', 2297C * (TWO*T2AM(NDK,NAJ)-T2AM(NDJ,NAK)) 2298C write(lupri,*) ' ' 2299C endif 2300COMMENT COMMENT 2301COMMENT COMMENT 2302C 2303 endif 2304C 2305 ENDDO 2306C 2307 DO L = 1,NRHFT 2308C 2309 NAL = NVIRT*(L-1) + A 2310 NBL = NVIRT*(L-1) + B 2311 NCL = NVIRT*(L-1) + C 2312C 2313 if (L3IAJK1) then 2314 OCC(I,A,J,K) = OCC(I,A,J,K) + TWO* 2315 * ( 2316 * + T3AM(NBJ,NAL,NCI) 2317 * - TWO*T3AM(NBJ,NAI,NCL) 2318 * + T3AM(NBI,NAJ,NCL) 2319 * ) 2320 * *(TWO*T2AM(NCL,NBK)-T2AM(NCK,NBL)) 2321C 2322 endif 2323C 2324 if (L3IAJK2) then 2325 OCC(I,A,J,J) = OCC(I,A,J,J) + TWO*TWO* 2326 * (T3AM(NCL,NBK,NAI)-T3AM(NCL,NBI,NAK))* 2327 * (TWO*T2AM(NCL,NBK)-T2AM(NCK,NBL)) 2328C 2329 endif 2330C 2331 if (L3IAJK3) then 2332 OCC(I,A,J,I) = OCC(I,A,J,I) - TWO* 2333 * (T3AM(NCL,NBK,NAJ)-T3AM(NCL,NBJ,NAK))* 2334 * (TWO*T2AM(NCL,NBK)-T2AM(NCK,NBL)) 2335C 2336 endif 2337C 2338 ENDDO 2339C 2340 ENDDO 2341 ENDDO 2342C 2343 ENDDO 2344 ENDDO 2345C 2346 ENDDO 2347 ENDDO 2348C 2349C---------------------------------------------------------------------- 2350C Symmetrize the omega2 density with both P^{2CME} and normal P 2351C---------------------------------------------------------------------- 2352C 2353 IF (L3IAJB) THEN 2354 DO NAI = 1, NT1AMX 2355 DO NBJ = 1, NAI 2356 XAIBJ = OMEGA2(NAI,NBJ) + OMEGA2(NBJ,NAI) 2357 OMEGA2(NAI,NBJ) = XAIBJ 2358 OMEGA2(NBJ,NAI) = XAIBJ 2359 ENDDO 2360 ENDDO 2361C 2362 CALL DZERO(OMEGA22,NT1AMX*NT1AMX) 2363C 2364 DO A = 1,NVIRT 2365 DO I = 1,NRHFT 2366 DO B = 1,NVIRT 2367 DO J = 1,NRHFT 2368C 2369 NAI = NVIRT*(I-1) + A 2370 NAJ = NVIRT*(J-1) + A 2371 NBI = NVIRT*(I-1) + B 2372 NBJ = NVIRT*(J-1) + B 2373 OMEGA22(NAI,NBJ) = OMEGA22(NAI,NBJ) + TWO*OMEGA2(NAI,NBJ) 2374 OMEGA22(NAJ,NBI) = OMEGA22(NAJ,NBI) - OMEGA2(NAI,NBJ) 2375C 2376 ENDDO 2377 ENDDO 2378 ENDDO 2379 ENDDO 2380 ENDIF 2381C 2382C-------------------------------------------- 2383C Print the calculated terms : 2384C-------------------------------------------- 2385C 2386 if (L3ABIC .AND. LDEBUG) then 2387 do a = 1, nvirt 2388 do b = 1, nvirt 2389 do i = 1, nrhft 2390 do c = 1, nvirt 2391 if (abs(virt(a,b,i,c)) .gt. 1.0D-9) then 2392 write(lupri,'(A,I3,A,I3,A,I3,A,I3,A,e20.10)') 2393 * 'd2-3-abic(',a,',',b,',',i,',',c,') = ', 2394 * virt(a,b,i,c) 2395 endif 2396 enddo 2397 enddo 2398 enddo 2399 enddo 2400 write(lupri,*) ' ' 2401 endif 2402C 2403 if (L3IAJB .AND. LDEBUG) then 2404 do a = 1, nvirt 2405 do i = 1, nrhft 2406 NAI = NVIRT*(I-1) + A 2407 do b = 1, nvirt 2408 do j = 1, nrhft 2409 NBJ = NVIRT*(J-1) + B 2410 if (abs(omega22(nai,nbj)) .gt. 1.0D-11) then 2411 write(lupri,'(A,I3,A,I3,A,I3,A,I3,A,e20.10)') 2412 * 'd2-3-iajb(',i,',',a,',',j,',',b,') = ', 2413 * omega22(nai,nbj) 2414 endif 2415 enddo 2416 enddo 2417 enddo 2418 enddo 2419 write(lupri,*) ' ' 2420 call around('Two electron density from NODDY') 2421 call output(omega22,1,nt1amx,1,nt1amx,nt1amx,nt1amx,1,lupri) 2422 write(lupri,*) ' ' 2423 endif 2424 2425C 2426 if ((L3IAJK1 .OR. L3IAJK2 .OR. L3IAJK3) .AND. LDEBUG) then 2427 do i = 1, nrhft 2428 do a = 1, nvirt 2429 do j = 1, nrhft 2430 do k = 1, nrhft 2431 if (abs(occ(i,a,j,k)) .gt. 1.0D-9) then 2432 write(lupri,'(A,I3,A,I3,A,I3,A,I3,A,e20.10)') 2433 * 'd2-3-iajk(',i,',',a,',',j,',',k,') = ', 2434 * occ(i,a,j,k) 2435 endif 2436 enddo 2437 enddo 2438 enddo 2439 enddo 2440 write(lupri,*) ' ' 2441 endif 2442C 2443C--------------------------------------------- 2444C Now calculate the T3BAR terms 2445C--------------------------------------------- 2446C 2447 CALL DZERO(VIRT,NVIRT*NVIRT*NRHFT*NVIRT) 2448 CALL DZERO(OCC,NRHFT*NVIRT*NRHFT*NRHFT) 2449C 2450 DO I = 1,NRHFT 2451 DO A = 1,NVIRT 2452 NAI = NVIRT*(I-1) + A 2453C 2454 DO J = 1,NRHFT 2455 DO B = 1,NVIRT 2456 NAJ = NVIRT*(J-1) + A 2457 NBJ = NVIRT*(J-1) + B 2458C 2459 DO K = 1,NRHFT 2460 NBK = NVIRT*(K-1) + B 2461 NAK = NVIRT*(K-1) + A 2462 DO C = 1,NVIRT 2463C 2464 NCK = NVIRT*(K-1) + C 2465 NCJ = NVIRT*(J-1) + C 2466 NCI = NVIRT*(I-1) + C 2467C 2468 DO D = 1,NVIRT 2469C 2470 NDJ = NVIRT*(J-1) + D 2471 NDK = NVIRT*(K-1) + D 2472 NDI = NVIRT*(I-1) + D 2473C 2474 if (LBARABCI) then 2475C 2476C Store it as a,b,i,c though it is a,b,c,i 2477C 2478 VIRT(A,B,I,C) = VIRT(A,B,I,C) - 2479 * T3BAR(NAJ,NCI,NDK)*T2AM(NDK,NBJ) 2480COMMENT COMMENT 2481COMMENT COMMENT 2482C if (abs(T3BAR(NAJ,NCI,NDK)*T2AM(NDK,NBJ)) .gt. 1.0D-12) then 2483C write(lupri,*) 'ABCI-bar contribution : ', 2484C * T3BAR(NAJ,NCI,NDK)*T2AM(NDK,NBJ) 2485C write(lupri,*) 'T3, T2 : ',T3BAR(NAJ,NCI,NDK),T2AM(NDK,NBJ) 2486C endif 2487COMMENT COMMENT 2488COMMENT COMMENT 2489C 2490 endif 2491C 2492 ENDDO 2493C 2494 DO L = 1,NRHFT 2495C 2496 NBL = NVIRT*(L-1) + B 2497 NCL = NVIRT*(L-1) + C 2498 NAL = NVIRT*(L-1) + A 2499C 2500 if (LBARAIJK) then 2501C 2502C Store as though it is i,a,j,k though it is a,i,j,k 2503C 2504 OCC(I,A,J,K) = OCC(I,A,J,K) - 2505 * T3BAR(NAI,NBK,NCL)*T2AM(NCL,NBJ) 2506C 2507 endif 2508C 2509 ENDDO 2510C 2511 ENDDO 2512 ENDDO 2513C 2514 ENDDO 2515 ENDDO 2516C 2517 ENDDO 2518 ENDDO 2519C 2520C---------------------------------------------------- 2521C Print the contributions from t3bar. 2522C---------------------------------------------------- 2523C 2524 if (LBARABCI .AND. LDEBUG) then 2525 do a = 1, nvirt 2526 do b = 1, nvirt 2527 do i = 1, nrhft 2528 do c = 1, nvirt 2529 if (abs(virt(a,b,i,c)) .gt. 1.0D-9) then 2530 write(lupri,'(A,I3,A,I3,A,I3,A,I3,A,e20.10)') 2531 * 'd2-bar-abci(',a,',',b,',',c,',',i,') = ', 2532 * virt(a,b,i,c) 2533 endif 2534 enddo 2535 enddo 2536 enddo 2537 enddo 2538 write(lupri,*) ' ' 2539 endif 2540C 2541 if (LBARAIJK .AND. LDEBUG) then 2542 do a = 1, nvirt 2543 do i = 1, nrhft 2544 do j = 1, nrhft 2545 do k = 1, nrhft 2546 if (abs(occ(i,a,j,k)) .gt. 1.0D-9) then 2547 write(lupri,'(A,I3,A,I3,A,I3,A,I3,A,e20.10)') 2548 * 'd2-bar-aijk(',a,',',i,',',j,',',k,') = ', 2549 * occ(i,a,j,k) 2550 endif 2551 enddo 2552 enddo 2553 enddo 2554 enddo 2555 write(lupri,*) ' ' 2556 endif 2557C 2558C------------------------------- 2559C Sum to final result. 2560C------------------------------- 2561C 2562 CALL DAXPY(NT1AMX*NT1AMX,ONE,OMEGA22,1,OMEGA2,1) 2563C 2564 RETURN 2565 END 2566C /* Deck debug_kappadiag */ 2567 SUBROUTINE DEBUG_KAPPADIAG(XKAPAA,XKAPII,T3AM,T3BAR) 2568C 2569C 2570#include "implicit.h" 2571#include "priunit.h" 2572#include "ccsdsym.h" 2573#include "ccorb.h" 2574C 2575 DIMENSION XKAPAA(NVIRT), XKAPII(NRHFT) 2576 DIMENSION T3AM(NT1AMX,NT1AMX,NT1AMX) 2577 DIMENSION T3BAR(NT1AMX,NT1AMX,NT1AMX) 2578C 2579 PARAMETER (HALF = 0.5D0) 2580C 2581C---------------------------------------------- 2582C Kappa_{aa} 2583C---------------------------------------------- 2584C 2585 DO A = 1, NVIRT 2586 DO I = 1, NRHFT 2587 NAI = NVIRT*(I-1) + A 2588 DO NBJ = 1, NT1AMX 2589 DO NCK = 1, NT1AMX 2590 XKAPAA(A) = XKAPAA(A) 2591 * + T3AM(NAI,NBJ,NCK)*T3BAR(NAI,NBJ,NCK) 2592 ENDDO 2593 ENDDO 2594 ENDDO 2595 ENDDO 2596C 2597C---------------------------------------------- 2598C Kappa_{ii} 2599C---------------------------------------------- 2600C 2601 DO I = 1, NRHFT 2602 DO A = 1, NVIRT 2603 NAI = NVIRT*(I-1) + A 2604 DO NBJ = 1, NT1AMX 2605 DO NCK = 1, NT1AMX 2606 XKAPII(I) = XKAPII(I) 2607 * + T3AM(NAI,NBJ,NCK)*T3BAR(NAI,NBJ,NCK) 2608 ENDDO 2609 ENDDO 2610 ENDDO 2611 ENDDO 2612C 2613C------------------------------ 2614C Scale the kappas 2615C------------------------------ 2616C 2617 CALL DSCAL(NVIRT,HALF,XKAPAA,1) 2618 CALL DSCAL(NRHFT,HALF,XKAPII,1) 2619C 2620 RETURN 2621 END 2622C /* Deck print_integrals */ 2623 SUBROUTINE PRINT_INTEGRALS(XINT1,XINT2) 2624C 2625#include "implicit.h" 2626C 2627 DIMENSION XINT1(NT1AMX,NVIRT,NVIRT) 2628 DIMENSION XINT2(NT1AMX,NRHFT,NRHFT) 2629 LOGICAL LDEBUG 2630C 2631#include "priunit.h" 2632#include "ccorb.h" 2633#include "ccsdinp.h" 2634#include "ccsdsym.h" 2635C 2636 LDEBUG = .FALSE. 2637C 2638 IF (LDEBUG) THEN 2639 DO NAI = 1, NT1AMX 2640 DO NB = 1, NVIRT 2641 DO NC = 1, NVIRT 2642 IF (ABS(XINT1(NAI,NB,NC)) .GT. 1.0D-9) THEN 2643 WRITE(LUPRI,'(A,3I3,A,E20.10)') 2644 * 'INT(NAI,B,C) = INT(',NAI,NB,NC,') = ', 2645 * XINT1(NAI,NB,NC) 2646 ENDIF 2647 ENDDO 2648 ENDDO 2649 ENDDO 2650 WRITE(LUPRI,*) ' ' 2651 DO NAI = 1, NT1AMX 2652 DO NL = 1, NVIRT 2653 DO NK = 1, NVIRT 2654 IF (ABS(XINT2(NAI,NL,NK)) .GT. 1.0D-9) THEN 2655 WRITE(LUPRI,'(A,3I3,A,E20.10)') 2656 * 'INT(NAI,K,L) = INT(',NAI,NL,NK,') = ', 2657 * XINT2(NAI,NL,NK) 2658 ENDIF 2659 ENDDO 2660 ENDDO 2661 ENDDO 2662 ENDIF 2663C 2664 RETURN 2665 END 2666C 2667C /* Deck ccsdt_unre */ 2668 SUBROUTINE CCSDT_UNRE(T2AM,T3AM,T3BAR,DENS,FCKFIELD) 2669C 2670#include "implicit.h" 2671#include "priunit.h" 2672#include "ccorb.h" 2673#include "ccsdsym.h" 2674C 2675 DIMENSION T2AM(NT1AMX,NT1AMX) 2676 DIMENSION T3AM(NT1AMX,NT1AMX,NT1AMX) 2677 DIMENSION T3BAR(NT1AMX,NT1AMX,NT1AMX) 2678 DIMENSION DENS(NT1AMX) 2679 DIMENSION FCKFIELD(NORBT,NORBT) 2680COMMENT COMMENT 2681 dimension sumtemp(nt1amx,nrhft,nrhft) 2682 dimension omeg1(nvirt,nrhft) 2683COMMENT COMMENT 2684 PARAMETER (HALF = 0.5D0) 2685C 2686COMMENT COMMENT 2687 call dzero(sumtemp,nt1amx*nrhft*nrhft) 2688 call dzero(omeg1,nvirt*nrhft) 2689COMMENT COMMENT 2690 SUM = 0.0D0 2691 SUM1 = 0.0D0 2692 SUM2 = 0.0D0 2693 SUM3 = 0.0D0 2694C 2695 DO NA = 1, NVIRT 2696 DO NI = 1, NRHFT 2697 NAI = NVIRT*(NI-1) + NA 2698 SUM = SUM + DENS(NAI)*FCKFIELD(NI,NRHFT+NA) 2699 ENDDO 2700 ENDDO 2701C 2702 DO NCK = 1,NT1AMX 2703C 2704 DO J = 1,NRHFT 2705 DO B = 1,NVIRT 2706C 2707 NBJ = NVIRT*(J-1) + B 2708C 2709 DO NAI = 1,NT1AMX 2710C 2711 DO D = 1, NVIRT 2712 NDJ = NVIRT*(J-1) + D 2713 DO L = 1, NRHFT 2714 NBL = NVIRT*(L-1) + B 2715 SUM1 = SUM1 2716 * - (T2AM(NAI,NBL) 2717 * *T2AM(NCK,NDJ) 2718 * *FCKFIELD(L,NRHFT+D))*T3BAR(NAI,NBJ,NCK) 2719 ENDDO 2720 ENDDO 2721C 2722 DO D = 1, NVIRT 2723 NDJ = NVIRT*(J-1) + D 2724 SUM2 = SUM2 2725 * + (HALF*T3AM(NAI,NDJ,NCK) 2726 * *FCKFIELD(NRHFT+B,NRHFT+D) 2727 * )*T3BAR(NAI,NBJ,NCK) 2728 ENDDO 2729C 2730 DO L = 1, NRHFT 2731 NBL = NVIRT*(L-1) + B 2732 SUM3 = SUM3 2733 * - (HALF*T3AM(NAI,NBL,NCK) 2734 * *FCKFIELD(L,J))*T3BAR(NAI,NBJ,NCK) 2735 ENDDO 2736 ENDDO 2737 ENDDO 2738 ENDDO 2739 ENDDO 2740C 2741 CALL AROUND('Contribution to the unrelaxed density :') 2742 write(lupri,*) 't^{*}*H(1)T3 = ',sum 2743 write(lupri,*) 't3bar*H(1)*T2*T2 = ',sum1 2744 write(lupri,*) 't3bar*H(1)*T3 (virt) = ',sum2 2745 write(lupri,*) 't3bar*H(1)*T3 (occ) = ',sum3 2746 write(lupri,*) 'Total contribution = ',sum+sum1+sum2+sum3 2747C 2748 RETURN 2749 END 2750C -- end of cc_fop3.F -- 2751