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 cc3_triplet */ 20 SUBROUTINE CC3_TRIPLET(ECURR,ISYMTR,OMEGA1,OMEGA2P,T2AM, 21 * C2AMP,C2AMM,FOCK, 22 * XLAMDP,XLAMDH,XLAMPC,XLAMHC,WORK,LWORK, 23 * LUFR2,FRHO2,IV,ITR,C1AM) 24! 25! Written by Kasper Hald 26! 27! Purpose: Calculate the CC3 corrections to the CCSD 28! triplet excitation energies. 29! 30! NB: The T2 amplitudes are assumed to be a full square, 31! C2AMP is (C2+) and C2AMM is (C2-). 32! The C2AMP and C2AMM are also assumed to be full squares 33! 34! NB: It is assumed that the vectors are allocated in the following 35! order: 36! OMEGA1(*), OMEGA2(*), T2AM(*), SCR(*), WRK(*). 37! 38! NB: This is ONLY a noddy code ... no intermediates etc. 39! 40 IMPLICIT NONE 41! 42#include "priunit.h" 43#include "dummy.h" 44#include "maxorb.h" 45#include "maxash.h" 46#include "mxcent.h" 47#include "aovec.h" 48#include "iratdef.h" 49#include "eritap.h" 50#include "ccorb.h" 51#include "inftap.h" 52#include "ccisao.h" 53#include "r12int.h" 54! 55 INTEGER INDEX, LWORK, KCMO, KSCR1, KFOCKD, KINT1T 56 INTEGER KINT2T, KINT3T, KINT4T, KINT1S, KINT2S, KINT3S, KINT4S 57 INTEGER KINT5S, KINT6S, KT3AM, KOME1, KOME2P, KOME2M, KIAJB 58 INTEGER KYIAJB, KEND1, LWRK1, NTOSYM 59 INTEGER ISYMD1, NTOT, ILLL, NUMDIS, IDEL2, IDEL, ISYMD 60 INTEGER ISYDIS, KXINT, KEND2, LWRK2, LUFR2, KOME2, IV, ITR 61 INTEGER IJ, NIJ, KRECNR, ISYMTR 62 INTEGER INDEXA(MXCORB_CC) 63! 64 INTEGER KOME22, IERRR1, IERRR2 65 INTEGER LUCC3R1, LUCC3R2, TEMP 66 INTEGER KINTTOT, KINTTOT2 67 INTEGER KOFF1, KOFF2 68 INTEGER KINTTMP1, KINTTMP2, KTRTMP1 69 INTEGER KTRTMP2, KTRTMP3, KTRTMP4 70! 71 INTEGER IDUM 72! 73#if defined (SYS_CRAY) 74 REAL ZERO, HALF, ONE, TWO, DTIME ! , TIMHER1, TIMHER2 75 REAL TIMRDAO, OMEGA1(*), OMEGA2P(*), T2AM(*) 76 REAL FOCK(*), XLAMDP(*), XLAMDH(*), XLAMPC(*) 77 REAL XLAMHC(*), WORK(LWORK), C2AMP(*), C2AMM(*) 78 REAL DDOT, RHO2N 79 REAL XONE, XHALF, ECURR 80 REAL C1AM(*) 81#else 82 DOUBLE PRECISION ZERO, HALF, ONE, TWO, DTIME ! , TIMHER1, TIMHER2 83 DOUBLE PRECISION TIMRDAO, OMEGA1(*), OMEGA2P(*), T2AM(*) 84 DOUBLE PRECISION FOCK(*), XLAMDP(*), XLAMDH(*), XLAMPC(*) 85 DOUBLE PRECISION XLAMHC(*), WORK(LWORK), C2AMP(*), C2AMM(*) 86 DOUBLE PRECISION DDOT, RHO2N 87 DOUBLE PRECISION XONE, XHALF, ECURR 88 DOUBLE PRECISION C1AM(*) 89#endif 90! 91 LOGICAL LOCDBG 92! 93 CHARACTER*8 FRHO2 94 CHARACTER*8 RHO1FIL, RHO2FIL 95! 96 PARAMETER (ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, TWO = 2.0D0) 97 PARAMETER (XONE = -1.0D00, XHALF = -0.5D0) 98 PARAMETER (LOCDBG = .FALSE.) 99! 100! 101#include "infind.h" 102#include "blocks.h" 103#include "ccsdinp.h" 104#include "ccsdsym.h" 105#include "ccsdio.h" 106#include "second.h" 107! 108 INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 109! 110 IF (NSYM .NE. 1) THEN 111 WRITE(LUPRI,*)'CC3 Triplet excitation energies using the' 112 WRITE(LUPRI,*)'noddy code can only be performed for NSYM=1' 113 CALL QUIT('No symmetry in this part yet') 114 ENDIF 115! 116 IF (LOCDBG) THEN 117 WRITE(LUPRI,*)'ENTERING THE CC3 FILE.' 118 WRITE(LUPRI,*)'The energy (ECURR) in this iteration is ',ECURR 119 RHO2N = DDOT(NT2SQ(ISYMTR),C2AMP(1),1,C2AMP(1),1) 120 WRITE(LUPRI,*) 121 * 'The norm of the C2(+) vector in CC3 :',RHO2N 122 RHO2N = DDOT(NT2SQ(ISYMTR),C2AMM(1),1,C2AMM(1),1) 123 WRITE(LUPRI,*) 124 * 'The norm of the C2(-) vector in CC3 :',RHO2N 125 RHO2N = DDOT(NT2SQ(1),T2AM(1),1,T2AM(1),1) 126 WRITE(LUPRI,*) 127 * 'The norm of the T2 vector in CC3 :',RHO2N 128 WRITE(LUPRI,*)' ' 129 CALL FLSHFO(LUPRI) 130 ENDIF 131! 132!------------------------ 133! Dynamic Allocation. 134!------------------------ 135! 136 KCMO = 1 137 KSCR1 = KCMO + NLAMDS 138 KFOCKD = KSCR1 + NT1AMX 139 KINT1T = KFOCKD + NORBTS 140 KINT2T = KINT1T + NT1AMX*NVIRT*NVIRT 141 KINT3T = KINT2T + NRHFT*NRHFT*NT1AMX 142 KINT4T = KINT3T + NT1AMX*NVIRT*NVIRT 143 KINT1S = KINT4T + NRHFT*NRHFT*NT1AMX 144 KINT2S = KINT1S + NT1AMX*NVIRT*NVIRT 145 KINT3S = KINT2S + NRHFT*NRHFT*NT1AMX 146 KINT4S = KINT3S + NT1AMX*NVIRT*NVIRT 147 KINT5S = KINT4S + NRHFT*NRHFT*NT1AMX 148 KINT6S = KINT5S + NT1AMX*NVIRT*NVIRT 149 KT3AM = KINT6S + NRHFT*NRHFT*NT1AMX 150 KOME1 = KT3AM + NT1AMX*NT1AMX*NT1AMX 151 KOME2P = KOME1 + NT1AMX 152 KOME2M = KOME2P + NT1AMX*NT1AMX 153 KIAJB = KOME2M + NT1AMX*NT1AMX 154 KYIAJB = KIAJB + NT1AMX*NT1AMX 155 KINTTOT= KYIAJB + NT1AMX*NT1AMX 156 KINTTOT2= KINTTOT+ NORBT*NORBT*NORBT*NORBT 157! 158! This is intermediates needed for the integrals. 159! Not needed at the same time 160! 161 KINTTMP1= KINTTOT2 + NORBT*NORBT*NORBT*NORBT 162 KINTTMP2= KINTTMP1 + NBAST*NBAST*NBAST*NBAST 163! 164 KTRTMP1= KINTTOT2 + NORBT*NORBT*NORBT*NORBT 165 KTRTMP2= KTRTMP1 + NRHFT*NORBT*NORBT*NORBT 166 KTRTMP3= KTRTMP2 + NVIRT*NORBT*NORBT*NORBT 167 KTRTMP4= KTRTMP3 + NRHFT*NORBT*NORBT*NORBT 168! 169 IF (NBAST*NBAST*NBAST*NBAST .GT. 170 * (NRHFT+NVIRT)*NORBT*NORBT*NORBT) THEN 171! 172 KEND1 = KINTTMP2 + NBAST*NBAST*NBAST*NBAST 173 ELSE 174 KEND1 = KTRTMP4 + NVIRT*NORBT*NORBT*NORBT 175 ENDIF 176! 177 LWRK1 = LWORK - KEND1 178! 179!------------------------------------------------------ 180! We sum the minus up in KOME2, 181! which is put at the same place as KIAJB. 182!------------------------------------------------------ 183! 184 KOME2 = KIAJB 185! 186 IF (LWRK1 .LT. 0) THEN 187 CALL QUIT('Insufficient space in CC3_TRIPLET') 188 ENDIF 189! 190!-------------------------------- 191! Initialize integral arrays. 192!-------------------------------- 193! 194 CALL DZERO(WORK(KINT1T),NT1AMX*NVIRT*NVIRT) 195 CALL DZERO(WORK(KINT2T),NT1AMX*NRHFT*NRHFT) 196 CALL DZERO(WORK(KINT3T),NT1AMX*NVIRT*NVIRT) 197 CALL DZERO(WORK(KINT4T),NT1AMX*NRHFT*NRHFT) 198 CALL DZERO(WORK(KINT1S),NT1AMX*NVIRT*NVIRT) 199 CALL DZERO(WORK(KINT2S),NT1AMX*NRHFT*NRHFT) 200 CALL DZERO(WORK(KINT3S),NT1AMX*NVIRT*NVIRT) 201 CALL DZERO(WORK(KINT4S),NT1AMX*NRHFT*NRHFT) 202 CALL DZERO(WORK(KINT5S),NT1AMX*NVIRT*NVIRT) 203 CALL DZERO(WORK(KINT6S),NT1AMX*NRHFT*NRHFT) 204 CALL DZERO(WORK(KT3AM),NT1AMX*NT1AMX*NT1AMX) 205 CALL DZERO(WORK(KOME1),NT1AMX) 206 CALL DZERO(WORK(KOME2P),NT1AMX*NT1AMX) 207 CALL DZERO(WORK(KOME2M),NT1AMX*NT1AMX) 208 CALL DZERO(WORK(KIAJB),NT1AMX*NT1AMX) 209 CALL DZERO(WORK(KYIAJB),NT1AMX*NT1AMX) 210 CALL DZERO(WORK(KINTTOT),NORBT*NORBT*NORBT*NORBT) 211 CALL DZERO(WORK(KINTTOT2),NORBT*NORBT*NORBT*NORBT) 212! 213 IF (NBAST*NBAST*NBAST*NBAST .GT. 214 * (NRHFT+NVIRT)*NORBT*NORBT*NORBT) THEN 215! 216 CALL DZERO(WORK(KINTTMP1),NBAST*NBAST*NBAST*NBAST) 217 CALL DZERO(WORK(KINTTMP2),NBAST*NBAST*NBAST*NBAST) 218 ELSE 219 CALL DZERO(WORK(KTRTMP1),NRHFT*NORBT*NORBT*NORBT) 220 CALL DZERO(WORK(KTRTMP2),NVIRT*NORBT*NORBT*NORBT) 221 CALL DZERO(WORK(KTRTMP3),NRHFT*NORBT*NORBT*NORBT) 222 CALL DZERO(WORK(KTRTMP4),NVIRT*NORBT*NORBT*NORBT) 223 ENDIF 224! 225!---------------------------------------------- 226! Read caconical orbital energies & 227! MO-coefficients from interface file. 228!---------------------------------------------- 229! 230 IF (LUSIFC .LE. 0) THEN 231 LUSIFC = -1 232 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED', 233 * IDUMMY,.FALSE.) 234 ENDIF 235! 236 REWIND (LUSIFC) 237! 238 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 239! 240 READ (LUSIFC) 241 READ (LUSIFC) (WORK(KFOCKD+I-1), I=1,NORBTS) 242 READ (LUSIFC) (WORK(KCMO+I-1), I=1,NLAMDS) 243! 244 CALL GPCLOSE(LUSIFC,'KEEP') 245! 246!=============================================== 247! Frozen orbital stuff. 248!=============================================== 249! 250 CALL CMO_REORDER(WORK(KCMO),WORK(KEND1),LWRK1) 251! 252 IF (FROIMP .OR. FROEXP) 253 * CALL CCSD_DELFRO(WORK(KFOCKD),WORK(KEND1),LWRK1) 254! 255!==================================================== 256! Start the loop over distributions of integrals. 257!==================================================== 258! 259 IF (DIRECT) THEN 260 NTOSYM = 1 261! DTIME = SECOND() 262!!!! CALL HERDI1(WORK(KEND1),LWRK1,IPRINT) 263 CALL QUIT('Direct not implemented in CC3_TRIPLET') 264! DTIME = SECOND() - DTIME 265! TIMHER1 = TIMHER1 + DTIME 266 ELSE 267 NTOSYM = NSYM 268 ENDIF 269! 270 DO 100 ISYMD1 = 1,NTOSYM 271! 272 IF (DIRECT) THEN 273 NTOT = MAXSHL 274 ELSE 275 NTOT = NBAS(ISYMD1) 276 ENDIF 277! 278 DO 110 ILLL = 1,NTOT 279! 280!----------------------------------------------------------------- 281! If direct calculate the integrals and transposed t2am. 282!----------------------------------------------------------------- 283! 284 IF (DIRECT) THEN 285! DTIME = SECOND() 286!!!! CALL HERDI2(WORK(KEND1),LWRK1,INDEXA,ILLL,NUMDIS,IPRINT) 287 CALL QUIT('Direct not implemented') 288! DTIME = SECOND() - DTIME 289! TIMHER2 = TIMHER2 + DTIME 290! 291 KRECNR = KEND1 292 KEND1 = KRECNR + (NBUFX(0) - 1)/IRAT + 1 293 LWRK1 = LWORK - KEND1 294 IF (LWRK1 .LT. 0) THEN 295 CALL QUIT('Insufficient space in cc3_trip.F (KRECNR)') 296 ENDIF 297 ELSE 298 NUMDIS = 1 299 ENDIF 300! 301!----------------------------------------------------- 302! Loop over number of distributions in disk. 303!----------------------------------------------------- 304! 305 DO 120 IDEL2 = 1,NUMDIS 306! 307 IF (DIRECT) THEN 308 IDEL = INDEXA(IDEL2) 309 IF (NOAUXB) THEN 310 IDUM = 1 311 CALL IJKAUX(IDEL,IDUM,IDUM,IDUM) 312 END IF 313 ISYMD = ISAO(IDEL) 314 ELSE 315 IDEL = IBAS(ISYMD1) + ILLL 316 ISYMD = ISYMD1 317 ENDIF 318! 319 ISYDIS = MULD2H(ISYMD,ISYMOP) 320! 321!------------------------------------------ 322! Work space allocation no. 2. 323!------------------------------------------ 324! 325 KXINT = KEND1 326 KEND2 = KXINT + NDISAO(ISYDIS) 327 LWRK2 = LWORK - KEND2 328! 329 IF (LWRK2 .LT. 0) THEN 330 WRITE(LUPRI,*) 'Need : ',KEND2, 331 * 'Available : ',LWORK 332 CALL QUIT('Insufficient space in CC3_TRIPLET') 333 ENDIF 334! 335!----------------------------------------- 336! Read in batch of integrals. 337!----------------------------------------- 338! 339 DTIME = SECOND() 340 CALL CCRDAO(WORK(KXINT),IDEL,IDEL2,WORK(KEND2),LWRK2, 341 * WORK(KRECNR),DIRECT) 342 DTIME = SECOND() - DTIME 343 TIMRDAO = TIMRDAO + DTIME 344! 345!---------------------------------------------------- 346! Calculate the integrals that are 347! needed in CC3. INTTOT contains T_{1} 348! while INTTOT2 doesnt. 349!---------------------------------------------------- 350! 351 CALL CC3EXCI_INT1(WORK(KCMO),WORK(KXINT), 352 * IDEL,WORK(KINTTOT2), 353 * WORK(KINTTMP1),WORK(KINTTMP2)) 354! 355 CALL CC3EXCI_INT2(WORK(KINTTOT),XLAMDP,XLAMDH, 356 * WORK(KXINT),IDEL, 357 * WORK(KINTTMP1),WORK(KINTTMP2)) 358! 359! 360 120 CONTINUE 361 110 CONTINUE 362 100 CONTINUE 363! 364 IF (LOCDBG) THEN 365! 366 DO I = 1, NORBT*NORBT*NORBT*NORBT 367 IF (ABS(WORK(KINTTOT - 1 + I)) .GT. 1.0D-8) THEN 368 WRITE(LUPRI,*)'Hat{g}(I) ',WORK(KINTTOT - 1 + I) 369 ENDIF 370 ENDDO 371! 372 WRITE(LUPRI,*)' ' 373! 374 DO I = 1, NORBT*NORBT*NORBT*NORBT 375 IF (ABS(WORK(KINTTOT2 - 1 + I)) .GT. 1.0D-8) THEN 376 WRITE(LUPRI,*)'g(I) ',WORK(KINTTOT2 - 1 + I) 377 ENDIF 378 ENDDO 379 ENDIF 380! 381! 382 IF (LOCDBG) THEN 383 WRITE(LUPRI,*)'After the integral section' 384 ENDIF 385! 386!------------------------------------------------------- 387! Transform the integrals from g(p,q,r,s) to 388! whatever is needed 389!------------------------------------------------------- 390! 391 CALL CC33_TRANSTEMP(WORK(KINTTOT),C1AM,WORK(KINT1S), 392 * WORK(KINT2S),WORK(KINT3S), 393 * WORK(KINT4S),WORK(KINT5S), 394 * WORK(KINT6S),WORK(KINTTOT2), 395 * WORK(KIAJB),WORK(KYIAJB), 396 * WORK(KINT1T),WORK(KINT2T), 397 * WORK(KINT3T),WORK(KINT4T), 398 * WORK(KTRTMP1),WORK(KTRTMP2), 399 * WORK(KTRTMP3),WORK(KTRTMP4)) 400! 401! 402!--------------------------------------- 403! Calculate the triple amplitudes. 404!--------------------------------------- 405! 406 IF (LOCDBG) THEN 407 WRITE(LUPRI,*)'Before _T3AM' 408 CALL FLSHFO(LUPRI) 409 ENDIF 410! 411 CALL CCSDT_T3AM(WORK(KT3AM),WORK(KINT1S),WORK(KINT2S),T2AM, 412 * WORK(KSCR1),WORK(KFOCKD)) 413! 414!--------------------------------------------------------- 415! Calculate the corrections from the triple ampl.. 416!--------------------------------------------------------- 417! 418 IF (LOCDBG) THEN 419 WRITE(LUPRI,*)'Before OMEGA2T3AM' 420 CALL FLSHFO(LUPRI) 421 ENDIF 422! 423 CALL CC33_OMEGA2T3(WORK(KOME2P),WORK(KOME2M),WORK(KINT3T), 424 * WORK(KINT4T),WORK(KT3AM)) 425! 426! 427 IF (LOCDBG) THEN 428 WRITE(LUPRI,*)'After corr. with t3; norm of 2 omegas :' 429 TEMP = NT1AMX*NT1AMX 430 RHO2N = DDOT(TEMP,WORK(KOME2P),1,WORK(KOME2P),1) 431 WRITE(LUPRI,*)'The (+) :',RHO2N 432 RHO2N = DDOT(TEMP,WORK(KOME2M),1,WORK(KOME2M),1) 433 WRITE(LUPRI,*)'The (-) :',RHO2N 434 ENDIF 435! 436!-------------------------------------------------------------- 437! Reset the array and then 438! calculate the triplet triple trialvector amplitudes. 439!-------------------------------------------------------------- 440! 441 IF (LOCDBG) THEN 442 WRITE(LUPRI,*)'Before _C3AM' 443 CALL FLSHFO(6) 444 ENDIF 445! 446 CALL DZERO(WORK(KT3AM),NT1AMX*NT1AMX*NT1AMX) 447! 448 CALL CC33_C3AM(WORK(KT3AM),WORK(KINT1S),WORK(KINT2S), 449 * WORK(KINT3S),WORK(KINT4S),WORK(KINT5S), 450 * WORK(KINT6S),T2AM,C2AMP,C2AMM,WORK(KSCR1), 451 * WORK(KFOCKD)) 452! 453 IF (LOCDBG) THEN 454 TEMP = NT1AMX*NT1AMX*NT1AMX 455 RHO2N = DDOT(TEMP,WORK(KT3AM),1,WORK(KT3AM),1) 456 WRITE(LUPRI,*)'The norm of the C3 vector is:',RHO2N 457 ENDIF 458! 459!------------------------------------------------------------- 460! Calculate the corrections from the triple trialvector 461! to the 1 excited result vector. 462!------------------------------------------------------------- 463! 464 IF (LOCDBG) THEN 465 RHO2N = DDOT(NT1AMX,WORK(KOME1),1,WORK(KOME1),1) 466 WRITE(LUPRI,*)'The norm of omega1 at the beginning.:',RHO2N 467 WRITE(LUPRI,*)'Before OMEGA1C3AM' 468 CALL FLSHFO(LUPRI) 469 ENDIF 470! 471 CALL CC33_OMEGA1C3(WORK(KOME1),WORK(KIAJB), 472 * WORK(KYIAJB),WORK(KT3AM)) 473! 474 IF (LOCDBG) THEN 475 RHO2N = DDOT(NT1AMX,WORK(KOME1),1,WORK(KOME1),1) 476 WRITE(LUPRI,*)'The norm of omega1 after the corr. :',RHO2N 477 ENDIF 478! 479!------------------------------------------------------------ 480! Sum the T3-results into the CCSD results and go 481! on to calculate the C3 results. 482!------------------------------------------------------------ 483! 484! 485 CALL DZERO(WORK(KOME2),NT1AMX*NT1AMX) 486! 487 CALL CC_RVEC3(LUFR2,FRHO2,NT2AM(ISYMTR)+NT2AMA(ISYMTR), 488 * NT2AMA(ISYMTR),IV+ITR-1,NT2AM(ISYMTR), 489 * WORK(KOME2)) 490! 491 IF (LOCDBG) THEN 492 TEMP = NT2AM(ISYMTR) 493 RHO2N = DDOT(TEMP,WORK(KOME2),1,WORK(KOME2),1) 494 WRITE(LUPRI,*)'Norm of (-)-result from CCSD',RHO2N 495 RHO2N = DDOT(TEMP,OMEGA2P(1),1,OMEGA2P(1),1) 496 WRITE(LUPRI,*)'Norm of (+)-result from CCSD',RHO2N 497 WRITE(LUPRI,*)'The Omega2 after the readin of Rho2-.' 498 TEMP = NT1AMX*NT1AMX 499 RHO2N = DDOT(TEMP,WORK(KOME2P),1,WORK(KOME2P),1) 500 WRITE(LUPRI,*)'The norm of the (+)-vector:',RHO2N 501 RHO2N = DDOT(TEMP,WORK(KOME2M),1,WORK(KOME2M),1) 502 WRITE(LUPRI,*)'The norm of the (-)-vector:',RHO2N 503 ENDIF 504! 505!------------------------------------ 506! Sum Omega2 (+ & -) for T3. 507!------------------------------------ 508! 509 DO 310 I = 1,NT1AMX 510 DO 320 J = 1,I 511 IJ = NT1AMX*(I-1) + J 512 NIJ = INDEX(I,J) 513! 514 OMEGA2P(NIJ) = OMEGA2P(NIJ) + WORK(KOME2P+IJ-1) 515! 516 WORK(KOME2-1+NIJ) = WORK(KOME2-1+NIJ) + WORK(KOME2M+IJ-1) 517! 518 320 CONTINUE 519 310 CONTINUE 520! 521 IF (LOCDBG) THEN 522 TEMP = NT2AM(ISYMTR) 523 RHO2N = DDOT(TEMP,OMEGA2P,1,OMEGA2P,1) 524 WRITE(LUPRI,*) 525 * 'The norm of the (+) vector after CCSD & T3',RHO2N 526 TEMP = NT2AMA(ISYMTR) 527 RHO2N = DDOT(TEMP,WORK(KOME2),1,WORK(KOME2),1) 528 WRITE(LUPRI,*) 529 * 'The norm of the (-) vector after CCSD & T3',RHO2N 530 ENDIF 531! 532!------------------------------------------------------------ 533! Proceed with the calculation omega2. Calculate the 534! C3 dependent terms, but first zero the results vectors. 535!------------------------------------------------------------ 536! 537 IF (LOCDBG) THEN 538 WRITE(LUPRI,*)'Before OMEGA2C3AM' 539 CALL FLSHFO(6) 540 ENDIF 541! 542 CALL DZERO(WORK(KOME2P),NT1AMX*NT1AMX) 543 CALL DZERO(WORK(KOME2M),NT1AMX*NT1AMX) 544! 545 CALL CC33_OMEGA2C3(WORK(KOME2P),WORK(KOME2M),WORK(KINT1T), 546 * WORK(KINT2T),WORK(KT3AM),FOCK) 547! 548 IF (LOCDBG) THEN 549 WRITE(LUPRI,*)'The Omega2 after the C3 corr.' 550 TEMP = NT1AMX*NT1AMX 551 RHO2N = DDOT(TEMP,WORK(KOME2P),1,WORK(KOME2P),1) 552 WRITE(LUPRI,*)'The norm of the (+)-vector:',RHO2N 553 RHO2N = DDOT(TEMP,WORK(KOME2M),1,WORK(KOME2M),1) 554 WRITE(LUPRI,*)'The norm of the (-)-vector:',RHO2N 555 ENDIF 556! 557!-------------------------------- 558! Sum Omega1 559!-------------------------------- 560! 561 DO 400 I = 1,NT1AMX 562 OMEGA1(I) = OMEGA1(I) + WORK(KOME1+I-1) 563 400 CONTINUE 564! 565 IF (LOCDBG) THEN 566 RHO2N = DDOT(NT1AMX,WORK(KOME1),1,WORK(KOME1),1) 567 WRITE(LUPRI,*) 568 * 'The norm of omega1 after the summation. :',RHO2N 569 ENDIF 570! 571!------------------------------------------------------ 572! Sum Omega2 (+ & -) for the C3 dependent terms 573!------------------------------------------------------ 574! 575 IF (LOCDBG) THEN 576 WRITE(LUPRI,*)'Before the final summation' 577 WRITE(LUPRI,*)' ' 578 CALL FLSHFO(LUPRI) 579 ENDIF 580! 581 DO 410 I = 1,NT1AMX 582 DO 420 J = 1,I 583 IJ = NT1AMX*(I-1) + J 584 NIJ = INDEX(I,J) 585! 586 OMEGA2P(NIJ) = OMEGA2P(NIJ) + WORK(KOME2P+IJ-1) 587! 588 WORK(KOME2-1+NIJ) = WORK(KOME2-1+NIJ) + WORK(KOME2M+IJ-1) 589! 590 420 CONTINUE 591 410 CONTINUE 592! 593 IF (LOCDBG) THEN 594 TEMP = NT2AM(ISYMTR) 595 RHO2N = DDOT(TEMP,OMEGA2P,1,OMEGA2P,1) 596 WRITE(LUPRI,*) 597 * 'The norm of the (+) vector at the end',RHO2N 598 TEMP = NT2AMA(ISYMTR) 599 RHO2N = DDOT(TEMP,WORK(KOME2),1,WORK(KOME2),1) 600 WRITE(LUPRI,*) 601 * 'The norm of the (-) vector at the end',RHO2N 602 WRITE(LUPRI,*)'End of triple part of the CC3 calc.' 603 ENDIF 604! 605!------------------------------- 606! Write Omega2- to file. 607!------------------------------- 608! 609 CALL CC_WVEC3(LUFR2,FRHO2,NT2AM(ISYMTR)+NT2AMA(ISYMTR), 610 * NT2AMA(ISYMTR),IV + ITR -1,NT2AM(ISYMTR), 611 * WORK(KOME2)) 612! 613!-------------------------- 614! End the calc. 615!-------------------------- 616! 617 RETURN 618 END 619C /* Deck cc3_c3am */ 620 SUBROUTINE CC33_C3AM(ECURR,C3AM, 621 * XINT1S,XINT2S,XINT3S,XINT4S,XINT5S, 622 * XINT6S,T2AM,C2AMP,C2AMM,SCR1,FOCKD) 623! 624! Written by Kasper Hald Jan 2000. 625! 626 IMPLICIT NONE 627! 628#include "priunit.h" 629#include "ccorb.h" 630#include "ccsdinp.h" 631#include "ccsdsym.h" 632! 633 INTEGER NAI, NCK, NCJ, NBK, NBJ, NDI, NDJ, NDK 634 INTEGER NBL, NAL, NCL, NAK, NAJ, NCI, NBI 635! 636#if defined (SYS_CRAY) 637 REAL XINT1S(NT1AMX,NVIRT,NVIRT), XINT2S(NT1AMX,NRHFT,NRHFT) 638 REAL XINT3S(NT1AMX,NVIRT,NVIRT), XINT4S(NT1AMX,NRHFT,NRHFT) 639 REAL XINT5S(NT1AMX,NVIRT,NVIRT), XINT6S(NT1AMX,NRHFT,NRHFT) 640 REAL C3AM(NT1AMX,NT1AMX,NT1AMX),SCR1(NT1AMX),FOCKD(*) 641 REAL T2AM(NT1AMX,NT1AMX), HALF, QUART, ONE, ECURR 642 REAL AIBJCK, C2AMP(NT1AMX,NT1AMX), C2AMM(NT1AMX,NT1AMX) 643#else 644 DOUBLE PRECISION XINT1S(NT1AMX,NVIRT,NVIRT) 645 DOUBLE PRECISION XINT2S(NT1AMX,NRHFT,NRHFT) 646 DOUBLE PRECISION XINT3S(NT1AMX,NVIRT,NVIRT) 647 DOUBLE PRECISION XINT4S(NT1AMX,NRHFT,NRHFT) 648 DOUBLE PRECISION XINT5S(NT1AMX,NVIRT,NVIRT) 649 DOUBLE PRECISION XINT6S(NT1AMX,NRHFT,NRHFT) 650 DOUBLE PRECISION C3AM(NT1AMX,NT1AMX,NT1AMX),SCR1(NT1AMX),FOCKD(*) 651 DOUBLE PRECISION T2AM(NT1AMX,NT1AMX), HALF, QUART, AIBJCK 652 DOUBLE PRECISION C2AMP(NT1AMX,NT1AMX), C2AMM(NT1AMX,NT1AMX) 653 DOUBLE PRECISION ONE, ECURR 654#endif 655! 656 PARAMETER (QUART = 0.25D00, HALF=0.5D00, ONE=1.0D00) 657! 658 DO I = 1,NRHFT 659 DO A = 1,NVIRT 660 NAI = NVIRT*(I-1) + A 661 SCR1(NAI) = FOCKD(NRHFT+A) - FOCKD(I) 662 ENDDO 663 ENDDO 664! 665! 666 DO 100 K = 1, NRHFT 667 DO 110 C = 1, NVIRT 668 NCK = NVIRT*(K-1) + C 669! 670 DO 120 J = 1,NRHFT 671! 672 NCJ = NVIRT*(J-1) + C 673 DO 130 B = 1,NVIRT 674! 675 NBJ = NVIRT*(J-1) + B 676 NBK = NVIRT*(K-1) + B 677! 678 DO 140 I = 1,NRHFT 679 NBI = NVIRT*(I-1) + B 680 NCI = NVIRT*(I-1) + C 681 DO 150 A = 1,NVIRT 682 NAI = NVIRT*(I-1) + A 683 NAJ = NVIRT*(J-1) + A 684 NAK = NVIRT*(K-1) + A 685! 686 AIBJCK = 0.0D0 687 DO 160 D = 1,NVIRT 688! 689 NDI = NVIRT*(I-1) + D 690 NDJ = NVIRT*(J-1) + D 691 NDK = NVIRT*(K-1) + D 692! 693 AIBJCK = AIBJCK + HALF*XINT1S(NAI,C,D)* 694 * C2AMP(NDK,NBJ) 695! 696 AIBJCK = AIBJCK + ONE*XINT1S(NBJ,A,D)* 697 * C2AMM(NCK,NDI) 698! 699 AIBJCK = AIBJCK + ONE*XINT1S(NCJ,B,D)* 700 * C2AMM(NAI,NDK) 701! 702 AIBJCK = AIBJCK + HALF*XINT3S(NCJ,B,D)*T2AM(NAI,NDK) 703! 704 AIBJCK = AIBJCK + HALF*XINT5S(NAI,C,D)*T2AM(NBJ,NDK) 705! 706 AIBJCK = AIBJCK + HALF*XINT5S(NCJ,A,D)*T2AM(NBK,NDI) 707! 708 160 CONTINUE 709! 710 DO 170 L = 1,NRHFT 711! 712 NAL = NVIRT*(L-1) + A 713 NBL = NVIRT*(L-1) + B 714 NCL = NVIRT*(L-1) + C 715! 716 AIBJCK = AIBJCK - HALF*XINT2S(NAI,L,K)* 717 * C2AMP(NCL,NBJ) 718! 719 AIBJCK = AIBJCK + ONE*XINT2S(NBJ,L,I)* 720 * C2AMM(NAL,NCK) 721! 722 AIBJCK = AIBJCK + ONE*XINT2S(NCJ,L,K)* 723 * C2AMM(NBL,NAI) 724! 725 AIBJCK = AIBJCK - HALF*XINT4S(NBK,L,J)*T2AM(NAI,NCL) 726! 727 AIBJCK = AIBJCK - HALF*XINT6S(NAI,L,K)*T2AM(NBJ,NCL) 728! 729 AIBJCK = AIBJCK - HALF*XINT6S(NCJ,L,I)*T2AM(NBK,NAL) 730! 731 170 CONTINUE 732! 733 AIBJCK = AIBJCK/(ECURR-SCR1(NAI)-SCR1(NBJ)-SCR1(NCK)) 734! 735! 736 C3AM(NAI,NBJ,NCK) = C3AM(NAI,NBJ,NCK) + AIBJCK 737 C3AM(NAI,NBK,NCJ) = C3AM(NAI,NBK,NCJ) - AIBJCK 738 C3AM(NAI,NCK,NBJ) = C3AM(NAI,NCK,NBJ) + AIBJCK 739 C3AM(NAI,NCJ,NBK) = C3AM(NAI,NCJ,NBK) - AIBJCK 740! 741 150 CONTINUE 742 140 CONTINUE 743 130 CONTINUE 744 120 CONTINUE 745 110 CONTINUE 746 100 CONTINUE 747! 748! 749 RETURN 750 END 751C /* Deck cc3_omega1c3 */ 752 SUBROUTINE CC33_OMEGA1C3(OMEG1,XIAJB,YIAJB,C3AM) 753! 754! Written by Kasper Hald Jan. 2000 755! 756 IMPLICIT NONE 757! 758#include "ccorb.h" 759#include "ccsdinp.h" 760#include "ccsdsym.h" 761! 762 INTEGER NAI, NAJ, NBJ, NBK, NCK, NCI 763 INTEGER VIR1, VIR2, OCC1, OCC2 764! 765#if defined (SYS_CRAY) 766 REAL TWO 767 REAL OMEG1(NT1AMX),XIAJB(NT1AMX,NT1AMX) 768 REAL YIAJB(NT1AMX,NT1AMX) 769 REAL C3AM(NT1AMX,NT1AMX,NT1AMX) 770#else 771 DOUBLE PRECISION TWO 772 DOUBLE PRECISION OMEG1(NT1AMX),XIAJB(NT1AMX,NT1AMX) 773 DOUBLE PRECISION YIAJB(NT1AMX,NT1AMX) 774 DOUBLE PRECISION C3AM(NT1AMX,NT1AMX,NT1AMX) 775#endif 776! 777 PARAMETER (TWO = 2.0D0) 778! 779! 780 DO 100 I = 1,NRHFT 781 DO 110 A = 1,NVIRT 782 NAI = NVIRT*(I-1) + A 783! 784 DO 120 J = 1,NRHFT 785 NAJ = NVIRT*(J-1) + A 786 DO 130 B = 1,NVIRT 787! 788 NBJ = NVIRT*(J-1) + B 789! 790 DO 140 K = 1,NRHFT 791 NBK = NVIRT*(K-1) + B 792 DO 150 C = 1,NVIRT 793 NCK = NVIRT*(K-1) + C 794 NCI = NVIRT*(I-1) + C 795! 796! 797 OMEG1(NAI) = OMEG1(NAI) + TWO*XIAJB(NCK,NBJ)* 798 * C3AM(NBJ,NAI,NCK) 799 * + TWO*YIAJB(NCK,NBJ)* 800 * (C3AM(NAJ,NBK,NCI)+ 801 * C3AM(NCI,NBK,NAJ)) 802! 803! 804 150 CONTINUE 805 140 CONTINUE 806 130 CONTINUE 807 120 CONTINUE 808 110 CONTINUE 809 100 CONTINUE 810! 811 RETURN 812 END 813C /* Deck cc33_omega2t3 */ 814 SUBROUTINE CC33_OMEGA2T3(OMEGA2P,OMEGA2M,XINT3T,XINT4T,T3AM) 815! 816! Written by Kasper Hald Jan. 2000. 817! 818! Calculate the T3,R1 contribution to omega2+: 819! .5*P(ab,ij) (tilde)P(ij) 820! [sum(dlm) t(ai,dj,bl)*g(lm-bar,md) + 821! sum(dlm)(t(am,bl,di)+t(ai,dl,bm)-2t(ai,dm,bl))*g(lj-bar,md) + 822! sum(del)(2t(ai,dl,ej)-t(al,di,ej)-t(ai,el,dj))*g(b-bar e,ld) 823! ] 824! N.B. Do not calculate with (Tilde)P(ij). 825! 826! Calculate the T3,R1 contribution to omega2-: 827! .5*(tilde)P(ab,ij) 828! [sum(dlm) t(ai,dj,bl)*g(lm-bar,md) + 829! sum(dlm)(t(am,bl,di)+t(ai,dl,bm)-2t(ai,dm,bl))*g(lj-bar,md) + 830! sum(del)(2t(ai,dl,ej)-t(al,di,ej)-t(ai,el,dj))*g(b-bar e,ld) 831! ] 832! 833 IMPLICIT NONE 834! 835#include "priunit.h" 836#include "ccorb.h" 837#include "ccsdinp.h" 838#include "ccsdsym.h" 839! 840 INTEGER NBK, NAK, NAI, NBJ, NCK, NCJ, NCI, NDJ 841 INTEGER NDK, NDI, NBL, NCL, NAL 842 INTEGER NAJ, NBI 843! 844#if defined (SYS_CRAY) 845 REAL XINT3T(NT1AMX,NVIRT,NVIRT) 846 REAL XINT4T(NT1AMX,NRHFT,NRHFT) 847 REAL OMEGA2P(NT1AMX,NT1AMX) 848 REAL OMEGA2M(NT1AMX,NT1AMX) 849 REAL T3AM(NT1AMX,NT1AMX,NT1AMX) 850 REAL HALF,TWO, XAIBJ, ZERO 851#else 852 DOUBLE PRECISION XINT3T(NT1AMX,NVIRT,NVIRT) 853 DOUBLE PRECISION XINT4T(NT1AMX,NRHFT,NRHFT) 854 DOUBLE PRECISION OMEGA2P(NT1AMX,NT1AMX) 855 DOUBLE PRECISION OMEGA2M(NT1AMX,NT1AMX) 856 DOUBLE PRECISION T3AM(NT1AMX,NT1AMX,NT1AMX) 857 DOUBLE PRECISION HALF,TWO, XAIBJ, ZERO, XHALF 858#endif 859! 860 PARAMETER (ZERO = 0.0D00, HALF = 0.5D00, TWO = 2.0D0) 861 PARAMETER (XHALF=-0.5D0) 862! 863! 864 DO A = 1, NVIRT 865 DO I = 1, NRHFT 866 NAI = NVIRT*(I-1) + A 867! 868 DO J = 1, NRHFT 869 NAJ = NVIRT*(J-1) + A 870 DO K = 1, NRHFT 871 NAK = NVIRT*(K-1) + A 872 T3AM(NAI,NAJ,NAK) = 0.0D0 873 T3AM(NAI,NAK,NAJ) = 0.0D0 874 T3AM(NAJ,NAI,NAK) = 0.0D0 875 T3AM(NAJ,NAK,NAI) = 0.0D0 876 T3AM(NAK,NAJ,NAI) = 0.0D0 877 T3AM(NAK,NAI,NAJ) = 0.0D0 878 ENDDO 879 ENDDO 880! 881 DO B = 1, NVIRT 882 NBI = NVIRT*(I-1)+B 883 DO C = 1, NVIRT 884 NCI = NVIRT*(I-1)+C 885! 886 T3AM(NAI,NBI,NCI) = 0.0D0 887 T3AM(NAI,NCI,NBI) = 0.0D0 888 T3AM(NBI,NAI,NCI) = 0.0D0 889 T3AM(NBI,NCI,NAI) = 0.0D0 890 T3AM(NCI,NAI,NBI) = 0.0D0 891 T3AM(NCI,NBI,NAI) = 0.0D0 892 ENDDO 893 ENDDO 894 ENDDO 895 ENDDO 896! 897! 898 DO 100 I = 1,NRHFT 899 DO 110 A = 1,NVIRT 900 NAI = NVIRT*(I-1) + A 901! 902 DO 120 J = 1,NRHFT 903 DO 130 B = 1,NVIRT 904 NBJ = NVIRT*(J-1) + B 905! 906 IF (NAI .NE. NBJ) THEN 907! 908 DO 140 K = 1,NRHFT 909 NBK = NVIRT*(K-1) + B 910 NAK = NVIRT*(K-1) + A 911 DO 150 C = 1,NVIRT 912! 913 NCK = NVIRT*(K-1) + C 914 NCJ = NVIRT*(J-1) + C 915 NCI = NVIRT*(I-1) + C 916! 917 DO 160 D = 1,NVIRT 918! 919 NDJ = NVIRT*(J-1) + D 920 NDK = NVIRT*(K-1) + D 921 NDI = NVIRT*(I-1) + D 922! 923 XAIBJ = (TWO*T3AM(NAI,NDK,NCJ) - 924 * T3AM(NAK,NDI,NCJ) - 925 * T3AM(NAI,NCK,NDJ))*XINT3T(NDK,B,C) 926! 927 OMEGA2M(NAI,NBJ) = OMEGA2M(NAI,NBJ) + XAIBJ 928! 929 IF ((A .NE. B) .AND. (I .NE. J)) THEN 930 OMEGA2P(NAI,NBJ) = OMEGA2P(NAI,NBJ) + XAIBJ 931 ENDIF 932! 933 160 CONTINUE 934! 935 DO 170 L = 1,NRHFT 936! 937 NBL = NVIRT*(L-1) + B 938 NCL = NVIRT*(L-1) + C 939 NAL = NVIRT*(L-1) + A 940! 941 XAIBJ = T3AM(NAI,NCJ,NBL)*XINT4T(NCK,L,K) 942 * +(T3AM(NAK,NBL,NCI) + 943 * T3AM(NAI,NCL,NBK) - 944 * TWO*T3AM(NAI,NCK,NBL))*XINT4T(NCK,L,J) 945! 946 OMEGA2M(NAI,NBJ) = OMEGA2M(NAI,NBJ) + XAIBJ 947! 948 IF ((A .NE. B) .AND. (I .NE. J)) THEN 949 OMEGA2P(NAI,NBJ) = OMEGA2P(NAI,NBJ) + XAIBJ 950 ENDIF 951! 952 170 CONTINUE 953! 954 150 CONTINUE 955 140 CONTINUE 956! 957 ENDIF 958! 959 130 CONTINUE 960 120 CONTINUE 961! 962 110 CONTINUE 963 100 CONTINUE 964! 965 DO 200 NAI = 1,NT1AMX 966 DO 210 NBJ = 1,NAI 967! 968! The scaling of (+) with 1/2 is done in the end 969! of the CCSD routine. 970! 971 IF (NAI .NE. NBJ) THEN 972 XAIBJ = OMEGA2P(NAI,NBJ) + OMEGA2P(NBJ,NAI) 973 OMEGA2P(NAI,NBJ) = XAIBJ 974 OMEGA2P(NBJ,NAI) = XAIBJ 975 XAIBJ = HALF*(OMEGA2M(NAI,NBJ) - OMEGA2M(NBJ,NAI)) 976 OMEGA2M(NBJ,NAI) = XAIBJ 977 OMEGA2M(NAI,NBJ) = -XAIBJ 978 ELSE 979 OMEGA2M(NAI,NBJ) = ZERO 980 ENDIF 981! 982 210 CONTINUE 983 200 CONTINUE 984! 985 RETURN 986 END 987C /* Deck cc33_omega2c3 */ 988 SUBROUTINE CC33_OMEGA2C3(OMEGA2P,OMEGA2M,XINT1T,XINT2T,C3AM,FOCK) 989! 990! Written by Kasper Hald Jan. 2000. 991! 992! Calc. the cont. : 993! 994! (+): .5*P(ab,ij)[sum(dl)( c3(dl,bj,ai)+2*c3(ai,bj,dl)+ 995! c3(di,bl,aj)+c3(bl,di,aj) )*F(ld) 996! +sum(del)( 2*c3(dl,ei,bj)-c3(el,di,bj)+ 997! 2*c3(bj,ei,dl)+2*c3(ei,bj,dl)+ 998! 2*c3(dj,bi,el)-c3(bl,dj,ei) )*g(ld,ae) 999! +sum(dlm)( 2*c3(dl,am,bi)-c3(dm,al,bi)+ 1000! 2*c3(bi,am,dl)+2*c3(am,bi,dl) + 1001! 2*c3(al,di,bm)-c3(di,bl,am) )*g(ld,mj) 1002! ] 1003! 1004! (-): (Tilde)P(ab,ij)[sum(ld) c3((ai,bj,dl)*F(ld)+ 1005! sum(del) (c3(bj,di,el)+c3(ei,bj,dl))*g(ae,ld) 1006! sum(dlm) (c3(ai,dm,bl)+c3(bm,ai,dl))*g(ld,mj) 1007! ] 1008! 1009 IMPLICIT NONE 1010! 1011#include "ccorb.h" 1012#include "ccsdinp.h" 1013#include "ccsdsym.h" 1014! 1015 INTEGER NAI,NAJ,NBI,NBJ,NBK,NCK,NCI,NDJ,NDK,NDI 1016 INTEGER NBL, NAL, NCL, NAK 1017! 1018#if defined (SYS_CRAY) 1019 REAL XINT1T(NT1AMX,NVIRT,NVIRT) 1020 REAL XINT2T(NT1AMX,NRHFT,NRHFT) 1021 REAL OMEGA2P(NT1AMX,NT1AMX), HALF, ZERO 1022 REAL OMEGA2M(NT1AMX,NT1AMX), TWO, XAIBJ 1023 REAL C3AM(NT1AMX,NT1AMX,NT1AMX),FOCK(NORBT,NORBT) 1024#else 1025 DOUBLE PRECISION XINT1T(NT1AMX,NVIRT,NVIRT) 1026 DOUBLE PRECISION XINT2T(NT1AMX,NRHFT,NRHFT) 1027 DOUBLE PRECISION OMEGA2P(NT1AMX,NT1AMX), HALF, ZERO 1028 DOUBLE PRECISION OMEGA2M(NT1AMX,NT1AMX), TWO, XAIBJ 1029 DOUBLE PRECISION C3AM(NT1AMX,NT1AMX,NT1AMX),FOCK(NORBT,NORBT) 1030#endif 1031! 1032 PARAMETER (ZERO = 0.0D00, HALF = 0.5D00, TWO = 2.0D0) 1033! 1034! 1035 DO 100 I = 1,NRHFT 1036 DO 110 A = 1,NVIRT 1037 NAI = NVIRT*(I-1) + A 1038! 1039 DO 120 J = 1,NRHFT 1040 NAJ = NVIRT*(J-1) + A 1041 DO 130 B = 1,NVIRT 1042 NBI = NVIRT*(I-1) + B 1043 NBJ = NVIRT*(J-1) + B 1044! 1045 IF (NAI .NE. NBJ) THEN 1046! 1047 DO 140 K = 1,NRHFT 1048 NAK = NVIRT*(K-1) + A 1049 NBK = NVIRT*(K-1) + B 1050 DO 150 C = 1,NVIRT 1051! 1052 NCK = NVIRT*(K-1) + C 1053 NCI = NVIRT*(I-1) + C 1054! 1055!=================================== 1056! The plus part of omega2 1057!=================================== 1058! 1059 IF ((A .NE. B) .AND. (I .NE. J)) THEN 1060! 1061 OMEGA2P(NAI,NBJ) = OMEGA2P(NAI,NBJ) + 1062 * (C3AM(NCK,NBJ,NAI)+TWO*C3AM(NAI,NBJ,NCK)+ 1063 * C3AM(NCI,NBK,NAJ)+C3AM(NBK,NCI,NAJ))*FOCK(K,NRHFT+C) 1064! 1065 ENDIF 1066! 1067!=================================== 1068! The minus part of omega2 1069!=================================== 1070! 1071 OMEGA2M(NAI,NBJ) = OMEGA2M(NAI,NBJ) + 1072 * C3AM(NAI,NBJ,NCK)*FOCK(K,NRHFT+C) 1073! 1074!----------------------------------- 1075!----------------------------------- 1076! 1077 DO 160 D = 1,NVIRT 1078! 1079 NDJ = NVIRT*(J-1) + D 1080 NDK = NVIRT*(K-1) + D 1081 NDI = NVIRT*(I-1) + D 1082! 1083!=================================== 1084! The plus part of omega2 1085!=================================== 1086! 1087 IF ((A .NE. B) .AND. (I .NE. J)) THEN 1088! 1089 OMEGA2P(NAI,NBJ) = OMEGA2P(NAI,NBJ) + 1090 * (TWO*C3AM(NDK,NCI,NBJ)-C3AM(NCK,NDI,NBJ)+TWO*C3AM(NBJ,NCI,NDK) 1091 * +TWO*C3AM(NCI,NBJ,NDK)+TWO*C3AM(NDJ,NBI,NCK)-C3AM(NBK,NDJ,NCI)) 1092 * *XINT1T(NDK,A,C) 1093! 1094 ENDIF 1095! 1096!=================================== 1097! The minus part of omega2 1098!=================================== 1099! 1100 OMEGA2M(NAI,NBJ) = OMEGA2M(NAI,NBJ) + 1101 * (C3AM(NBJ,NDI,NCK)+C3AM(NCI,NBJ,NDK))*XINT1T(NDK,A,C) 1102! 1103 160 CONTINUE 1104! 1105!----------------------------------- 1106!----------------------------------- 1107! 1108 DO 170 L = 1,NRHFT 1109! 1110 NBL = NVIRT*(L-1) + B 1111 NCL = NVIRT*(L-1) + C 1112 NAL = NVIRT*(L-1) + A 1113! 1114!=================================== 1115! The plus part of omega2 1116!=================================== 1117! 1118 IF ((A .NE. B) .AND. (I .NE. J)) THEN 1119! 1120 OMEGA2P(NAI,NBJ) = OMEGA2P(NAI,NBJ) + 1121 * (TWO*C3AM(NCL,NAK,NBI)-C3AM(NCK,NAL,NBI)+TWO*C3AM(NBI,NAK,NCL) 1122 * +TWO*C3AM(NAK,NBI,NCL)+TWO*C3AM(NAL,NCI,NBK)-C3AM(NCI,NBL,NAK)) 1123 * *XINT2T(NCL,K,J) 1124! 1125 ENDIF 1126! 1127!=================================== 1128! The minus part of omega2 1129!=================================== 1130! 1131 OMEGA2M(NAI,NBJ) = OMEGA2M(NAI,NBJ) + 1132 * (C3AM(NAI,NCK,NBL)+C3AM(NBK,NAI,NCL))*XINT2T(NCL,K,J) 1133! 1134 170 CONTINUE 1135! 1136 150 CONTINUE 1137 140 CONTINUE 1138! 1139 ENDIF 1140! 1141 130 CONTINUE 1142 120 CONTINUE 1143! 1144 110 CONTINUE 1145 100 CONTINUE 1146! 1147 DO 200 NAI = 1,NT1AMX 1148 DO 210 NBJ = 1,NAI 1149! 1150! (Anti)symmetrisation 1151! 1152! The Half for + is done in the CCSD part. 1153! 1154 IF (NAI .NE. NBJ) THEN 1155 XAIBJ = OMEGA2P(NAI,NBJ)+OMEGA2P(NBJ,NAI) 1156 OMEGA2P(NAI,NBJ) = XAIBJ 1157 OMEGA2P(NBJ,NAI) = XAIBJ 1158 XAIBJ = OMEGA2M(NAI,NBJ) - OMEGA2M(NBJ,NAI) 1159 OMEGA2M(NAI,NBJ) = -XAIBJ 1160 OMEGA2M(NBJ,NAI) = XAIBJ 1161 ELSE 1162 OMEGA2P(NAI,NBJ) = ZERO 1163 OMEGA2M(NAI,NBJ) = ZERO 1164 ENDIF 1165! 1166 210 CONTINUE 1167 200 CONTINUE 1168! 1169 RETURN 1170 END 1171C /* Deck ccsdt_t3am */ 1172 SUBROUTINE CCSDT_T3AM(T3AM,XINT1,XINT2,T2AM,SCR1,FOCKD) 1173! 1174! Written by Henrik Koch. 1175! 1176! Calculates : T3 1177! 1178! N.B. 1179! This routine is also present in the file ccsd_triple.F 1180! which however is not a normal part of the cc-program. 1181! 1182 IMPLICIT NONE 1183! 1184#include "priunit.h" 1185#include "ccorb.h" 1186#include "ccsdinp.h" 1187#include "ccsdsym.h" 1188! 1189 INTEGER NAI, NCK, NBJ, NDJ, NBL 1190! 1191#if defined (SYS_CRAY) 1192 REAL XINT1(NT1AMX,NVIRT,NVIRT), XINT2(NT1AMX,NRHFT,NRHFT) 1193 REAL T3AM(NT1AMX,NT1AMX,NT1AMX),SCR1(NT1AMX),FOCKD(*) 1194 REAL T2AM(NT1AMX,NT1AMX), AIBJCK 1195#else 1196 DOUBLE PRECISION XINT1(NT1AMX,NVIRT,NVIRT) 1197 DOUBLE PRECISION XINT2(NT1AMX,NRHFT,NRHFT) 1198 DOUBLE PRECISION T3AM(NT1AMX,NT1AMX,NT1AMX) 1199 DOUBLE PRECISION SCR1(NT1AMX),FOCKD(*) 1200 DOUBLE PRECISION T2AM(NT1AMX,NT1AMX), AIBJCK 1201#endif 1202! 1203! 1204 DO 50 I = 1,NRHFT 1205 DO 60 A = 1,NVIRT 1206 NAI = NVIRT*(I-1) + A 1207 SCR1(NAI) = FOCKD(NRHFT+A) - FOCKD(I) 1208 60 CONTINUE 1209 50 CONTINUE 1210! 1211 DO 100 NCK = 1,NT1AMX 1212! 1213 DO 110 J = 1,NRHFT 1214 DO 120 B = 1,NVIRT 1215! 1216 NBJ = NVIRT*(J-1) + B 1217! 1218 DO 130 NAI = 1,NT1AMX 1219! 1220 AIBJCK = 0.0D0 1221 DO 140 D = 1,NVIRT 1222! 1223 NDJ = NVIRT*(J-1) + D 1224! 1225 AIBJCK = AIBJCK + XINT1(NCK,B,D)*T2AM(NDJ,NAI) 1226! 1227 140 CONTINUE 1228! 1229 DO 150 L = 1,NRHFT 1230! 1231 NBL = NVIRT*(L-1) + B 1232! 1233 AIBJCK = AIBJCK - XINT2(NCK,L,J)*T2AM(NBL,NAI) 1234! 1235 150 CONTINUE 1236! 1237 AIBJCK = AIBJCK/(SCR1(NAI) + SCR1(NBJ) + SCR1(NCK)) 1238! 1239 T3AM(NAI,NBJ,NCK) = T3AM(NAI,NBJ,NCK) - AIBJCK 1240 T3AM(NAI,NCK,NBJ) = T3AM(NAI,NCK,NBJ) - AIBJCK 1241 T3AM(NBJ,NAI,NCK) = T3AM(NBJ,NAI,NCK) - AIBJCK 1242 T3AM(NCK,NAI,NBJ) = T3AM(NCK,NAI,NBJ) - AIBJCK 1243 T3AM(NBJ,NCK,NAI) = T3AM(NBJ,NCK,NAI) - AIBJCK 1244 T3AM(NCK,NBJ,NAI) = T3AM(NCK,NBJ,NAI) - AIBJCK 1245! 1246 130 CONTINUE 1247 120 CONTINUE 1248 110 CONTINUE 1249 100 CONTINUE 1250! 1251 RETURN 1252 END 1253C /* Deck cc3exci_int2 */ 1254 SUBROUTINE CC3EXCI_INT2(XINT,XLAMDP,XLAMDH,AOINT,IDEL, 1255 * XINTTEMP,XINTTEMP2) 1256! 1257! Written by Kasper Hald, April 2000 1258! 1259 IMPLICIT NONE 1260! 1261#include "priunit.h" 1262#include "ccorb.h" 1263#include "ccsdinp.h" 1264#include "ccsdsym.h" 1265! 1266 INTEGER IB, INDEX, IDEL, NAB 1267! 1268#if defined (SYS_CRAY) 1269 REAL XINT(NORBT,NORBT,NORBT,NORBT) 1270 REAL XLAMDP(NBAST,NORBT), XLAMDH(NBAST,NORBT) 1271 REAL AOINT((NBAST*(NBAST+1))/2,NBAST), ZERO 1272 REAL XINTTEMP(NBAST,NBAST,NBAST,NBAST) 1273 REAL XINTTEMP2(NBAST,NBAST,NBAST,NBAST) 1274#else 1275 DOUBLE PRECISION XINT(NORBT,NORBT,NORBT,NORBT) 1276 DOUBLE PRECISION XLAMDP(NBAST,NORBT), XLAMDH(NBAST,NORBT) 1277 DOUBLE PRECISION AOINT((NBAST*(NBAST+1))/2,NBAST), ZERO 1278 DOUBLE PRECISION XINTTEMP(NBAST,NBAST,NBAST,NBAST) 1279 DOUBLE PRECISION XINTTEMP2(NBAST,NBAST,NBAST,NBAST) 1280#endif 1281! 1282 PARAMETER (ZERO = 0.0D0) 1283! 1284! 1285 INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 1286! 1287 CALL DZERO(XINTTEMP,NBAST*NBAST*NBAST*NBAST) 1288 CALL DZERO(XINTTEMP2,NBAST*NBAST*NBAST*NBAST) 1289! 1290!============================================ 1291! Calculate the needed integrals. 1292!============================================ 1293! 1294 DO G = 1,NBAST 1295 DO IB = 1,NBAST 1296 DO A = 1,NBAST 1297 NAB = INDEX(A,IB) 1298! 1299 IF (.NOT. (AOINT(NAB,G) .EQ. ZERO)) THEN 1300! 1301 DO P = 1, NORBT 1302 XINTTEMP(P,IB,G,IDEL) = XINTTEMP(P,IB,G,IDEL) 1303 * + AOINT(NAB,G)*XLAMDP(A,P) 1304 ENDDO 1305 ENDIF 1306 ENDDO 1307 ENDDO 1308 ENDDO 1309! 1310 DO P = 1, NORBT 1311 DO G = 1, NBAST 1312 DO IB = 1, NBAST 1313 IF (.NOT. (XINTTEMP(P,IB,G,IDEL) .EQ. ZERO)) THEN 1314 DO Q = 1, NORBT 1315 XINTTEMP2(P,Q,G,IDEL) = XINTTEMP2(P,Q,G,IDEL) 1316 * + XINTTEMP(P,IB,G,IDEL)*XLAMDH(IB,Q) 1317 ENDDO 1318 ENDIF 1319 ENDDO 1320 ENDDO 1321 ENDDO 1322! 1323 CALL DZERO(XINTTEMP,NBAST*NBAST*NBAST*NBAST) 1324! 1325 DO P = 1, NORBT 1326 DO Q = 1, NORBT 1327 DO G = 1, NBAST 1328 IF (.NOT. (XINTTEMP2(P,Q,G,IDEL) .EQ. ZERO)) THEN 1329 DO R = 1, NORBT 1330 XINTTEMP(P,Q,R,IDEL) = XINTTEMP(P,Q,R,IDEL) 1331 * + XINTTEMP2(P,Q,G,IDEL)*XLAMDP(G,R) 1332 ENDDO 1333 ENDIF 1334 ENDDO 1335 ENDDO 1336 ENDDO 1337! 1338 DO P = 1, NORBT 1339 DO Q = 1, NORBT 1340 DO R = 1, NORBT 1341 IF (.NOT. (XINTTEMP(P,Q,R,IDEL) .EQ. ZERO)) THEN 1342 DO S = 1, NORBT 1343 XINT(P,Q,R,S) = XINT(P,Q,R,S) + 1344 * XINTTEMP(P,Q,R,IDEL)*XLAMDH(IDEL,S) 1345 ENDDO 1346 ENDIF 1347 ENDDO 1348 ENDDO 1349 ENDDO 1350! 1351 RETURN 1352 END 1353C /* Deck cc33_transtemp */ 1354 SUBROUTINE CC33_TRANSTEMP(XINT,C1AM,XINT1S,XINT2S,XINT3S, 1355 * XINT4S,XINT5S,XINT6S,XINT2, 1356 * XIAJB,YIAJB,XINT1T,XINT2T, 1357 * XINT3T,XINT4T,XITRAN1,XITRAN2, 1358 * XITRAN3,XITRAN4) 1359! 1360! Written by Kasper Hald January 2000. 1361! 1362! 1363 IMPLICIT NONE 1364! 1365#include "priunit.h" 1366#include "maxorb.h" 1367#include "maxash.h" 1368#include "mxcent.h" 1369#include "aovec.h" 1370#include "ccsdsym.h" 1371#include "ccorb.h" 1372! 1373 INTEGER NCK, NCL, NDK, NAI, NBJ, NBK 1374! 1375#if defined (SYS_CRAY) 1376 REAL XINT(NORBT,NORBT,NORBT,NORBT) 1377 REAL XINT2(NORBT,NORBT,NORBT,NORBT) 1378 REAL XITRAN1(NRHFT,NORBT,NORBT,NORBT) 1379 REAL XITRAN2(NVIRT,NORBT,NORBT,NORBT) 1380 REAL XITRAN3(NORBT,NRHFT,NORBT,NORBT) 1381 REAL XITRAN4(NORBT,NVIRT,NORBT,NORBT) 1382 REAL XINT1S(NT1AMX,NVIRT,NVIRT) 1383 REAL XINT2S(NT1AMX,NRHFT,NRHFT) 1384 REAL XINT3S(NT1AMX,NVIRT,NVIRT) 1385 REAL XINT4S(NT1AMX,NRHFT,NRHFT) 1386 REAL XINT5S(NT1AMX,NVIRT,NVIRT) 1387 REAL XINT6S(NT1AMX,NRHFT,NRHFT) 1388 REAL XINT1T(NT1AMX,NVIRT,NVIRT) 1389 REAL XINT2T(NT1AMX,NRHFT,NRHFT) 1390 REAL XINT3T(NT1AMX,NVIRT,NVIRT) 1391 REAL XINT4T(NT1AMX,NRHFT,NRHFT) 1392 REAL TEMP 1393 REAL C1AM(NT1AMX), TWO 1394 REAL XIAJB(NT1AMX,NT1AMX), YIAJB(NT1AMX,NT1AMX) 1395#else 1396 DOUBLE PRECISION XINT(NORBT,NORBT,NORBT,NORBT) 1397 DOUBLE PRECISION XINT2(NORBT,NORBT,NORBT,NORBT) 1398 DOUBLE PRECISION XITRAN1(NRHFT,NORBT,NORBT,NORBT) 1399 DOUBLE PRECISION XITRAN2(NVIRT,NORBT,NORBT,NORBT) 1400 DOUBLE PRECISION XITRAN3(NORBT,NRHFT,NORBT,NORBT) 1401 DOUBLE PRECISION XITRAN4(NORBT,NVIRT,NORBT,NORBT) 1402 DOUBLE PRECISION XINT1S(NT1AMX,NVIRT,NVIRT) 1403 DOUBLE PRECISION XINT2S(NT1AMX,NRHFT,NRHFT) 1404 DOUBLE PRECISION XINT3S(NT1AMX,NVIRT,NVIRT) 1405 DOUBLE PRECISION XINT4S(NT1AMX,NRHFT,NRHFT) 1406 DOUBLE PRECISION XINT5S(NT1AMX,NVIRT,NVIRT) 1407 DOUBLE PRECISION XINT6S(NT1AMX,NRHFT,NRHFT) 1408 DOUBLE PRECISION XINT1T(NT1AMX,NVIRT,NVIRT) 1409 DOUBLE PRECISION XINT2T(NT1AMX,NRHFT,NRHFT) 1410 DOUBLE PRECISION XINT3T(NT1AMX,NVIRT,NVIRT) 1411 DOUBLE PRECISION XINT4T(NT1AMX,NRHFT,NRHFT) 1412 DOUBLE PRECISION TEMP 1413 DOUBLE PRECISION C1AM(NT1AMX), TWO 1414 DOUBLE PRECISION XIAJB(NT1AMX,NT1AMX), YIAJB(NT1AMX,NT1AMX) 1415#endif 1416! 1417 LOGICAL LOCDBG 1418! 1419 PARAMETER (TWO = 2.0D0) 1420 PARAMETER (LOCDBG = .FALSE.) 1421! 1422! 1423 DO A = 1, NVIRT 1424 DO I = 1, NRHFT 1425 DO B = 1, NVIRT 1426 DO C = 1, NVIRT 1427 NAI = NVIRT*(I-1)+A 1428 XINT1S(NAI,B,C) = XINT1S(NAI,B,C) + 1429 * XINT(NRHFT+A,I,NRHFT+B,NRHFT+C) 1430 XINT1T(NAI,B,C) = XINT1T(NAI,B,C) + 1431 * XINT(I,NRHFT+A,NRHFT+B,NRHFT+C) 1432 ENDDO 1433 ENDDO 1434 ENDDO 1435 ENDDO 1436! 1437 DO I = 1, NRHFT 1438 DO A = 1, NVIRT 1439 DO J = 1, NRHFT 1440 DO K = 1, NRHFT 1441 NAI = NVIRT*(I-1)+A 1442 XINT2S(NAI,J,K) = XINT2S(NAI,J,K) + 1443 * XINT(NRHFT+A,I,J,K) 1444 XINT2T(NAI,J,K) = XINT2T(NAI,J,K) + 1445 * XINT(I,NRHFT+A,J,K) 1446 ENDDO 1447 ENDDO 1448 ENDDO 1449 ENDDO 1450! 1451 DO P = 1, NORBT 1452! 1453 DO Q = 1, NORBT 1454! 1455 DO R = 1, NORBT 1456! 1457 DO K = 1, NRHFT 1458! 1459 TEMP = 0.0D0 1460 DO D = 1, NVIRT 1461! 1462 NDK = NVIRT*(K-1) + D 1463 TEMP = TEMP + C1AM(NDK)*XINT(P,D+NRHFT,Q,R) 1464! 1465 ENDDO ! The D-loop stops here. 1466! 1467 XITRAN3(P,K,Q,R) = TEMP 1468! 1469 ENDDO ! The K-loop stops here. 1470! 1471 DO C = 1, NVIRT 1472! 1473 TEMP = 0.0D0 1474 DO L = 1, NRHFT 1475! 1476 NCL = NVIRT*(L-1) + C 1477 TEMP = TEMP + C1AM(NCL)*XINT(L,P,Q,R) 1478! 1479 ENDDO ! The L-loop stops here 1480! 1481 XITRAN2(C,P,Q,R) = -TEMP 1482! 1483 ENDDO ! The C-loop ends here. 1484! 1485 ENDDO ! The Q-loop stops here. 1486 ENDDO ! The R-loop stops here. 1487 ENDDO ! The S-loop stops here. 1488! 1489!==================================== 1490! Debug information. 1491!==================================== 1492! 1493 IF (LOCDBG) THEN 1494 DO A = 1, NVIRT 1495 DO K = 1, NORBT 1496 DO B = 1, NORBT 1497 DO C = 1, NORBT 1498 IF (ABS(XITRAN2(A,K,B,C)) .GT. 1.0D-8) THEN 1499 WRITE(LUPRI,*)'C, P, R, S :',A,K,B,C 1500 WRITE(LUPRI,*)'XITRAN2(C,P,R,S) = ', 1501 * XITRAN2(A,K,B,C) 1502 ENDIF 1503 ENDDO 1504 ENDDO 1505 ENDDO 1506 ENDDO 1507 WRITE(LUPRI,*)' ' 1508! 1509 DO A = 1, NORBT 1510 DO K = 1, NRHFT 1511 DO B = 1, NORBT 1512 DO C = 1, NORBT 1513 IF (ABS(XITRAN3(A,K,B,C)) .GT. 1.0D-8) THEN 1514 WRITE(LUPRI,*)'P, K, R, S :',A,K,B,C 1515 WRITE(LUPRI,*)'XITRAN3(P,K,R,S) = ', 1516 * XITRAN3(A,K,B,C) 1517 ENDIF 1518 ENDDO 1519 ENDDO 1520 ENDDO 1521 ENDDO 1522 WRITE(LUPRI,*)' ' 1523 ENDIF 1524! 1525!======================= 1526! 1527 DO A = 1, NVIRT 1528 DO B = 1, NVIRT 1529 DO C = 1, NVIRT 1530 DO K = 1, NRHFT 1531! 1532 NCK = NVIRT*(K-1) + C 1533 XINT3S(NCK,A,B) = XINT3S(NCK,A,B) 1534 * + XITRAN2(C,B+NRHFT,A+NRHFT,K) 1535 * - XITRAN2(C,K,A+NRHFT,B+NRHFT) 1536 * - XITRAN3(C+NRHFT,K,A+NRHFT,B+NRHFT) 1537 XINT5S(NCK,A,B) = XINT5S(NCK,A,B) 1538 * + XITRAN2(A,B+NRHFT,C+NRHFT,K) 1539 * - XITRAN2(C,K,A+NRHFT,B+NRHFT) 1540 * - XITRAN3(C+NRHFT,K,A+NRHFT,B+NRHFT) 1541! 1542 ENDDO ! The K-loop ends here 1543 ENDDO ! The C-loop ends here. 1544 ENDDO ! The B-loop ends here. 1545 ENDDO ! The A-loop ends here. 1546! 1547 DO I = 1, NRHFT 1548 DO J = 1, NRHFT 1549 DO C = 1, NVIRT 1550 DO K = 1, NRHFT 1551! 1552 NCK = NVIRT*(K-1) + C 1553 XINT4S(NCK,I,J) = XINT4S(NCK,I,J) 1554 * + XITRAN3(I,K,C+NRHFT,J) 1555 * - XITRAN3(C+NRHFT,K,I,J) 1556 * - XITRAN2(C,K,I,J) 1557! 1558 XINT6S(NCK,I,J) = XINT6S(NCK,I,J) 1559 * + XITRAN3(I,J,C+NRHFT,K) 1560 * - XITRAN3(C+NRHFT,K,I,J) 1561 * - XITRAN2(C,K,I,J) 1562! 1563 ENDDO ! The K-loop ends here 1564 ENDDO ! The C-loop ends here 1565 ENDDO ! The J-loop ends here 1566 ENDDO ! The I-loop ends here 1567! 1568!============================== 1569! More debug information. 1570!============================== 1571! 1572 IF (LOCDBG) THEN 1573 DO A = 1, NT1AMX 1574 DO I = 1, NRHFT 1575 DO J = 1, NRHFT 1576 IF (ABS(XINT3S(A,I,J)) .GT. 1.0D-8) THEN 1577 WRITE(LUPRI,*)'XINT3S(NAI,J,K) ',XINT3S(A,I,J) 1578 WRITE(LUPRI,*)'NAI, J, K ',A,I,J 1579 ENDIF 1580 IF (ABS(XINT5S(A,I,J)) .GT. 1.0D-8) THEN 1581 WRITE(LUPRI,*)'XINT5S(NAI,J,K) ',XINT5S(A,I,J) 1582 WRITE(LUPRI,*)'NAI, J, K ',A,I,J 1583 ENDIF 1584 ENDDO 1585 ENDDO 1586 ENDDO 1587! 1588 DO A = 1, NT1AMX 1589 DO C = 1, NVIRT 1590 DO D = 1, NVIRT 1591 IF (ABS(XINT4S(A,C,D)) .GT. 1.0D-8) THEN 1592 WRITE(LUPRI,*)'XINT4S(NAI,B,C) ',XINT4S(A,C,D) 1593 WRITE(LUPRI,*)'NAI, B, C ',A,C,D 1594 ENDIF 1595 IF (ABS(XINT6S(A,C,D)) .GT. 1.0D-8) THEN 1596 WRITE(LUPRI,*)'XINT6S(NAI,B,C) ',XINT6S(A,C,D) 1597 WRITE(LUPRI,*)'NAI, B, C ',A,C,D 1598 ENDIF 1599 ENDDO 1600 ENDDO 1601 ENDDO 1602 ENDIF 1603! 1604!=========================== 1605! 1606 DO A = 1, NVIRT 1607 DO I = 1, NRHFT 1608 NAI = NVIRT*(I-1) + A 1609! 1610 DO B = 1, NVIRT 1611 DO J = 1, NRHFT 1612 NBJ = NVIRT*(J-1) + B 1613! 1614 XIAJB(NAI,NBJ) = TWO*XINT2(A+NRHFT,I,B+NRHFT,J) 1615 * - XINT2(A+NRHFT,J,B+NRHFT,I) 1616 YIAJB(NAI,NBJ) = XINT2(A+NRHFT,I,B+NRHFT,J) 1617! 1618 ENDDO 1619 ENDDO 1620 ENDDO 1621 ENDDO 1622! 1623! 1624 DO A = 1, NVIRT 1625 DO I = 1, NRHFT 1626 NAI = NVIRT*(I-1) + A 1627 DO B = 1, NVIRT 1628 DO C = 1, NVIRT 1629 TEMP = 0.0D0 1630 DO K = 1, NRHFT 1631 NCK = NVIRT*(K-1) + C 1632 TEMP = TEMP + C1AM(NCK)*XINT(I,A+NRHFT,K,B+NRHFT) 1633 ENDDO 1634 XINT3T(NAI,C,B) = -TEMP 1635 ENDDO 1636 ENDDO 1637 ENDDO 1638 ENDDO 1639! 1640! 1641 DO A = 1, NVIRT 1642 DO I = 1, NRHFT 1643 NAI = NVIRT*(I-1) + A 1644 DO J = 1, NRHFT 1645 DO K = 1, NRHFT 1646 TEMP = 0.0D0 1647 DO B = 1, NVIRT 1648! 1649 NBK = NVIRT*(K-1) + B 1650 TEMP = TEMP + C1AM(NBK)*XINT(I,A+NRHFT,J,B+NRHFT) 1651! 1652 ENDDO 1653! 1654 XINT4T(NAI,J,K) = TEMP 1655! 1656 ENDDO 1657 ENDDO 1658 ENDDO 1659 ENDDO 1660! 1661! 1662 RETURN 1663 END 1664C /* Deck intcompare */ 1665 SUBROUTINE INTCOMPARE(XIAJB,XIAJB2,XIAJB3,YIAJB,YIAJB2,YIAJB3) 1666! 1667! Written by K. Hald, Jan. 2000. 1668! 1669! 1670 IMPLICIT NONE 1671! 1672#include "priunit.h" 1673#include "ccsdsym.h" 1674! 1675#if defined (SYS_CRAY) 1676 REAL XIAJB(NT1AMX,NT1AMX), XIAJB2(NT1AMX,NT1AMX) 1677 REAL YIAJB(NT1AMX,NT1AMX), YIAJB2(NT1AMX,NT1AMX) 1678 REAL XIAJB3(NT1AMX,NT1AMX), YIAJB3(NT1AMX,NT1AMX) 1679#else 1680 DOUBLE PRECISION XIAJB(NT1AMX,NT1AMX), XIAJB2(NT1AMX,NT1AMX) 1681 DOUBLE PRECISION YIAJB(NT1AMX,NT1AMX), YIAJB2(NT1AMX,NT1AMX) 1682 DOUBLE PRECISION XIAJB3(NT1AMX,NT1AMX), YIAJB3(NT1AMX,NT1AMX) 1683#endif 1684! 1685 DO A = 1, NT1AMX 1686 DO I = 1, NT1AMX 1687 IF (ABS(XIAJB(A,I)) .GT. 1.0D-8) THEN 1688 WRITE(LUPRI,*)'IAJB (',A,',',I,') = ', 1689 * XIAJB(A,I) 1690 ENDIF 1691 IF (ABS(XIAJB2(A,I)) .GT. 1.0D-8) THEN 1692 WRITE(LUPRI,*)'IAJBTEMP (',A,',',I,') = ', 1693 * XIAJB2(A,I) 1694 ENDIF 1695 IF (ABS(XIAJB3(A,I)) .GT. 1.0D-8) THEN 1696 WRITE(LUPRI,*)'IAJBTEMP(CMO)(',A,',',I,') = ', 1697 * XIAJB3(A,I) 1698 ENDIF 1699 IF ((ABS(XIAJB(A,I)-XIAJB2(A,I)) .GT. 1.0D-8)) THEN 1700 WRITE(LUPRI,*)'WARNING THE TWO FIRST VALUES DIFFER!!' 1701 ENDIF 1702 IF ((ABS(XIAJB(A,I)-XIAJB3(A,I)) .GT. 1.0D-8)) THEN 1703 WRITE(LUPRI,*)'WARNING THE 1.ST AND 3.RD VALUE DIFFER!' 1704 ENDIF 1705 IF ((ABS(XIAJB2(A,I)-XIAJB3(A,I)) .GT. 1.0D-8)) THEN 1706 WRITE(LUPRI,*)'WARNING THE 2.ND AND 3.RD VALUE DIFFER!' 1707 ENDIF 1708 ENDDO 1709 ENDDO 1710! 1711 DO A = 1, NT1AMX 1712 DO I = 1, NT1AMX 1713 IF (ABS(YIAJB(A,I)) .GT. 1.0D-8) THEN 1714 WRITE(LUPRI,*)'YIAJB (',A,',',I,') = ', 1715 * YIAJB(A,I) 1716 ENDIF 1717 IF (ABS(YIAJB2(A,I)) .GT. 1.0D-8) THEN 1718 WRITE(LUPRI,*)'YIAJBTEMP (',A,',',I,') = ', 1719 * YIAJB2(A,I) 1720 ENDIF 1721 IF (ABS(YIAJB3(A,I)) .GT. 1.0D-8) THEN 1722 WRITE(LUPRI,*)'YIAJBTEMP(CMO)(',A,',',I,') = ', 1723 * YIAJB3(A,I) 1724 ENDIF 1725! 1726 CALL FLSHFO(LUPRI) 1727! 1728 IF (ABS(YIAJB(A,I)-YIAJB2(A,I)) .GT. 1.0D-8) THEN 1729 WRITE(LUPRI,*)'WARNING 1. AND 2. VALUES DIFFER!!' 1730 ENDIF 1731 IF (ABS(YIAJB(A,I)-YIAJB3(A,I)) .GT. 1.0D-8) THEN 1732 WRITE(LUPRI,*)'WARNING 1. AND 3. VALUES DIFFER!!' 1733 ENDIF 1734 IF (ABS(YIAJB2(A,I)-YIAJB3(A,I)) .GT. 1.0D-8) THEN 1735 WRITE(LUPRI,*)'WARNING 2. AND 3. VALUES DIFFER!!' 1736 ENDIF 1737 ENDDO 1738 ENDDO 1739! 1740 RETURN 1741 END 1742C /* Deck gcompare */ 1743 SUBROUTINE GCOMPARE(INT1,INT2) 1744! 1745! Written by K. Hald, April 2000. 1746! 1747! 1748 IMPLICIT NONE 1749! 1750#include "priunit.h" 1751#include "ccsdsym.h" 1752#include "ccorb.h" 1753! 1754#if defined (SYS_CRAY) 1755 REAL INT1(NORBT,NORBT,NORBT,NORBT) 1756 REAL INT2(NORBT,NORBT,NORBT,NORBT) 1757#else 1758 DOUBLE PRECISION INT1(NORBT,NORBT,NORBT,NORBT) 1759 DOUBLE PRECISION INT2(NORBT,NORBT,NORBT,NORBT) 1760#endif 1761! 1762 DO A = 1, NORBT 1763 DO I = 1, NORBT 1764 DO B = 1, NORBT 1765 DO J = 1, NORBT 1766 IF (ABS(INT1(A,I,B,J)) .GT. 1.0D-8) THEN 1767 WRITE(LUPRI,*)'G(lambda) = ',INT1(A,I,B,J) 1768 ENDIF 1769 IF (ABS(INT2(A,I,B,J)) .GT. 1.0D-8) THEN 1770 WRITE(LUPRI,*)'G(CMO) = ',INT2(A,I,B,J) 1771 ENDIF 1772 IF (ABS(INT1(A,I,B,J)-INT2(A,I,B,J)) .GT. 1.0D-8) THEN 1773 WRITE(LUPRI,*)'THE TWO VALUES DIFFER' 1774 ENDIF 1775 ENDDO 1776 ENDDO 1777 ENDDO 1778 ENDDO 1779! 1780 RETURN 1781 END 1782C /* Deck intcompare2 */ 1783 SUBROUTINE INTCOMPARE2(XINT1S,XINT1STEMP,XINT2S,XINT2STEMP, 1784 * XINT3S,XINT3STEMP,XINT4S,XINT4STEMP, 1785 * XINT5S,XINT5STEMP,XINT6S,XINT6STEMP) 1786! 1787! Written by K. Hald, April 2000. 1788! 1789! 1790 IMPLICIT NONE 1791! 1792#include "priunit.h" 1793#include "ccsdsym.h" 1794#include "ccorb.h" 1795! 1796#if defined (SYS_CRAY) 1797 REAL XINT1S(NT1AMX,NVIRT,NVIRT) 1798 REAL XINT1STEMP(NT1AMX,NVIRT,NVIRT) 1799 REAL XINT2S(NT1AMX,NRHFT,NRHFT) 1800 REAL XINT2STEMP(NT1AMX,NRHFT,NRHFT) 1801 REAL XINT3S(NT1AMX,NVIRT,NVIRT) 1802 REAL XINT3STEMP(NT1AMX,NVIRT,NVIRT) 1803 REAL XINT4S(NT1AMX,NRHFT,NRHFT) 1804 REAL XINT4STEMP(NT1AMX,NRHFT,NRHFT) 1805 REAL XINT5S(NT1AMX,NVIRT,NVIRT) 1806 REAL XINT5STEMP(NT1AMX,NVIRT,NVIRT) 1807 REAL XINT6S(NT1AMX,NRHFT,NRHFT) 1808 REAL XINT6STEMP(NT1AMX,NRHFT,NRHFT) 1809#else 1810 DOUBLE PRECISION XINT1S(NT1AMX,NVIRT,NVIRT) 1811 DOUBLE PRECISION XINT1STEMP(NT1AMX,NVIRT,NVIRT) 1812 DOUBLE PRECISION XINT2S(NT1AMX,NRHFT,NRHFT) 1813 DOUBLE PRECISION XINT2STEMP(NT1AMX,NRHFT,NRHFT) 1814 DOUBLE PRECISION XINT3S(NT1AMX,NVIRT,NVIRT) 1815 DOUBLE PRECISION XINT3STEMP(NT1AMX,NVIRT,NVIRT) 1816 DOUBLE PRECISION XINT4S(NT1AMX,NRHFT,NRHFT) 1817 DOUBLE PRECISION XINT4STEMP(NT1AMX,NRHFT,NRHFT) 1818 DOUBLE PRECISION XINT5S(NT1AMX,NVIRT,NVIRT) 1819 DOUBLE PRECISION XINT5STEMP(NT1AMX,NVIRT,NVIRT) 1820 DOUBLE PRECISION XINT6S(NT1AMX,NRHFT,NRHFT) 1821 DOUBLE PRECISION XINT6STEMP(NT1AMX,NRHFT,NRHFT) 1822#endif 1823! 1824 DO A = 1, NT1AMX 1825 DO B = 1, NVIRT 1826 DO I = 1, NVIRT 1827 IF (ABS(XINT1S(A,B,I)) .GT. 1.0D-8) THEN 1828 WRITE(LUPRI,*)'XINT1S = ',XINT1S(A,B,I) 1829 ENDIF 1830 IF (ABS(XINT1STEMP(A,B,I)) .GT. 1.0D-8) THEN 1831 WRITE(LUPRI,*)'XINT1STEMP = ',XINT1STEMP(A,B,I) 1832 ENDIF 1833 IF (ABS(XINT1S(A,B,I)-XINT1STEMP(A,B,I)) .GT. 1.0D-8) THEN 1834 WRITE(LUPRI,*)'For AI, B, C',A,B,I 1835 WRITE(LUPRI,*)'XINT1S & XINT1STEMP',XINT1S(A,B,I), 1836 * XINT1STEMP(A,B,I) 1837 WRITE(LUPRI,*)'XINT1S AND XINT1STEMP DIFFER' 1838 ENDIF 1839 ENDDO 1840 ENDDO 1841 ENDDO 1842! 1843 DO A = 1, NT1AMX 1844 DO B = 1, NRHFT 1845 DO I = 1, NRHFT 1846 IF (ABS(XINT2S(A,B,I)) .GT. 1.0D-8) THEN 1847 WRITE(LUPRI,*)'XINT2S = ',XINT2S(A,B,I) 1848 ENDIF 1849 IF (ABS(XINT2STEMP(A,B,I)) .GT. 1.0D-8) THEN 1850 WRITE(LUPRI,*)'XINT2STEMP = ',XINT2STEMP(A,B,I) 1851 ENDIF 1852 IF (ABS(XINT2S(A,B,I)-XINT2STEMP(A,B,I)) .GT. 1.0D-8) 1853 * THEN 1854 WRITE(LUPRI,*)'For AI, B, C',A,B,I 1855 WRITE(LUPRI,*)'XINT2S & XINT2STEMP',XINT2S(A,B,I), 1856 * XINT2STEMP(A,B,I) 1857 WRITE(LUPRI,*)'XINT2S AND XINT2STEMP DIFFER' 1858 ENDIF 1859 ENDDO 1860 ENDDO 1861 ENDDO 1862! 1863 DO A = 1, NT1AMX 1864 DO B = 1, NVIRT 1865 DO I = 1, NVIRT 1866 IF (ABS(XINT3S(A,B,I)) .GT. 1.0D-8) THEN 1867 WRITE(LUPRI,*)'XINT3S = ',XINT3S(A,B,I) 1868 ENDIF 1869 IF (ABS(XINT3STEMP(A,B,I)) .GT. 1.0D-8) THEN 1870 WRITE(LUPRI,*)'XINT3STEMP = ',XINT3STEMP(A,B,I) 1871 ENDIF 1872 IF (ABS(XINT3S(A,B,I)-XINT3STEMP(A,B,I)) .GT. 1.0D-8) 1873 * THEN 1874 WRITE(LUPRI,*)'For AI, B, C',A,B,I 1875 WRITE(LUPRI,*)'XINT3S & XINT3STEMP',XINT3S(A,B,I), 1876 * XINT3STEMP(A,B,I) 1877 WRITE(LUPRI,*)'XINT3S AND XINT3STEMP DIFFER' 1878 ENDIF 1879 ENDDO 1880 ENDDO 1881 ENDDO 1882! 1883 DO A = 1, NT1AMX 1884 DO B = 1, NRHFT 1885 DO I = 1, NRHFT 1886 IF (ABS(XINT4S(A,B,I)) .GT. 1.0D-8) THEN 1887 WRITE(LUPRI,*)'XINT4S = ',XINT4S(A,B,I) 1888 ENDIF 1889 IF (ABS(XINT4STEMP(A,B,I)) .GT. 1.0D-8) THEN 1890 WRITE(LUPRI,*)'XINT4STEMP = ',XINT4STEMP(A,B,I) 1891 ENDIF 1892 IF (ABS(XINT4S(A,B,I)-XINT4STEMP(A,B,I)) .GT. 1.0D-8) 1893 * THEN 1894 WRITE(LUPRI,*)'For AI, B, C',A,B,I 1895 WRITE(LUPRI,*)'XINT4S & XINT4STEMP',XINT4S(A,B,I), 1896 * XINT4STEMP(A,B,I) 1897 WRITE(LUPRI,*)'XINT4S AND XINT4STEMP DIFFER' 1898 ENDIF 1899 ENDDO 1900 ENDDO 1901 ENDDO 1902! 1903 DO A = 1, NT1AMX 1904 DO B = 1, NVIRT 1905 DO I = 1, NVIRT 1906 IF (ABS(XINT5S(A,B,I)) .GT. 1.0D-8) THEN 1907 WRITE(LUPRI,*)'XINT5S = ',XINT5S(A,B,I) 1908 ENDIF 1909 IF (ABS(XINT5STEMP(A,B,I)) .GT. 1.0D-8) THEN 1910 WRITE(LUPRI,*)'XINT5STEMP = ',XINT5STEMP(A,B,I) 1911 ENDIF 1912 IF (ABS(XINT5S(A,B,I)-XINT5STEMP(A,B,I)) .GT. 1.0D-8) 1913 * THEN 1914 WRITE(LUPRI,*)'For AI, B, C',A,B,I 1915 WRITE(LUPRI,*)'XINT5S & XINT5STEMP',XINT5S(A,B,I), 1916 * XINT5STEMP(A,B,I) 1917 WRITE(LUPRI,*)'XINT5S AND XINT5STEMP DIFFER' 1918 ENDIF 1919 ENDDO 1920 ENDDO 1921 ENDDO 1922! 1923 DO A = 1, NT1AMX 1924 DO B = 1, NRHFT 1925 DO I = 1, NRHFT 1926 IF (ABS(XINT6S(A,B,I)) .GT. 1.0D-8) THEN 1927 WRITE(LUPRI,*)'XINT6S = ',XINT6S(A,B,I) 1928 ENDIF 1929 IF (ABS(XINT6STEMP(A,B,I)) .GT. 1.0D-8) THEN 1930 WRITE(LUPRI,*)'XINT6STEMP = ',XINT6STEMP(A,B,I) 1931 ENDIF 1932 IF (ABS(XINT6S(A,B,I)-XINT6STEMP(A,B,I)) .GT. 1.0D-8) 1933 * THEN 1934 WRITE(LUPRI,*)'For AI, B, C',A,B,I 1935 WRITE(LUPRI,*)'XINT6S & XINT6STEMP',XINT6S(A,B,I), 1936 * XINT6STEMP(A,B,I) 1937 WRITE(LUPRI,*)'XINT6S AND XINT6STEMP DIFFER' 1938 ENDIF 1939 ENDDO 1940 ENDDO 1941 ENDDO 1942! 1943 RETURN 1944 END 1945C /* Deck tcompare */ 1946 SUBROUTINE TCOMPARE(XINT1T,XINT1TTEMP,XINT2T,XINT2TTEMP, 1947 * XINT3T,XINT3TTEMP,XINT4T,XINT4TTEMP) 1948! 1949! Written by K. Hald, April 2000. 1950! 1951! 1952 IMPLICIT NONE 1953! 1954#include "priunit.h" 1955#include "ccsdsym.h" 1956#include "ccorb.h" 1957! 1958#if defined (SYS_CRAY) 1959 REAL XINT1T(NT1AMX,NVIRT,NVIRT) 1960 REAL XINT1TTEMP(NT1AMX,NVIRT,NVIRT) 1961 REAL XINT2T(NT1AMX,NRHFT,NRHFT) 1962 REAL XINT2TTEMP(NT1AMX,NRHFT,NRHFT) 1963 REAL XINT3T(NT1AMX,NVIRT,NVIRT) 1964 REAL XINT3TTEMP(NT1AMX,NVIRT,NVIRT) 1965 REAL XINT4T(NT1AMX,NRHFT,NRHFT) 1966 REAL XINT4TTEMP(NT1AMX,NRHFT,NRHFT) 1967#else 1968 DOUBLE PRECISION XINT1T(NT1AMX,NVIRT,NVIRT) 1969 DOUBLE PRECISION XINT1TTEMP(NT1AMX,NVIRT,NVIRT) 1970 DOUBLE PRECISION XINT2T(NT1AMX,NRHFT,NRHFT) 1971 DOUBLE PRECISION XINT2TTEMP(NT1AMX,NRHFT,NRHFT) 1972 DOUBLE PRECISION XINT3T(NT1AMX,NVIRT,NVIRT) 1973 DOUBLE PRECISION XINT3TTEMP(NT1AMX,NVIRT,NVIRT) 1974 DOUBLE PRECISION XINT4T(NT1AMX,NRHFT,NRHFT) 1975 DOUBLE PRECISION XINT4TTEMP(NT1AMX,NRHFT,NRHFT) 1976#endif 1977! 1978 DO A = 1, NT1AMX 1979 DO B = 1, NVIRT 1980 DO D = 1, NVIRT 1981 IF (ABS(XINT1T(A,B,D)) .GT. 1.0D-8) THEN 1982 WRITE(LUPRI,*)'XINT1T = ',XINT1T(A,B,D) 1983 ENDIF 1984 IF (ABS(XINT1TTEMP(A,B,D)) .GT. 1.0D-8) THEN 1985 WRITE(LUPRI,*)'XINT1TTEMP = ',XINT1TTEMP(A,B,D) 1986 ENDIF 1987 IF (ABS(XINT1T(A,B,D)-XINT1TTEMP(A,B,D)) .GT. 1.0D-8) THEN 1988 WRITE(LUPRI,*)'For NAI, B, D',A,B,D 1989 WRITE(LUPRI,*)'XINT1T & XINT1TTEMP DIFFER!!!' 1990 ENDIF 1991 IF (ABS(XINT3T(A,B,D)) .GT. 1.0D-8) THEN 1992 WRITE(LUPRI,*)'XINT3T = ',XINT3T(A,B,D) 1993 ENDIF 1994 IF (ABS(XINT3TTEMP(A,B,D)) .GT. 1.0D-8) THEN 1995 WRITE(LUPRI,*)'XINT3TTEMP = ',XINT3TTEMP(A,B,D) 1996 ENDIF 1997 IF (ABS(XINT3T(A,B,D)-XINT3TTEMP(A,B,D)) .GT. 1.0D-8) THEN 1998 WRITE(LUPRI,*)'For NAI, B, D',A,B,D 1999 WRITE(LUPRI,*)'XINT3T = ',XINT3T(A,B,D) 2000 WRITE(LUPRI,*)'XINT3TTEMP = ',XINT3TTEMP(A,B,D) 2001 WRITE(LUPRI,*)'XINT3T & XINT3TTEMP DIFFER!!!' 2002 ENDIF 2003 ENDDO 2004 ENDDO 2005 ENDDO 2006! 2007 DO A = 1, NT1AMX 2008 DO B = 1, NRHFT 2009 DO D = 1, NRHFT 2010 IF (ABS(XINT2T(A,B,D)) .GT. 1.0D-8) THEN 2011 WRITE(LUPRI,*)'XINT2T = ',XINT2T(A,B,D) 2012 ENDIF 2013 IF (ABS(XINT2TTEMP(A,B,D)) .GT. 1.0D-8) THEN 2014 WRITE(LUPRI,*)'XINT2TTEMP = ',XINT2TTEMP(A,B,D) 2015 ENDIF 2016 IF (ABS(XINT2T(A,B,D)-XINT2TTEMP(A,B,D)) .GT. 1.0D-8) THEN 2017 WRITE(LUPRI,*)'For NAI, B, D',A,B,D 2018 WRITE(LUPRI,*)'XINT2T = ',XINT2T(A,B,D) 2019 WRITE(LUPRI,*)'XINT2TTEMP = ',XINT2TTEMP(A,B,D) 2020 WRITE(LUPRI,*)'XINT2T & XINT2TTEMP DIFFER!!!' 2021 ENDIF 2022 IF (ABS(XINT4T(A,B,D)) .GT. 1.0D-8) THEN 2023 WRITE(LUPRI,*)'XINT4T = ',XINT4T(A,B,D) 2024 ENDIF 2025 IF (ABS(XINT4TTEMP(A,B,D)) .GT. 1.0D-8) THEN 2026 WRITE(LUPRI,*)'XINT4TTEMP = ',XINT4TTEMP(A,B,D) 2027 ENDIF 2028 IF (ABS(XINT4T(A,B,D)-XINT4TTEMP(A,B,D)) .GT. 1.0D-8) THEN 2029 WRITE(LUPRI,*)'For NAI, B, D',A,B,D 2030 WRITE(LUPRI,*)'XINT4T = ',XINT4T(A,B,D) 2031 WRITE(LUPRI,*)'XINT4TTEMP = ',XINT4TTEMP(A,B,D) 2032 WRITE(LUPRI,*)'XINT4T & XINT4TTEMP DIFFER!!!' 2033 ENDIF 2034 ENDDO 2035 ENDDO 2036 ENDDO 2037! 2038 RETURN 2039 END 2040C /* Deck cc3exci_int1 */ 2041 SUBROUTINE CC3EXCI_INT1(CMO,AOINT,IDEL,INTTOT, 2042 * INTTOTTEMP,INTTOTTEMP2) 2043! 2044! Written by K. Hald, April 2000 2045! 2046! Calculate g(p,q,r,s) 2047! 2048 IMPLICIT NONE 2049! 2050#include "priunit.h" 2051#include "ccorb.h" 2052#include "ccsdinp.h" 2053#include "ccsdsym.h" 2054 INTEGER INDEX, NAB, NDL, NCK, IDEL 2055! 2056#if defined (SYS_CRAY) 2057 REAL TWO, AOINT((NBAST*(NBAST+1))/2,NBAST) 2058 REAL CMO(NBAST,NORBT), ZERO 2059 REAL INTTOT(NORBT,NORBT,NORBT,NORBT) 2060 REAL INTTOTTEMP(NBAST,NBAST,NBAST,NBAST) 2061 REAL INTTOTTEMP2(NBAST,NBAST,NBAST,NBAST) 2062#else 2063 DOUBLE PRECISION TWO, AOINT((NBAST*(NBAST+1))/2,NBAST) 2064 DOUBLE PRECISION CMO(NBAST,NORBT), ZERO 2065 DOUBLE PRECISION INTTOT(NORBT,NORBT,NORBT,NORBT) 2066 DOUBLE PRECISION INTTOTTEMP(NBAST,NBAST,NBAST,NBAST) 2067 DOUBLE PRECISION INTTOTTEMP2(NBAST,NBAST,NBAST,NBAST) 2068#endif 2069! 2070 PARAMETER (TWO = 2.0D0, ZERO = 0.0D0) 2071! 2072! 2073 INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 2074! 2075! 2076 CALL DZERO(INTTOTTEMP,NBAST*NBAST*NBAST*NBAST) 2077 CALL DZERO(INTTOTTEMP2,NBAST*NBAST*NBAST*NBAST) 2078! 2079 DO G = 1,NBAST 2080 DO B = 1,NBAST 2081 DO A = 1,NBAST 2082 NAB = INDEX(A,B) 2083! 2084 IF (.NOT.(AOINT(NAB,G) .EQ. ZERO)) THEN 2085 DO C = 1, NORBT 2086 INTTOTTEMP(C,B,G,IDEL) = INTTOTTEMP(C,B,G,IDEL) 2087 * + AOINT(NAB,G)*CMO(A,C) 2088 ENDDO 2089 ENDIF 2090 ENDDO 2091 ENDDO 2092 ENDDO 2093! 2094 DO C = 1, NORBT 2095 DO G = 1,NBAST 2096 DO B = 1,NBAST 2097 IF (.NOT.(INTTOTTEMP(C,B,G,IDEL) .EQ. ZERO)) THEN 2098 DO K = 1, NORBT 2099 INTTOTTEMP2(C,K,G,IDEL) = INTTOTTEMP2(C,K,G,IDEL) 2100 * + INTTOTTEMP(C,B,G,IDEL)*CMO(B,K) 2101 ENDDO 2102 ENDIF 2103 ENDDO 2104 ENDDO 2105 ENDDO 2106! 2107 CALL DZERO(INTTOTTEMP,NBAST*NBAST*NBAST*NBAST) 2108! 2109 DO C = 1, NORBT 2110 DO K = 1, NORBT 2111 DO G = 1, NBAST 2112 IF (.NOT.(INTTOTTEMP2(C,K,G,IDEL) .EQ. ZERO)) THEN 2113 DO D = 1, NORBT 2114 INTTOTTEMP(C,K,D,IDEL) = INTTOTTEMP(C,K,D,IDEL) 2115 * + INTTOTTEMP2(C,K,G,IDEL)*CMO(G,D) 2116 ENDDO 2117 ENDIF 2118 ENDDO 2119 ENDDO 2120 ENDDO 2121! 2122 DO C = 1, NORBT 2123 DO K = 1, NORBT 2124 DO D = 1, NORBT 2125 IF (.NOT. (INTTOTTEMP(C,K,D,IDEL) .EQ. ZERO)) THEN 2126 DO L = 1, NORBT 2127 INTTOT(C,K,D,L) = INTTOT(C,K,D,L) + 2128 * INTTOTTEMP(C,K,D,IDEL)*CMO(IDEL,L) 2129 ENDDO 2130 ENDIF 2131 ENDDO 2132 ENDDO 2133 ENDDO 2134! 2135! 2136! DO C = 1, NORBT 2137! DO K = 1, NORBT 2138! WRITE(LUPRI,*)'Output for C,K = ',C,K 2139! CALL OUTPUT(INTTOT(C,K,1,1),1,NORBT,1,NORBT,NORBT,NORBT,1,6) 2140! ENDDO 2141! ENDDO 2142! 2143! 2144 RETURN 2145 END 2146