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 cc_inifro */ 20 SUBROUTINE CC_INIFRO(WORK,LWORK) 21C 22C Written by Asger Halkier 22/5 - 1998. 23C 24C Set up index arrays needed for frozen core gradients. 25C 26#include "implicit.h" 27C 28 DIMENSION WORK(LWORK) 29C 30 EXTERNAL INIDAT 31C 32#include "priunit.h" 33#include "maxorb.h" 34#include "ccsdinp.h" 35#include "ccorb.h" 36#include "ccsdsym.h" 37#include "ccfro.h" 38#include "symsq.h" 39#include "ccisao.h" 40#include "r12int.h" 41C 42C---------------------------------------- 43C Calculate index arrays for lengths. 44C---------------------------------------- 45C 46 DO 100 ISYMAI = 1,NSYM 47 NT1FRO(ISYMAI) = 0 48 NF1FRO(ISYMAI) = 0 49 NDSFRO(ISYMAI) = 0 50 NCOFRO(ISYMAI) = 0 51 NFROIJ(ISYMAI) = 0 52 NFROFR(ISYMAI) = 0 53 NV1FRO(ISYMAI) = 0 54 DO 110 ISYMI = 1,NSYM 55 ISYMA = MULD2H(ISYMAI,ISYMI) 56 NT1FRO(ISYMAI) = NT1FRO(ISYMAI) + NVIR(ISYMA)*NRHFFR(ISYMI) 57 NF1FRO(ISYMAI) = NF1FRO(ISYMAI) + NRHF(ISYMA)*NRHFFR(ISYMI) 58 NDSFRO(ISYMAI) = NDSFRO(ISYMAI) + NNBST(ISYMA)*NRHFFR(ISYMI) 59 NCOFRO(ISYMAI) = NCOFRO(ISYMAI) + NRHF(ISYMA)*NRHFFR(ISYMI) 60 NFROIJ(ISYMAI) = NFROIJ(ISYMAI) + NRHFS(ISYMA)*NRHFS(ISYMI) 61 NFROFR(ISYMAI) = NFROFR(ISYMAI) 62 * + NRHFFR(ISYMA)*NRHFFR(ISYMI) 63celena 64 NV1FRO(ISYMAI) = NV1FRO(ISYMAI) + NRHFFR(ISYMI)* 65 & (NORB1(ISYMA)-NRHFFR(ISYMA)) 66celena 67 68 110 CONTINUE 69 100 CONTINUE 70C 71 DO 120 ISYMAI = 1,NSYM 72 NT2FRO(ISYMAI) = 0 73 NF2FRO(ISYMAI) = 0 74 NALLAI(ISYMAI) = 0 75 NOFROO(ISYMAI) = 0 76 NF2IJG(ISYMAI) = 0 77 NFRIJK(ISYMAI) = 0 78 NCOFRF(ISYMAI) = 0 79celena 80 NFROVF(ISYMAI) = 0 81 NFROVR(ISYMAI) = 0 82celena 83 DO 130 ISYMI = 1,NSYM 84 ISYMA = MULD2H(ISYMAI,ISYMI) 85 NT2FRO(ISYMAI) = NT2FRO(ISYMAI) + NT1FRO(ISYMA)*NT1AM(ISYMI) 86 NF2FRO(ISYMAI) = NF2FRO(ISYMAI) 87 * + NT1FRO(ISYMA)*NF1FRO(ISYMI) 88 NALLAI(ISYMAI) = NALLAI(ISYMAI) + NVIRS(ISYMA)*NRHFS(ISYMI) 89 NOFROO(ISYMAI) = NOFROO(ISYMAI) + NCOFRO(ISYMA)*NRHF(ISYMI) 90 NF2IJG(ISYMAI) = NF2IJG(ISYMAI) + NFROIJ(ISYMA)*NBAS(ISYMI) 91 NFRIJK(ISYMAI) = NFRIJK(ISYMAI) + NFROIJ(ISYMA)*NRHFS(ISYMI) 92 NCOFRF(ISYMAI) = NCOFRF(ISYMAI) 93 * + NCOFRO(ISYMA)*NRHFFR(ISYMI) 94celena 95 NFROVF(ISYMAI) = NFROVF(ISYMAI)+NFROFR(ISYMI)*NT1FRO(ISYMA) 96 NFROVR(ISYMAI) = NFROVR(ISYMAI)+NV1FRO(ISYMI)*NT1FRO(ISYMA) 97celena 98 130 CONTINUE 99 120 CONTINUE 100C 101C-------------------------------------------- 102C Calculate symmetry offset index arrays. 103C-------------------------------------------- 104C 105 DO 140 ISYMAI = 1,NSYM 106 ICOUN1 = 0 107 ICOUN2 = 0 108 ICOUN3 = 0 109 ICOUN4 = 0 110 ICOUN5 = 0 111 ICOUN6 = 0 112 ICOUN7 = 0 113 ICOUN8 = 0 114 ICOUN9 = 0 115 ICOUN10 = 0 116 ICOUN11 = 0 117 ICOUN12 = 0 118 ICOUN13 = 0 119!sonia 120 ICOUN14 = 0 121!wim 122 ICOUN15 = 0 123 ICOUN16 = 0 124 125 DO 150 ISYMI = 1,NSYM 126 ISYMA = MULD2H(ISYMAI,ISYMI) 127 IT1FRO(ISYMA,ISYMI) = ICOUN1 128 IT2FRO(ISYMA,ISYMI) = ICOUN2 129 IDSFRO(ISYMA,ISYMI) = ICOUN3 130 ICOFRO(ISYMA,ISYMI) = ICOUN4 131 IALLAI(ISYMA,ISYMI) = ICOUN5 132 IOFROO(ISYMA,ISYMI) = ICOUN6 133 IF2IJG(ISYMA,ISYMI) = ICOUN7 134 IFROIJ(ISYMA,ISYMI) = ICOUN8 135 IFRIJK(ISYMA,ISYMI) = ICOUN9 136 IFROFR(ISYMA,ISYMI) = ICOUN10 137 ICOFRF(ISYMA,ISYMI) = ICOUN11 138 IOFOAO(ISYMA,ISYMI) = ICOUN12 139 IOFFAO(ISYMA,ISYMI) = ICOUN13 140!sonia 141 ICDKFR(ISYMA,ISYMI) = ICOUN14 142!wim 143 IF1FRO(ISYMA,ISYMI) = ICOUN15 144 IF2FRO(ISYMA,ISYMI) = ICOUN16 145 146 ICOUN1 = ICOUN1 + NVIR(ISYMA)*NRHFFR(ISYMI) 147 ICOUN2 = ICOUN2 + NT1FRO(ISYMA)*NT1AM(ISYMI) 148 ICOUN3 = ICOUN3 + NNBST(ISYMA)*NRHFFR(ISYMI) 149 ICOUN4 = ICOUN4 + NRHF(ISYMA)*NRHFFR(ISYMI) 150 ICOUN5 = ICOUN5 + NVIRS(ISYMA)*NRHFS(ISYMI) 151 ICOUN6 = ICOUN6 + NCOFRO(ISYMA)*NRHF(ISYMI) 152 ICOUN7 = ICOUN7 + NFROIJ(ISYMA)*NBAS(ISYMI) 153 ICOUN8 = ICOUN8 + NRHFS(ISYMA)*NRHFS(ISYMI) 154 ICOUN9 = ICOUN9 + NFROIJ(ISYMA)*NRHFS(ISYMI) 155 ICOUN10 = ICOUN10 + NRHFFR(ISYMA)*NRHFFR(ISYMI) 156 ICOUN11 = ICOUN11 + NCOFRO(ISYMA)*NRHFFR(ISYMI) 157 ICOUN12 = ICOUN12 + NOFROO(ISYMA)*NBAS(ISYMI) 158 ICOUN13 = ICOUN13 + NCOFRF(ISYMA)*NBAS(ISYMI) 159!sonia 160 ICOUN14 = ICOUN14 + NMATAB(ISYMA)*NRHF(ISYMI) 161!wim 162 ICOUN15 = ICOUN15 + NRHFFR(ISYMA)*NRHF(ISYMI) 163 ICOUN16 = ICOUN16 + NT1FRO(ISYMA)*NF1FRO(ISYMI) 164 165 150 CONTINUE 166 140 CONTINUE 167C 168 ICOUN1 = 0 169 ICOUN2 = NRHFTS 170 DO 160 ISYM = 1,NSYM 171 IRHFS(ISYM) = ICOUN1 172 IVIRS(ISYM) = ICOUN2 173 ICOUN1 = ICOUN1 + NRHFS(ISYM) 174 ICOUN2 = ICOUN2 + NVIRS(ISYM) 175 160 CONTINUE 176C 177C---------------------- 178C Print if desired. 179C---------------------- 180C 181 IF (IPRINT .GT. 6) THEN 182 CALL AROUND('Information from CC_INIFRO') 183 WRITE(LUPRI,1) 'NT1FRO :',(NT1FRO(I), I=1,NSYM) 184 WRITE(LUPRI,1) 'NT2FRO :',(NT2FRO(I), I=1,NSYM) 185 WRITE(LUPRI,1) 'NDSFRO :',(NDSFRO(I), I=1,NSYM) 186 WRITE(LUPRI,1) 'NCOFRO :',(NCOFRO(I), I=1,NSYM) 187 WRITE(LUPRI,1) 'NALLAI :',(NALLAI(I), I=1,NSYM) 188 WRITE(LUPRI,1) 'NOFROO :',(NOFROO(I), I=1,NSYM) 189 WRITE(LUPRI,1) 'NF2IJG :',(NF2IJG(I), I=1,NSYM) 190 WRITE(LUPRI,1) 'NFROIJ :',(NFROIJ(I), I=1,NSYM) 191 WRITE(LUPRI,1) 'NFRIJK :',(NFRIJK(I), I=1,NSYM) 192 WRITE(LUPRI,1) 'NFROFR :',(NFROFR(I), I=1,NSYM) 193 WRITE(LUPRI,1) 'NCOFRF :',(NCOFRF(I), I=1,NSYM) 194 WRITE(LUPRI,1) 'IRHFS :',(IRHFS(I), I=1,NSYM) 195 WRITE(LUPRI,1) 'IVIRS :',(IVIRS(I), I=1,NSYM) 196 ENDIF 197 IF (IPRINT .GT. 10) THEN 198 WRITE(LUPRI,*) 199 DO 11 I = 1,NSYM 200 WRITE(LUPRI,1) 'IT1FRO :',(IT1FRO(I,J), J=1,NSYM) 201 11 CONTINUE 202 DO 12 I = 1,NSYM 203 WRITE(LUPRI,1) 'IT2FRO :',(IT2FRO(I,J), J=1,NSYM) 204 12 CONTINUE 205 DO 13 I = 1,NSYM 206 WRITE(LUPRI,1) 'IDSFRO :',(IDSFRO(I,J), J=1,NSYM) 207 13 CONTINUE 208 DO 14 I = 1,NSYM 209 WRITE(LUPRI,1) 'ICOFRO :',(ICOFRO(I,J), J=1,NSYM) 210 14 CONTINUE 211 DO 15 I = 1,NSYM 212 WRITE(LUPRI,1) 'IALLAI :',(IALLAI(I,J), J=1,NSYM) 213 15 CONTINUE 214 DO 16 I = 1,NSYM 215 WRITE(LUPRI,1) 'IOFROO :',(IOFROO(I,J), J=1,NSYM) 216 16 CONTINUE 217 DO 17 I = 1,NSYM 218 WRITE(LUPRI,1) 'IF2IJG :',(IF2IJG(I,J), J=1,NSYM) 219 17 CONTINUE 220 DO 18 I = 1,NSYM 221 WRITE(LUPRI,1) 'IFROIJ :',(IFROIJ(I,J), J=1,NSYM) 222 18 CONTINUE 223 DO 19 I = 1,NSYM 224 WRITE(LUPRI,1) 'IFRIJK :',(IFRIJK(I,J), J=1,NSYM) 225 19 CONTINUE 226 DO 20 I = 1,NSYM 227 WRITE(LUPRI,1) 'IFROFR :',(IFROFR(I,J), J=1,NSYM) 228 20 CONTINUE 229 DO 21 I = 1,NSYM 230 WRITE(LUPRI,1) 'ICOFRF :',(ICOFRF(I,J), J=1,NSYM) 231 21 CONTINUE 232 DO 22 I = 1,NSYM 233 WRITE(LUPRI,1) 'IOFOAO :',(IOFOAO(I,J), J=1,NSYM) 234 22 CONTINUE 235 DO 23 I = 1,NSYM 236 WRITE(LUPRI,1) 'IOFFAO :',(IOFFAO(I,J), J=1,NSYM) 237 23 CONTINUE 238 WRITE(LUPRI,*) 239 END IF 240C 241 RETURN 242C 243 1 FORMAT(3X,A8,8I8) 244C 245 END 246C /* Deck cc_kabre */ 247 SUBROUTINE CC_KABRE(SOLU,ZKAM,WORK,LWORK) 248C 249C Written by Asger Halkier 22/5 - 1998. 250C 251C To reorder the solution vector from the abacus-response 252C solver to cc-standards. 253C 254#include "implicit.h" 255C 256 DIMENSION SOLU(*), ZKAM(*), WORK(LWORK) 257C 258#include "priunit.h" 259#include "maxorb.h" 260#include "ccsdinp.h" 261#include "ccorb.h" 262#include "ccsdsym.h" 263#include "ccfro.h" 264C 265 IF (FROIMP) THEN 266C 267 KOFF1 = 2*NT1AMX + NMATIJ(1) + NMATAB(1) + 1 268 CALL DZERO(ZKAM(KOFF1),2*NT1FRO(1)) 269C 270 DO 100 ISYMI = 1,NSYM 271 ISYMA = ISYMI 272C 273 IF (NRHFFR(ISYMI) .GT. 0) THEN 274 LEN1 = NVIR(ISYMA)*NRHFFR(ISYMI) 275 KOFF2 = IALLAI(ISYMA,ISYMI) + 1 276 KOFF3 = KOFF1 + IT1FRO(ISYMA,ISYMI) 277 KOFF4 = KOFF3 + NT1FRO(1) 278 CALL DCOPY(LEN1,SOLU(KOFF2),1,ZKAM(KOFF3),1) 279 CALL DCOPY(LEN1,SOLU(KOFF2),1,ZKAM(KOFF4),1) 280 ENDIF 281C 282 LEN2 = NVIR(ISYMA)*NRHF(ISYMI) 283 KOFF5 = IALLAI(ISYMA,ISYMI) + NVIR(ISYMA)*NRHFFR(ISYMI) + 1 284 KOFF6 = NMATAB(1) + NMATIJ(1) + IT1AM(ISYMA,ISYMI) + 1 285 KOFF7 = KOFF6 + NT1AMX 286 CALL DCOPY(LEN2,SOLU(KOFF5),1,ZKAM(KOFF6),1) 287 CALL DCOPY(LEN2,SOLU(KOFF5),1,ZKAM(KOFF7),1) 288C 289 IF (IPRINT .GT. 20) THEN 290C 291 WRITE(LUPRI,*) ' ' 292 WRITE(LUPRI,1) 'Symmetry block number:', ISYMI 293 1 FORMAT(3X,A22,2X,I1) 294C 295 KOFFTO = IALLAI(ISYMA,ISYMI) + 1 296 KOFFCP = KOFF6 297 KOFFFP = KOFF1 + IT1FRO(ISYMA,ISYMI) 298C 299 CALL AROUND('Total solution vector') 300 CALL OUTPUT(SOLU(KOFFTO),1,NVIR(ISYMA),1,NRHFS(ISYMI), 301 * NVIR(ISYMA),NRHFS(ISYMI),1,LUPRI) 302 CALL AROUND('Correlated part') 303 CALL OUTPUT(ZKAM(KOFFCP),1,NVIR(ISYMA),1,NRHF(ISYMI), 304 * NVIR(ISYMA),NRHF(ISYMI),1,LUPRI) 305 CALL AROUND('Frozen part') 306 CALL OUTPUT(ZKAM(KOFFFP),1,NVIR(ISYMA),1,NRHFFR(ISYMI), 307 * NVIR(ISYMA),NRHFFR(ISYMI),1,LUPRI) 308C 309 ENDIF 310C 311 100 CONTINUE 312 ELSE 313C 314 KOFF8 = NMATIJ(1) + NMATAB(1) + 1 315 KOFF9 = KOFF8 + NT1AMX 316 CALL DCOPY(NT1AMX,SOLU(1),1,ZKAM(KOFF8),1) 317 CALL DCOPY(NT1AMX,SOLU(1),1,ZKAM(KOFF9),1) 318C 319 ENDIF 320C 321 RETURN 322 END 323C /* Deck cc_etare */ 324 SUBROUTINE CC_ETARE(ETAAI,AFROI,WORK,LWORK) 325C 326C Written by Asger Halkier 22/5 - 1998. 327C 328C To reorder the eta-kappa-0 vector for the abacus-response 329C solver. 330C 331#include "implicit.h" 332C 333 DIMENSION ETAAI(*), AFROI(*), WORK(LWORK) 334C 335#include "priunit.h" 336#include "maxorb.h" 337#include "ccsdinp.h" 338#include "ccorb.h" 339#include "ccsdsym.h" 340#include "ccfro.h" 341C 342 IF (FROIMP) THEN 343C 344 IF (LWORK .LT. NT1AMX) THEN 345 WRITE(LUPRI,*) 'Needed:', NT1AMX, 'Available:', LWORK 346 CALL QUIT('Insufficient memory in cc_etare') 347 ENDIF 348C 349 CALL DCOPY(NT1AMX,ETAAI,1,WORK,1) 350 CALL DZERO(ETAAI,NT1AMX) 351C 352 DO 100 ISYMI = 1,NSYM 353 ISYMA = ISYMI 354C 355 IF (NRHFFR(ISYMI) .GT. 0) THEN 356 LEN1 = NVIR(ISYMA)*NRHFFR(ISYMI) 357 KOFF1 = IT1FRO(ISYMA,ISYMI) + 1 358 KOFF2 = IALLAI(ISYMA,ISYMI) + 1 359 CALL DCOPY(LEN1,AFROI(KOFF1),1,ETAAI(KOFF2),1) 360 ENDIF 361C 362 LEN2 = NVIR(ISYMA)*NRHF(ISYMI) 363 KOFF3 = IT1AM(ISYMA,ISYMI) + 1 364 KOFF4 = IALLAI(ISYMA,ISYMI) + NVIR(ISYMA)*NRHFFR(ISYMI) + 1 365 CALL DCOPY(LEN2,WORK(KOFF3),1,ETAAI(KOFF4),1) 366C 367 IF (IPRINT .GT. 20) THEN 368C 369 WRITE(LUPRI,*) ' ' 370 WRITE(LUPRI,1) 'Symmetry block number:', ISYMI 371 1 FORMAT(3X,A22,2X,I1) 372C 373 KOFFTO = IALLAI(ISYMA,ISYMI) + 1 374 KOFFCP = KOFF3 375 KOFFFP = IT1FRO(ISYMA,ISYMI) + 1 376C 377 CALL AROUND('Total RHS vector') 378 CALL OUTPUT(ETAAI(KOFFTO),1,NVIR(ISYMA),1,NRHFS(ISYMI), 379 * NVIR(ISYMA),NRHFS(ISYMI),1,LUPRI) 380 CALL AROUND('Correlated part') 381 CALL OUTPUT(WORK(KOFFCP),1,NVIR(ISYMA),1,NRHF(ISYMI), 382 * NVIR(ISYMA),NRHF(ISYMI),1,LUPRI) 383 CALL AROUND('Frozen part') 384 CALL OUTPUT(AFROI(KOFFFP),1,NVIR(ISYMA),1,NRHFFR(ISYMI), 385 * NVIR(ISYMA),NRHFFR(ISYMI),1,LUPRI) 386C 387 ENDIF 388C 389 100 CONTINUE 390 ENDIF 391C 392 RETURN 393 END 394C /* Deck mp2_zkfcb */ 395 SUBROUTINE MP2_ZKFCB(IPDD,R12PRP,RES,T2MAM,WORK,LWORK) 396C 397C Written by Asger Halkier 25/5 - 1998. 398C 399C To calculate the frozen-occupied block (iJ) of kappa-bar-0. 400C 401#include "implicit.h" 402#include "dummy.h" 403C 404 PARAMETER (ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0) 405 DIMENSION RES(*), T2MAM(*), WORK(LWORK) 406 INTEGER YIJIJ, YIJIJT 407 LOGICAL R12PRP 408C 409#include "priunit.h" 410#include "maxorb.h" 411#include "ccsdinp.h" 412#include "ccorb.h" 413#include "ccsdsym.h" 414#include "ccfro.h" 415C 416C--------------------------- 417C Work space allocation. 418C--------------------------- 419C 420 KINTFR = 1 421 KPIMIJ = KINTFR + NT2FRO(1) 422 KFOCKD = KPIMIJ + NCOFRO(1) 423 KEND1 = KFOCKD + NORBTS 424 LWRK1 = LWORK - KEND1 425 426celena 427C------------------------------------- 428C Read R12-contributions 429C------------------------------------- 430 431 IF (R12PRP .AND. (IPDD .EQ. 2 .OR. IPDD .EQ. 3 432 & .OR. IPDD .EQ. 5)) THEN 433 DO ISYMAJ = 1,NSYM 434 ISYMIJ = ISYMAJ 435 NCVIJ = 0 436 DO ISYMA = 1,NSYM 437 ISYMJ = MULD2H(ISYMAJ,ISYMA) 438 ISYMI = MULD2H(ISYMIJ,ISYMJ) 439 NCVIJ = NCVIJ + NRHFFR(ISYMA)*NRHF(ISYMI) 440 ENDDO 441 ENDDO 442 YIJIJ = KEND1 443 YIJIJT = YIJIJ + NCVIJ 444 KEND1 = YIJIJT+ NCVIJ 445 LWRK1 = LWORK - KEND1 446 LUVAJKL = -1 447 IF (IPDD .EQ. 2) THEN 448 CALL GPOPEN(LUVAJKL,'CCR12YIJIJ','UNKNOWN',' ', 449 * 'UNFORMATTED',IDUMMY,.FALSE.) 450 ELSEIF (IPDD .EQ. 3) THEN 451 CALL GPOPEN(LUVAJKL,'CCR12ZIJIJ','UNKNOWN',' ', 452 * 'UNFORMATTED',IDUMMY,.FALSE.) 453 ELSEIF (IPDD .EQ. 5) THEN 454 CALL GPOPEN(LUVAJKL,'CCR12XIJIJ','UNKNOWN',' ', 455 * 'UNFORMATTED',IDUMMY,.FALSE.) 456 ENDIF 457 REWIND(LUVAJKL) 458 READ(LUVAJKL) (WORK(YIJIJ+I-1),I=1,NCVIJ) 459 CALL GPCLOSE(LUVAJKL,'DELETE') 460 DO ISYMAJ = 1,NSYM 461 ISYMIJ = ISYMAJ 462 DO ISYMA = 1,NSYM 463 ISYMJ = MULD2H(ISYMA,ISYMAJ) 464 ISYMI = MULD2H(ISYMA,ISYMIJ) 465 do i =1,nrhf(isymi) 466 do j =1,nrhffr(isymj) 467 idxij = icofro(isymj,isymi)+ 468 & nrhffr(isymj)*(i-1) + j 469 idxji = icofro(isymi,isymj)+ 470 & nrhf(isymi)*(j-1) + i 471 work(YIJIJT+idxji-1) = work(YIJIJ+idxij-1) 472 enddo 473 enddo 474 ENDDO 475 ENDDO 476 ENDIF 477celena 478 479 480C 481 IF (LWRK1 .LT. 0) THEN 482 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 483 CALL QUIT('Insufficient work space for allocation '// 484 & 'in MP2_ZKFCB') 485 ENDIF 486C 487 CALL DZERO(WORK(KPIMIJ),NCOFRO(1)) 488 CALL DZERO(WORK(KFOCKD),NORBTS) 489C 490C-------------------------------------- 491C Read integrals (cJ|dk) from disk. 492C-------------------------------------- 493C 494 LUCJDK = -1 495 CALL GPOPEN(LUCJDK,'INCJDK','UNKNOWN',' ','UNFORMATTED',IDUMMY, 496 & .FALSE.) 497 REWIND(LUCJDK) 498 READ(LUCJDK) (WORK(KINTFR+I-1), I = 1,NT2FRO(1)) 499 CALL GPCLOSE(LUCJDK,'KEEP') 500 501C 502 IF (IPRINT. GT. 20) THEN 503 XFRNOR = DDOT(NT2FRO(1),WORK(KINTFR),1,WORK(KINTFR),1) 504 WRITE(LUPRI,*) 'Norm of integrals (cJdk):', XFRNOR 505 ENDIF 506C 507C--------------------------------------- 508C Contract integrals with amplitude. 509C--------------------------------------- 510C 511 DO 100 ISYMDK = 1,NSYM 512 ISYMCI = ISYMDK 513 ISYMCJ = ISYMDK 514 DO 110 NDK = 1,NT1AM(ISYMDK) 515 DO 120 ISYMC = 1,NSYM 516 ISYMI = MULD2H(ISYMC,ISYMCI) 517 ISYMJ = MULD2H(ISYMC,ISYMCJ) 518C 519 KOFF1 = IT2SQ(ISYMCI,ISYMDK) 520 * + NT1AM(ISYMCI)*(NDK - 1) 521 * + IT1AM(ISYMC,ISYMI) + 1 522 KOFF2 = KINTFR + IT2FRO(ISYMCJ,ISYMDK) 523 * + NT1FRO(ISYMCJ)*(NDK - 1) 524 * + IT1FRO(ISYMC,ISYMJ) 525 KOFF3 = KPIMIJ + ICOFRO(ISYMI,ISYMJ) 526C 527 NTOTC = MAX(NVIR(ISYMC),1) 528 NTOTI = MAX(NRHF(ISYMI),1) 529C 530 CALL DGEMM('T','N',NRHF(ISYMI),NRHFFR(ISYMJ), 531 * NVIR(ISYMC),ONE,T2MAM(KOFF1),NTOTC, 532 * WORK(KOFF2),NTOTC,ONE,WORK(KOFF3),NTOTI) 533cC 534 120 CONTINUE 535 110 CONTINUE 536 100 CONTINUE 537celena 538 IF (R12PRP) THEN 539 CALL DAXPY(NCVIJ,ONE,WORK(YIJIJT),1,WORK(KPIMIJ),1) 540 END IF 541celena 542C 543C-------------------------- 544C Get orbital energies. 545C-------------------------- 546C 547 CALL FOCK_ALL(WORK(KFOCKD),WORK(KEND1),LWRK1) 548C 549C--------------------------- 550C Calculate the results. 551C--------------------------- 552C 553 DO 130 ISYMI = 1,NSYM 554 ISYMJ = ISYMI 555 DO 140 J = 1,NRHFFR(ISYMJ) 556 DO 150 I = 1,NRHF(ISYMI) 557 NIJ = ICOFRO(ISYMI,ISYMJ) + NRHF(ISYMI)*(J - 1) + I 558 KOI = KFOCKD + IRHFS(ISYMI) + NRHFFR(ISYMI) + I - 1 559 KOJ = KFOCKD + IRHFS(ISYMJ) + J - 1 560C 561 DENOM = WORK(KOI) - WORK(KOJ) 562 RES(NIJ) = HALF*WORK(KPIMIJ+NIJ-1)/DENOM 563C 564 150 CONTINUE 565 140 CONTINUE 566 130 CONTINUE 567C 568 CALL DCOPY(NCOFRO(1),RES(1),1,RES(1+NCOFRO(1)),1) 569C 570 RETURN 571 END 572C /* Deck cmo_all */ 573 SUBROUTINE CMO_ALL(CMO,WORK,LWORK) 574C 575C Written by Asger Halkier 25/5 - 1998 576C 577C To read in and change the symmetry ordering of the FULL 578C MO coefficient matrix from Sirius to CC-ordering. 579C 580#include "implicit.h" 581#include "priunit.h" 582#include "dummy.h" 583 DIMENSION CMO(*),WORK(LWORK) 584#include "inftap.h" 585#include "ccorb.h" 586#include "ccsdinp.h" 587#include "ccsdsym.h" 588#include "r12int.h" 589C 590C---------------------------------------------- 591C Read MO-coefficients from interface file. 592C---------------------------------------------- 593C 594 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY, 595 & .FALSE.) 596 REWIND LUSIFC 597C 598C Use LABEL instead of 'TRCCINT ' (WK/UniKA/04-11-2002). 599 CALL MOLLAB(LABEL,LUSIFC,LUPRI) 600 READ (LUSIFC) 601C 602 READ (LUSIFC) 603 READ (LUSIFC) (CMO(I), I=1,NLAMDS) 604C 605 CALL GPCLOSE(LUSIFC,'KEEP') 606C 607C--------------------------- 608C Work space allocation. 609C--------------------------- 610C 611 KSCR1 = 1 612 KEND1 = KSCR1 + NLAMDS 613 LWRK1 = LWORK - KEND1 614C 615 IF (LWRK1 .LT. 0) THEN 616 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 617 CALL QUIT('Insufficient work space for allocation in CMO_ALL') 618 ENDIF 619C 620C---------------------------------- 621C Reorder all orbitals in work. 622C---------------------------------- 623C 624 ICRHF = KSCR1 625 ICVIR = KSCR1 + NLRHSI 626C 627 ICOUNT = 1 628 DO 100 ISYM = 1,NSYM 629C 630 CALL DCOPY(NBAS(ISYM)*NRHFS(ISYM),CMO(ICOUNT),1,WORK(ICRHF),1) 631 ICRHF = ICRHF + NBAS(ISYM)*NRHFS(ISYM) 632 ICOUNT = ICOUNT + NBAS(ISYM)*NRHFS(ISYM) 633C 634 CALL DCOPY(NBAS(ISYM)*NVIRS(ISYM),CMO(ICOUNT),1,WORK(ICVIR),1) 635 ICVIR = ICVIR + NBAS(ISYM)*NVIRS(ISYM) 636 ICOUNT = ICOUNT + NBAS(ISYM)*NVIRS(ISYM) 637C 638 100 CONTINUE 639C 640C----------------------------------------- 641C Copy reordered result back into CMO. 642C----------------------------------------- 643C 644 CALL DCOPY(NLAMDS,WORK(KSCR1),1,CMO,1) 645C 646 RETURN 647 END 648C /* Deck cc_frcoin */ 649 SUBROUTINE CC_FRCOIN(XFRIN,XINT,CMO,WORK,LWORK,IDEL,ISYDEL) 650C 651C Written by Asger Halkier 25/5 - 1998. 652C 653C To calculate the integrals (cJ|dk) needed for frozen-core gradients. 654C AO integrals are assumed totally symmetric. 655C 656#include "implicit.h" 657C 658 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0) 659 DIMENSION XFRIN(*), XINT(*), CMO(*), WORK(LWORK) 660C 661#include "priunit.h" 662#include "maxorb.h" 663#include "ccsdinp.h" 664#include "ccorb.h" 665#include "ccsdsym.h" 666#include "ccfro.h" 667#include "r12int.h" 668C 669C--------------------------------------- 670C Initial symmetries and allocation. 671C--------------------------------------- 672C 673 ISYDIS = ISYDEL 674 ISYMD = ISYDEL 675C 676 KDVEC = 1 677 KEND1 = KDVEC + NVIR(ISYMD) 678 LWRK1 = LWORK - KEND1 679C 680 IF (LWRK1 .LT. 0) THEN 681 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 682 CALL QUIT('Insufficient work space for allocation '// 683 & 'in CC_FRCOIN') 684 ENDIF 685C 686 CALL DZERO(WORK(KDVEC),NVIR(ISYMD)) 687C 688C----------------------------------- 689C Copy vector out of CMO matrix. 690C----------------------------------- 691C 692 KOFF1 = ILVISI(ISYMD) + IDEL - IBAS(ISYDEL) 693C 694 CALL DCOPY(NVIR(ISYMD),CMO(KOFF1),NBAS(ISYDEL),WORK(KDVEC),1) 695C 696C-------------------------------------------------------- 697C Outer symmetry loop and next work space allocation. 698C-------------------------------------------------------- 699C 700 DO 100 ISYMK = 1,NSYM 701C 702 ISYMG = ISYMK 703 ISYMDK = MULD2H(ISYMD,ISYMK) 704 ISYMCJ = ISYMDK 705 ISALBE = ISYMCJ 706C 707 KINAOK = KEND1 708 KSCR1 = KINAOK + NNBST(ISALBE)*NRHF(ISYMK) 709 KSCR2 = KSCR1 + N2BST(ISALBE) 710 KEND2 = KSCR2 + NT1FRO(ISYMCJ) 711 LWRK2 = LWORK - KEND2 712C 713 IF (LWRK2 .LT. 0) THEN 714 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND2 715 CALL QUIT('Insufficient space for allocation in CC_FRCOIN') 716 ENDIF 717C 718 CALL DZERO(WORK(KINAOK),NNBST(ISALBE)*NRHF(ISYMK)) 719C 720C-------------------------------------------- 721C Transform gamma-index to correlated. 722C-------------------------------------------- 723C 724 KOFF2 = IDSAOG(ISYMG,ISYDIS) + 1 725 KOFF3 = ILRHSI(ISYMG) + NBAS(ISYMG)*NRHFFR(ISYMK) + 1 726C 727 NTOTAB = MAX(NNBST(ISALBE),1) 728 NTOTG = MAX(NBAS(ISYMG),1) 729C 730 CALL DGEMM('N','N',NNBST(ISALBE),NRHF(ISYMK),NBAS(ISYMG),ONE, 731 * XINT(KOFF2),NTOTAB,CMO(KOFF3),NTOTG,ZERO, 732 * WORK(KINAOK),NTOTAB) 733 734 735 DO 110 K = 1,NRHF(ISYMK) 736C 737 CALL DZERO(WORK(KSCR1),N2BST(ISALBE)) 738 CALL DZERO(WORK(KSCR2),NT1FRO(ISYMCJ)) 739C 740C----------------------------------- 741C Square up the integrals. 742C----------------------------------- 743C 744 KOFF4 = KINAOK + NNBST(ISALBE)*(K - 1) 745 CALL CCSD_SYMSQ(WORK(KOFF4),ISALBE,WORK(KSCR1)) 746C 747C--------------------------------------------------------------- 748C Inner symmetry loop and final work space allocation. 749C--------------------------------------------------------------- 750C 751 DO 120 ISYMAL = 1,NSYM 752C 753 ISYMC = ISYMAL 754 ISYMBE = MULD2H(ISALBE,ISYMAL) 755 ISYMJ = ISYMBE 756C 757 KSCR3 = KEND2 758 KEND3 = KSCR3 + NBAS(ISYMAL)*NRHFFR(ISYMJ) 759 LWRK3 = LWORK - KEND3 760C 761 IF (LWRK3 .LT. 0) THEN 762 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND3 763 CALL QUIT('Insufficient space for allocation '// 764 & 'in CC_FRCOIN') 765 ENDIF 766C 767C------------------------------------------------ 768C Construct the integrals (cJ|kdel). 769C------------------------------------------------ 770C 771 KOFF5 = KSCR1 + IAODIS(ISYMAL,ISYMBE) 772 KOFF6 = ILRHSI(ISYMBE) + 1 773C 774 NTOTAL = MAX(NBAS(ISYMAL),1) 775 NTOTBE = MAX(NBAS(ISYMBE),1) 776C 777 CALL DGEMM('N','N',NBAS(ISYMAL),NRHFFR(ISYMJ), 778 * NBAS(ISYMBE),ONE,WORK(KOFF5),NTOTAL, 779 * CMO(KOFF6),NTOTBE,ZERO,WORK(KSCR3),NTOTAL) 780C 781 KOFF7 = ILVISI(ISYMAL) + 1 782 KOFF8 = KSCR2 + IT1FRO(ISYMC,ISYMJ) 783C 784 NTOTAL = MAX(NBAS(ISYMAL),1) 785 NTOTC = MAX(NVIR(ISYMC),1) 786C 787 CALL DGEMM('T','N',NVIR(ISYMC),NRHFFR(ISYMJ), 788 * NBAS(ISYMAL),ONE,CMO(KOFF7),NTOTAL, 789 * WORK(KSCR3),NTOTAL,ZERO,WORK(KOFF8),NTOTC) 790C 791 120 CONTINUE 792C 793C---------------------------------------------------------------- 794C Final scaling with CMO element and storage in result. 795C---------------------------------------------------------------- 796C 797 IF (R12TRA .AND. .NOT. R12PRP) THEN 798 DO 131 D = 1,NRHFFR(ISYMD) 799C 800 NDK = IF1FRO(ISYMD,ISYMK) + NRHFFR(ISYMD)*(K - 1) + D 801 KOFF9 = KDVEC + D - 1 802 KOFF10 = IF2FRO(ISYMCJ,ISYMDK) 803 * + NT1FRO(ISYMCJ)*(NDK - 1) + 1 804C 805 CALL DAXPY(NT1FRO(ISYMCJ),WORK(KOFF9),WORK(KSCR2),1, 806 * XFRIN(KOFF10),1) 807C 808 131 CONTINUE 809C 810 ELSE 811 DO 130 D = 1,NVIR(ISYMD) 812C 813 NDK = IT1AM(ISYMD,ISYMK) + NVIR(ISYMD)*(K - 1) + D 814 KOFF9 = KDVEC + D - 1 815 KOFF10 = IT2FRO(ISYMCJ,ISYMDK) 816 * + NT1FRO(ISYMCJ)*(NDK - 1) + 1 817C 818 CALL DAXPY(NT1FRO(ISYMCJ),WORK(KOFF9),WORK(KSCR2),1, 819 * XFRIN(KOFF10),1) 820 130 CONTINUE 821 END IF 822 110 CONTINUE 823 100 CONTINUE 824C 825 RETURN 826 END 827C /* Deck cc_frcogr */ 828 SUBROUTINE CC_FRCOGR(XFRIN,XINT,CMO,WORK,LWORK,IDEL,ISYDEL) 829C 830C 831C To calculate the integrals (pJ|Ik) needed for frozen-core gradients. 832C AO integrals are assumed totally symmetric. 833c based on cc_frcoin (Elena Vollmer 29/09/04) 834C 835#include "implicit.h" 836C 837 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0) 838 DIMENSION XFRIN(*), XINT(*), CMO(*), WORK(LWORK) 839 integer ilorbi(8),iv2fro(8,8),iv1fro(8,8) 840C 841#include "priunit.h" 842#include "maxorb.h" 843#include "ccsdinp.h" 844#include "ccorb.h" 845#include "ccsdsym.h" 846#include "ccfro.h" 847#include "r12int.h" 848C 849C--------------------------------------- 850C Initial symmetries and allocation. 851C--------------------------------------- 852C 853celena 854 ISYDIS = ISYDEL 855 ISYMD = ISYDEL 856 857 858 noffset = 0 859 icoun1 = 0 860 do isym = 1,nsym 861 ilorbi(isym) = icoun1 862 icoun1 = icoun1 + nbas(isym)*(nvir(isym)-nrhffr(isym)) 863 noffset = noffset + nbas(isym)*(nrhffr(isym)) 864 enddo 865 866 do isymk = 1, nsym 867 icoun3 = 0 868 icoun4 = 0 869 do isymj = 1,nsym 870 isymi = muld2h(isymk,isymj) 871 iv2fro(isymi,isymj) = icoun3 872 iv1fro(isymi,isymj) = icoun4 873 icoun3 = icoun3 + nt1fro(isymi)*nv1fro(isymj) 874 icoun4 = icoun4 + (nrhffr(isymi)) 875 & *(norb1(isymj)-nrhffr(isymj)) 876 enddo 877 enddo 878 879celena 880 881C 882 KDVEC = 1 883 KEND1 = KDVEC + NVIR(ISYMD) 884 LWRK1 = LWORK - KEND1 885C 886 IF (LWRK1 .LT. 0) THEN 887 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 888 CALL QUIT('Insufficient work space for allocation '// 889 & 'in CC_FRCOIN') 890 ENDIF 891C 892 CALL DZERO(WORK(KDVEC),NVIR(ISYMD)) 893C 894C----------------------------------- 895C Copy vector out of CMO matrix. 896C----------------------------------- 897C 898 KOFF1 = ILVISI(ISYMD) + IDEL - IBAS(ISYDEL) 899C 900 CALL DCOPY(NVIR(ISYMD),CMO(KOFF1),NBAS(ISYDEL),WORK(KDVEC),1) 901C 902C-------------------------------------------------------- 903C Outer symmetry loop and next work space allocation. 904C-------------------------------------------------------- 905C 906 DO 100 ISYMK = 1,NSYM 907C 908 ISYMG = ISYMK 909 ISYMDK = MULD2H(ISYMD,ISYMK) 910 ISYMCJ = ISYMDK 911 ISALBE = ISYMCJ 912C 913 KINAOK = KEND1 914 KSCR1 = KINAOK + NNBST(ISALBE)*(NORB1(ISYMK)-NRHFFR(ISYMK)) 915 KSCR2 = KSCR1 + N2BST(ISALBE) 916 KEND2 = KSCR2 + NT1FRO(ISYMCJ) 917 LWRK2 = LWORK - KEND2 918C 919 IF (LWRK2 .LT. 0) THEN 920 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND2 921 CALL QUIT('Insufficient space for allocation in CC_FRCOGR') 922 ENDIF 923C 924 CALL DZERO(WORK(KINAOK),NNBST(ISALBE)*(NORB1(ISYMK) 925 & -NRHFFR(ISYMK))) 926C 927C-------------------------------------------- 928C Transform gamma-index to correlated. 929C-------------------------------------------- 930C 931 KOFF2 = IDSAOG(ISYMG,ISYDIS) + 1 932 KOFF3 = NBAS(ISYMG)*NRHFFR(ISYMK)+ILVISI(ISYMG) + 1 933C 934 NTOTAB = MAX(NNBST(ISALBE),1) 935 NTOTG = MAX(NBAS(ISYMG),1) 936C 937 CALL DGEMM('N','N',NNBST(ISALBE),(NORB1(ISYMK)-NRHFFR(ISYMK)), 938 & NBAS(ISYMG),ONE, 939 * XINT(KOFF2),NTOTAB,CMO(KOFF3),NTOTG,ZERO, 940 * WORK(KINAOK),NTOTAB) 941 942 943 944 DO 110 K = 1,(NORB1(ISYMK)-NRHFFR(ISYMK)) 945C 946 CALL DZERO(WORK(KSCR1),N2BST(ISALBE)) 947 CALL DZERO(WORK(KSCR2),NT1FRO(ISYMCJ)) 948C 949C----------------------------------- 950C Square up the integrals. 951C----------------------------------- 952C 953 KOFF4 = KINAOK + NNBST(ISALBE)*(K - 1) 954 CALL CCSD_SYMSQ(WORK(KOFF4),ISALBE,WORK(KSCR1)) 955C 956C--------------------------------------------------------------- 957C Inner symmetry loop and final work space allocation. 958C--------------------------------------------------------------- 959C 960 DO 120 ISYMAL = 1,NSYM 961C 962 ISYMC = ISYMAL 963 ISYMBE = MULD2H(ISALBE,ISYMAL) 964 ISYMJ = ISYMBE 965C 966 KSCR3 = KEND2 967 KEND3 = KSCR3 + NBAS(ISYMAL)*NRHFFR(ISYMJ) 968 LWRK3 = LWORK - KEND3 969C 970 IF (LWRK3 .LT. 0) THEN 971 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND3 972 CALL QUIT('Insufficient space for allocation '// 973 & 'in CC_FRCOIN') 974 ENDIF 975C 976C------------------------------------------------ 977C Construct the integrals (cJ|kdel). 978C------------------------------------------------ 979C 980 KOFF5 = KSCR1 + IAODIS(ISYMAL,ISYMBE) 981 KOFF6 = ILRHSI(ISYMBE) + 1 982C 983 NTOTAL = MAX(NBAS(ISYMAL),1) 984 NTOTBE = MAX(NBAS(ISYMBE),1) 985C 986 CALL DGEMM('N','N',NBAS(ISYMAL),NRHFFR(ISYMJ), 987 * NBAS(ISYMBE),ONE,WORK(KOFF5),NTOTAL, 988 * CMO(KOFF6),NTOTBE,ZERO,WORK(KSCR3),NTOTAL) 989C 990 991 KOFF7 = ILVISI(ISYMAL) + 1 992 KOFF8 = KSCR2 + IT1FRO(ISYMC,ISYMJ) 993C 994 NTOTAL = MAX(NBAS(ISYMAL),1) 995 NTOTC = MAX(NVIR(ISYMC),1) 996C 997 CALL DGEMM('T','N',NVIR(ISYMC),NRHFFR(ISYMJ), 998 * NBAS(ISYMAL),ONE,CMO(KOFF7),NTOTAL, 999 * WORK(KSCR3),NTOTAL,ZERO,WORK(KOFF8),NTOTC) 1000C 1001 120 CONTINUE 1002C 1003C---------------------------------------------------------------- 1004C Final scaling with CMO element and storage in result. 1005C---------------------------------------------------------------- 1006C 1007 DO 131 D = 1,NRHFFR(ISYMD) 1008C 1009 NDK = IV1FRO(ISYMD,ISYMK) + NRHFFR(ISYMD)*(K - 1) + D 1010 KOFF9 = KDVEC + D - 1 1011 KOFF10 = IV2FRO(ISYMCJ,ISYMDK) 1012 * + NT1FRO(ISYMCJ)*(NDK - 1) + 1 1013C 1014 CALL DAXPY(NT1FRO(ISYMCJ),WORK(KOFF9),WORK(KSCR2),1, 1015 * XFRIN(KOFF10),1) 1016C 1017 131 CONTINUE 1018C 1019 110 CONTINUE 1020 100 CONTINUE 1021C 1022 RETURN 1023 END 1024C /* Deck cc_frcogr1 */ 1025 SUBROUTINE CC_FRCOGR1(XFRIN,XINT,CMO,WORK,LWORK,IDEL,ISYDEL) 1026C 1027C 1028C To calculate the integrals (pJ|IK) needed for frozen-core gradients. 1029C AO integrals are assumed totally symmetric. 1030c based on cc_frcoin (Elena Vollmer 29/09/04) 1031C 1032#include "implicit.h" 1033C 1034 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0) 1035 DIMENSION XFRIN(*), XINT(*), CMO(*), WORK(LWORK) 1036 integer iglmrf(8,8),ilorbi(8),iv2fro(8,8),iv1fro(8,8) 1037C 1038#include "priunit.h" 1039#include "maxorb.h" 1040#include "ccsdinp.h" 1041#include "ccorb.h" 1042#include "ccsdsym.h" 1043#include "ccfro.h" 1044#include "r12int.h" 1045C 1046C--------------------------------------- 1047C Initial symmetries and allocation. 1048C--------------------------------------- 1049C 1050celena 1051 ISYDIS = ISYDEL 1052 ISYMD = ISYDEL 1053 1054 do isymmua = 1,nsym 1055 icou2 = 0 1056 do isyma = 1,nsym 1057 isymmu = muld2h(isymmua,isyma) 1058 iglmrf(isymmu,isyma) = icou2 1059 icou2 = icou2 + nbas(isymmu)*norb1(isyma) 1060 enddo 1061 enddo 1062 1063 1064 noffset = 0 1065 icoun1 = 0 1066 do isym = 1,nsym 1067 ilorbi(isym) = icoun1 1068 icoun1 = icoun1 + nbas(isym)*(nrhfs(isym)) 1069 noffset = noffset + nbas(isym)*(nrhffr(isym)) 1070 enddo 1071 1072 do isymk = 1, nsym 1073 icoun3 = 0 1074 icoun4 = 0 1075 do isymj = 1,nsym 1076 isymi = muld2h(isymk,isymj) 1077 iv2fro(isymi,isymj) = icoun3 1078 iv1fro(isymi,isymj) = icoun4 1079 icoun3 = icoun3 + nt1fro(isymi)*nfrofr(isymj) 1080 icoun4 = icoun4 + (nrhffr(isymi)) 1081 & *(norb1(isymj)-nrhffr(isymj)) 1082 enddo 1083 enddo 1084 1085celena 1086 1087C 1088 KDVEC = 1 1089 KEND1 = KDVEC + NVIR(ISYMD) 1090 LWRK1 = LWORK - KEND1 1091C 1092 IF (LWRK1 .LT. 0) THEN 1093 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 1094 CALL QUIT('Insufficient work space for allocation '// 1095 & 'in CC_FRCOIN') 1096 ENDIF 1097C 1098 CALL DZERO(WORK(KDVEC),NVIR(ISYMD)) 1099C 1100C----------------------------------- 1101C Copy vector out of CMO matrix. 1102C----------------------------------- 1103C 1104 KOFF1 = ILVISI(ISYMD) + IDEL - IBAS(ISYDEL) 1105C 1106 CALL DCOPY(NVIR(ISYMD),CMO(KOFF1),NBAS(ISYDEL),WORK(KDVEC),1) 1107C 1108C-------------------------------------------------------- 1109C Outer symmetry loop and next work space allocation. 1110C-------------------------------------------------------- 1111C 1112 DO 100 ISYMK = 1,NSYM 1113C 1114 ISYMG = ISYMK 1115 ISYMDK = MULD2H(ISYMD,ISYMK) 1116 ISYMCJ = ISYMDK 1117 ISALBE = ISYMCJ 1118C 1119 KINAOK = KEND1 1120 KSCR1 = KINAOK + NNBST(ISALBE)*NRHFFR(ISYMK) 1121 KSCR2 = KSCR1 + N2BST(ISALBE) 1122 KEND2 = KSCR2 + NT1FRO(ISYMCJ) 1123 LWRK2 = LWORK - KEND2 1124C 1125 IF (LWRK2 .LT. 0) THEN 1126 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND2 1127 CALL QUIT('Insufficient space for allocation in CC_FRCOGR') 1128 ENDIF 1129C 1130 CALL DZERO(WORK(KINAOK),NNBST(ISALBE)*NRHFFR(ISYMK)) 1131C 1132C-------------------------------------------- 1133C Transform gamma-index to correlated. 1134C-------------------------------------------- 1135C 1136 KOFF2 = IDSAOG(ISYMG,ISYDIS) + 1 1137c KOFF3 = 1 + iglmrf(ISYMg,ISYMk) 1138 KOFF3 = 1 + ilorbi(ISYMg) 1139C 1140 NTOTAB = MAX(NNBST(ISALBE),1) 1141 NTOTG = MAX(NBAS(ISYMG),1) 1142C 1143 CALL DGEMM('N','N',NNBST(ISALBE),NRHFFR(ISYMK), 1144 & NBAS(ISYMG),ONE, 1145 * XINT(KOFF2),NTOTAB,CMO(KOFF3),NTOTG,ZERO, 1146 * WORK(KINAOK),NTOTAB) 1147 1148 1149 DO 110 K = 1,NRHFFR(ISYMK) 1150C 1151 CALL DZERO(WORK(KSCR1),N2BST(ISALBE)) 1152 CALL DZERO(WORK(KSCR2),NT1FRO(ISYMCJ)) 1153C 1154C----------------------------------- 1155C Square up the integrals. 1156C----------------------------------- 1157C 1158 KOFF4 = KINAOK + NNBST(ISALBE)*(K - 1) 1159 CALL CCSD_SYMSQ(WORK(KOFF4),ISALBE,WORK(KSCR1)) 1160C 1161C--------------------------------------------------------------- 1162C Inner symmetry loop and final work space allocation. 1163C--------------------------------------------------------------- 1164C 1165 DO 120 ISYMAL = 1,NSYM 1166C 1167 ISYMC = ISYMAL 1168 ISYMBE = MULD2H(ISALBE,ISYMAL) 1169 ISYMJ = ISYMBE 1170C 1171 KSCR3 = KEND2 1172 KEND3 = KSCR3 + NBAS(ISYMAL)*NRHFFR(ISYMJ) 1173 LWRK3 = LWORK - KEND3 1174C 1175 IF (LWRK3 .LT. 0) THEN 1176 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND3 1177 CALL QUIT('Insufficient space for allocation '// 1178 & 'in CC_FRCOIN') 1179 ENDIF 1180C 1181C------------------------------------------------ 1182C Construct the integrals (cJ|kdel). 1183C------------------------------------------------ 1184C 1185 KOFF5 = KSCR1 + IAODIS(ISYMAL,ISYMBE) 1186 KOFF6 = ILRHSI(ISYMBE) + 1 1187C 1188 NTOTAL = MAX(NBAS(ISYMAL),1) 1189 NTOTBE = MAX(NBAS(ISYMBE),1) 1190C 1191 CALL DGEMM('N','N',NBAS(ISYMAL),NRHFFR(ISYMJ), 1192 * NBAS(ISYMBE),ONE,WORK(KOFF5),NTOTAL, 1193 * CMO(KOFF6),NTOTBE,ZERO,WORK(KSCR3),NTOTAL) 1194C 1195 1196 KOFF7 = ILVISI(ISYMAL) + 1 1197 KOFF8 = KSCR2 + IT1FRO(ISYMC,ISYMJ) 1198C 1199 NTOTAL = MAX(NBAS(ISYMAL),1) 1200 NTOTC = MAX(NVIR(ISYMC),1) 1201C 1202 CALL DGEMM('T','N',NVIR(ISYMC),NRHFFR(ISYMJ), 1203 * NBAS(ISYMAL),ONE,CMO(KOFF7),NTOTAL, 1204 * WORK(KSCR3),NTOTAL,ZERO,WORK(KOFF8),NTOTC) 1205C 1206 120 CONTINUE 1207C 1208C---------------------------------------------------------------- 1209C Final scaling with CMO element and storage in result. 1210C---------------------------------------------------------------- 1211C 1212 DO 131 D = 1,NRHFFR(ISYMD) 1213C 1214 NDK = IFROFR(ISYMD,ISYMK) + NRHFFR(ISYMD)*(K - 1) + D 1215 KOFF9 = KDVEC + D - 1 1216 KOFF10 = Iv2FRO(ISYMCJ,ISYMDK) 1217 * + NT1FRO(ISYMCJ)*(NDK - 1) + 1 1218C 1219 CALL DAXPY(NT1FRO(ISYMCJ),WORK(KOFF9),WORK(KSCR2),1, 1220 * XFRIN(KOFF10),1) 1221C 1222 131 CONTINUE 1223C 1224 110 CONTINUE 1225 100 CONTINUE 1226C 1227 RETURN 1228 END 1229C /* Deck cc_frcr12 */ 1230 SUBROUTINE CC_FRCR12(XFRIN,XINT,CMO,WORK,LWORK,IDEL,ISYDEL) 1231C 1232C R12 version of CC_FRCOIN (WK/UniKA/20-08-2003) 1233C 1234C Transform (ab|cd) into (p'J|kL), with J and L frozen. 1235C 1236#include "implicit.h" 1237C 1238 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0) 1239 DIMENSION XFRIN(*), XINT(*), CMO(*), WORK(LWORK) 1240C 1241#include "priunit.h" 1242#include "maxorb.h" 1243#include "ccsdinp.h" 1244#include "ccorb.h" 1245#include "ccsdsym.h" 1246#include "ccfro.h" 1247#include "r12int.h" 1248C 1249C--------------------------------------- 1250C Initial symmetries and allocation. 1251C--------------------------------------- 1252C 1253 ISYDIS = ISYDEL 1254 ISYMD = ISYDEL 1255C 1256 KDVEC = 1 1257 KEND1 = KDVEC + NVIR(ISYMD) 1258 LWRK1 = LWORK - KEND1 1259C 1260 IF (LWRK1 .LT. 0) THEN 1261 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 1262 CALL QUIT('Insufficient work space for allocation '// 1263 & 'in CC_FRCR12') 1264 ENDIF 1265C 1266C----------------------------------- 1267C Copy vector out of CMO matrix. 1268C----------------------------------- 1269C 1270 KOFF1 = ILVISI(ISYMD) + IDEL - IBAS(ISYDEL) 1271C 1272 CALL DCOPY(NVIR(ISYMD),CMO(KOFF1),NBAS(ISYDEL),WORK(KDVEC),1) 1273C 1274C-------------------------------------------------------- 1275C Outer symmetry loop and next work space allocation. 1276C-------------------------------------------------------- 1277C 1278 DO 100 ISYMK = 1,NSYM 1279C 1280 ISYMG = ISYMK 1281 ISYMDK = MULD2H(ISYMD,ISYMK) 1282 ISYMCJ = ISYMDK 1283 ISALBE = ISYMCJ 1284C 1285 KINAOK = KEND1 1286 KSCR1 = KINAOK + NNBST(ISALBE)*NRHFFR(ISYMK) 1287 KSCR2 = KSCR1 + N2BST(ISALBE) 1288 KEND2 = KSCR2 + NF1FRO(ISYMCJ) 1289 LWRK2 = LWORK - KEND2 1290C 1291 IF (LWRK2 .LT. 0) THEN 1292 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND2 1293 CALL QUIT('Insufficient space for allocation in CC_FRCR12') 1294 ENDIF 1295C 1296 CALL DZERO(WORK(KINAOK),NNBST(ISALBE)*NRHFFR(ISYMK)) 1297C 1298C-------------------------------------------- 1299C Transform gamma-index to frozen. 1300C-------------------------------------------- 1301C 1302 KOFF2 = IDSAOG(ISYMG,ISYDIS) + 1 1303 KOFF3 = ILRHSI(ISYMG) + 1 1304C 1305 NTOTAB = MAX(NNBST(ISALBE),1) 1306 NTOTG = MAX(NBAS(ISYMG),1) 1307C 1308 CALL DGEMM('N','N',NNBST(ISALBE),NRHFFR(ISYMK),NBAS(ISYMG),ONE, 1309 * XINT(KOFF2),NTOTAB,CMO(KOFF3),NTOTG,ZERO, 1310 * WORK(KINAOK),NTOTAB) 1311C 1312 DO 110 K = 1,NRHFFR(ISYMK) 1313C 1314 CALL DZERO(WORK(KSCR1),N2BST(ISALBE)) 1315 CALL DZERO(WORK(KSCR2),NF1FRO(ISYMCJ)) 1316C 1317C----------------------------------- 1318C Square up the integrals. 1319C----------------------------------- 1320C 1321 KOFF4 = KINAOK + NNBST(ISALBE)*(K - 1) 1322 CALL CCSD_SYMSQ(WORK(KOFF4),ISALBE,WORK(KSCR1)) 1323C 1324C--------------------------------------------------------------- 1325C Inner symmetry loop and final work space allocation. 1326C--------------------------------------------------------------- 1327C 1328 DO 120 ISYMAL = 1,NSYM 1329C 1330 ISYMI = ISYMAL 1331 ISYMBE = MULD2H(ISALBE,ISYMAL) 1332 ISYMJ = ISYMBE 1333C 1334 KSCR3 = KEND2 1335 KEND3 = KSCR3 + NBAS(ISYMAL)*NRHF(ISYMJ) 1336 LWRK3 = LWORK - KEND3 1337C 1338 IF (LWRK3 .LT. 0) THEN 1339 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND3 1340 CALL QUIT('Insufficient space for allocation '// 1341 & 'in CC_FRCR12') 1342 ENDIF 1343C 1344C------------------------------------------------ 1345C Construct the integrals (iJ|Kd). 1346C------------------------------------------------ 1347C 1348 KOFF5 = KSCR1 + IAODIS(ISYMAL,ISYMBE) 1349 KOFF6 = ILRHSI(ISYMBE) + NBAS(ISYMBE)*NRHFFR(ISYMJ) + 1 1350C 1351 NTOTAL = MAX(NBAS(ISYMAL),1) 1352 NTOTBE = MAX(NBAS(ISYMBE),1) 1353C 1354 CALL DGEMM('N','N',NBAS(ISYMAL),NRHF(ISYMJ), 1355 * NBAS(ISYMBE),ONE,WORK(KOFF5),NTOTAL, 1356 * CMO(KOFF6),NTOTBE,ZERO,WORK(KSCR3),NTOTAL) 1357C 1358 KOFF7 = ILRHSI(ISYMAL) + 1 1359 KOFF8 = KSCR2 + IF1FRO(ISYMI,ISYMJ) 1360C 1361 NTOTAL = MAX(NBAS(ISYMAL),1) 1362 NTOTI = MAX(NRHFFR(ISYMI),1) 1363C 1364 CALL DGEMM('T','N',NRHFFR(ISYMI),NRHF(ISYMJ), 1365 * NBAS(ISYMAL),ONE,CMO(KOFF7),NTOTAL, 1366 * WORK(KSCR3),NTOTAL,ZERO,WORK(KOFF8),NTOTI) 1367C 1368 120 CONTINUE 1369C 1370C---------------------------------------------------------------- 1371C Final scaling with CMO element and storage in result. 1372C---------------------------------------------------------------- 1373C 1374 DO 130 D = 1,NVIR(ISYMD) 1375 NDK = IT1FRO(ISYMD,ISYMK) + NVIR(ISYMD)*(K - 1) + D 1376 KOFF9 = KDVEC + D - 1 1377 KOFF10 = IF2FRO(ISYMCJ,ISYMDK) + NDK 1378 CALL DAXPY(NF1FRO(ISYMCJ),WORK(KOFF9),WORK(KSCR2),1, 1379 * XFRIN(KOFF10),NT1FRO(ISYMDK)) 1380 130 CONTINUE 1381 110 CONTINUE 1382 100 CONTINUE 1383C 1384 RETURN 1385 END 1386C /* Deck fock_all */ 1387 SUBROUTINE FOCK_ALL(FOCKD,WORK,LWORK) 1388C 1389C Written by Asger Halkier 25/5 - 1998 1390C 1391C To read in and change the symmetry ordering of the FULL 1392C Fock matrix diagonal from Sirius to CC-ordering. 1393C 1394#include "implicit.h" 1395#include "dummy.h" 1396 DIMENSION FOCKD(*),WORK(LWORK) 1397#include "priunit.h" 1398#include "inftap.h" 1399#include "ccorb.h" 1400#include "ccsdinp.h" 1401#include "ccsdsym.h" 1402C 1403 IF (LWORK .LT. NORBTS) THEN 1404 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', NORBTS 1405 CALL QUIT('Insufficient work space for allocation in '// 1406 & 'FOCK_ALL') 1407 ENDIF 1408C 1409C------------------------------------- 1410C Read canonical orbital energies. 1411C------------------------------------- 1412C 1413 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY, 1414 & .FALSE.) 1415 REWIND (LUSIFC) 1416C 1417 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 1418 READ (LUSIFC) 1419 READ (LUSIFC) (WORK(I), I = 1,NORBTS) 1420C 1421 CALL GPCLOSE(LUSIFC,'KEEP') 1422C 1423C------------------------------ 1424C Do the actual reordering. 1425C------------------------------ 1426C 1427 IDRHF = 0 1428 IDVIR = NRHFTS 1429 ICOUNT = 0 1430C 1431 DO 100 ISYM = 1,NSYM 1432 DO 110 I = 1,NRHFS(ISYM) 1433 IDRHF = IDRHF + 1 1434 ICOUNT = ICOUNT + 1 1435 FOCKD(IDRHF) = WORK(ICOUNT) 1436 110 CONTINUE 1437C 1438 DO 120 A = 1,NVIRS(ISYM) 1439 IDVIR = IDVIR + 1 1440 ICOUNT = ICOUNT + 1 1441 FOCKD(IDVIR) = WORK(ICOUNT) 1442 120 CONTINUE 1443 100 CONTINUE 1444C 1445 IF (IPRINT .GT. 20) THEN 1446 CALL AROUND('Fock matrix diagonal in FOCK_ALL') 1447 WRITE(LUPRI,1) 1448 DO 200 I = 1,NORBTS 1449 WRITE(LUPRI,2) WORK(I), FOCKD(I) 1450 200 CONTINUE 1451 END IF 1452C 1453 RETURN 1454C 1455 1 FORMAT(7X,'Sirius order',5X,'CC order') 1456 2 FORMAT(6X,F14.10,3X,F14.10) 1457C 1458 END 1459C /* Deck cc_d1fcb */ 1460 SUBROUTINE CC_D1FCB(AODEN,DIJ,DJI,DAI,DIA,WORK,LWORK,ISDEN,ICON) 1461C 1462C Written by Asger Halkier 27/5 - 1998. 1463C 1464C To calculate the contributions to the AO one electron 1465C density AODEN from the frozen core blocks. 1466C 1467#include "implicit.h" 1468 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0) 1469 DIMENSION AODEN(*), DIJ(*), DJI(*), DAI(*), DIA(*), WORK(LWORK) 1470#include "priunit.h" 1471#include "maxorb.h" 1472#include "ccsdinp.h" 1473#include "ccorb.h" 1474#include "ccsdsym.h" 1475#include "ccfro.h" 1476C 1477 CHARACTER MODEL*5 1478C 1479C------------------------------- 1480C Work space allocation one. 1481C------------------------------- 1482C 1483 KCMO = 1 1484 KEND1 = KCMO + NLAMDS 1485 LWRK1 = LWORK - KEND1 1486C 1487 IF (LWRK1 .LT. 0) THEN 1488 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 1489 CALL QUIT('Insufficient work space for allocation '// 1490 & 'in CC_D1FCB') 1491 ENDIF 1492C 1493C---------------------------------------- 1494C Get the FULL MO coefficient matrix. 1495C---------------------------------------- 1496C 1497 CALL CMO_ALL(WORK(KCMO),WORK(KEND1),LWRK1) 1498C 1499C----------------------------------------------- 1500C First symmetry loop and second allocation. 1501C----------------------------------------------- 1502C 1503 DO 100 ISYMAL = 1,NSYM 1504C 1505 ISYMBE = MULD2H(ISDEN,ISYMAL) 1506 ISYMI1 = ISYMAL 1507 ISYMJ1 = ISYMBE 1508 ISYMI2 = ISYMBE 1509 ISYMJ2 = ISYMAL 1510C 1511 KSCR1 = KEND1 1512 KSCR2 = KSCR1 + NBAS(ISYMAL)*NRHFFR(ISYMJ1) 1513 KEND2 = KSCR2 + NBAS(ISYMBE)*NRHFFR(ISYMJ2) 1514 LWRK2 = LWORK - KEND2 1515C 1516 IF (LWRK2 .LT. 0) THEN 1517 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND2 1518 CALL QUIT('Insufficient work space for allocation '// 1519 & 'in CC_D1FCB') 1520 ENDIF 1521C 1522 CALL DZERO(WORK(KSCR1),KEND2-KEND1) 1523C 1524C----------------------------------------------------- 1525C Calculate the contribution from the iJ block. 1526C----------------------------------------------------- 1527C 1528 KOFF1 = KCMO + ILRHSI(ISYMAL) + NBAS(ISYMAL)*NRHFFR(ISYMI1) 1529 KOFF2 = ICOFRO(ISYMI1,ISYMJ1) + 1 1530C 1531 NTOTAL = MAX(NBAS(ISYMAL),1) 1532 NTOTI = MAX(NRHF(ISYMI1),1) 1533C 1534 CALL DGEMM('N','N',NBAS(ISYMAL),NRHFFR(ISYMJ1),NRHF(ISYMI1), 1535 * ONE,WORK(KOFF1),NTOTAL,DIJ(KOFF2),NTOTI,ZERO, 1536 * WORK(KSCR1),NTOTAL) 1537C 1538 KOFF3 = KCMO + ILRHSI(ISYMBE) 1539 KOFF4 = IAODIS(ISYMAL,ISYMBE) + 1 1540C 1541 NTOTAL = MAX(NBAS(ISYMAL),1) 1542 NTOTBE = MAX(NBAS(ISYMBE),1) 1543C 1544 CALL DGEMM('N','T',NBAS(ISYMAL),NBAS(ISYMBE),NRHFFR(ISYMJ1), 1545 * ONE,WORK(KSCR1),NTOTAL,WORK(KOFF3),NTOTBE,ONE, 1546 * AODEN(KOFF4),NTOTAL) 1547C 1548C----------------------------------------------------- 1549C Calculate the contribution from the Ji block. 1550C----------------------------------------------------- 1551C 1552 KOFF5 = KCMO + ILRHSI(ISYMBE) + NBAS(ISYMBE)*NRHFFR(ISYMI2) 1553 KOFF6 = ICOFRO(ISYMI2,ISYMJ2) + 1 1554C 1555 NTOTBE = MAX(NBAS(ISYMBE),1) 1556 NTOTI = MAX(NRHF(ISYMI2),1) 1557C 1558 CALL DGEMM('N','N',NBAS(ISYMBE),NRHFFR(ISYMJ2),NRHF(ISYMI2), 1559 * ONE,WORK(KOFF5),NTOTBE,DJI(KOFF6),NTOTI,ZERO, 1560 * WORK(KSCR2),NTOTBE) 1561C 1562 KOFF7 = KCMO + ILRHSI(ISYMAL) 1563 KOFF8 = IAODIS(ISYMAL,ISYMBE) + 1 1564C 1565 NTOTAL = MAX(NBAS(ISYMAL),1) 1566 NTOTBE = MAX(NBAS(ISYMBE),1) 1567C 1568 CALL DGEMM('N','T',NBAS(ISYMAL),NBAS(ISYMBE),NRHFFR(ISYMJ2), 1569 * ONE,WORK(KOFF7),NTOTAL,WORK(KSCR2),NTOTBE,ONE, 1570 * AODEN(KOFF8),NTOTAL) 1571C 1572 100 CONTINUE 1573C 1574C----------------------------------------------- 1575C Second symmetry loop and final allocation. 1576C----------------------------------------------- 1577C 1578 DO 110 ISYMAL = 1,NSYM 1579C 1580 ISYMBE = MULD2H(ISDEN,ISYMAL) 1581 ISYMA1 = ISYMAL 1582 ISYMI1 = ISYMBE 1583 ISYMA2 = ISYMBE 1584 ISYMI2 = ISYMAL 1585C 1586 KSCR1 = KEND1 1587 KSCR2 = KSCR1 + NBAS(ISYMAL)*NRHFFR(ISYMI1) 1588 KEND3 = KSCR2 + NBAS(ISYMBE)*NRHFFR(ISYMI2) 1589 LWRK3 = LWORK - KEND3 1590C 1591 IF (LWRK3 .LT. 0) THEN 1592 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND3 1593 CALL QUIT('Insufficient work space for allocation '// 1594 & 'in CC_D1FCB') 1595 ENDIF 1596C 1597 CALL DZERO(WORK(KSCR1),KEND3-KEND1) 1598C 1599C----------------------------------------------------- 1600C Calculate the contribution from the aI block. 1601C----------------------------------------------------- 1602C 1603 KOFF9 = KCMO + ILVISI(ISYMAL) 1604 KOFF10 = IT1FRO(ISYMA1,ISYMI1) + 1 1605C 1606 NTOTAL = MAX(NBAS(ISYMAL),1) 1607 NTOTA = MAX(NVIR(ISYMA1),1) 1608C 1609 CALL DGEMM('N','N',NBAS(ISYMAL),NRHFFR(ISYMI1),NVIR(ISYMA1), 1610 * ONE,WORK(KOFF9),NTOTAL,DAI(KOFF10),NTOTA,ZERO, 1611 * WORK(KSCR1),NTOTAL) 1612C 1613 KOFF11 = KCMO + ILRHSI(ISYMBE) 1614 KOFF12 = IAODIS(ISYMAL,ISYMBE) + 1 1615C 1616 NTOTAL = MAX(NBAS(ISYMAL),1) 1617 NTOTBE = MAX(NBAS(ISYMBE),1) 1618C 1619 CALL DGEMM('N','T',NBAS(ISYMAL),NBAS(ISYMBE),NRHFFR(ISYMI1), 1620 * ONE,WORK(KSCR1),NTOTAL,WORK(KOFF11),NTOTBE,ONE, 1621 * AODEN(KOFF12),NTOTAL) 1622C 1623C----------------------------------------------------- 1624C Calculate the contribution from the Ia block. 1625C----------------------------------------------------- 1626C 1627 KOFF13 = KCMO + ILVISI(ISYMBE) 1628 KOFF14 = IT1FRO(ISYMA2,ISYMI2) + 1 1629C 1630 NTOTBE = MAX(NBAS(ISYMBE),1) 1631 NTOTA = MAX(NVIR(ISYMA2),1) 1632C 1633 CALL DGEMM('N','N',NBAS(ISYMBE),NRHFFR(ISYMI2),NVIR(ISYMA2), 1634 * ONE,WORK(KOFF13),NTOTBE,DIA(KOFF14),NTOTA,ZERO, 1635 * WORK(KSCR2),NTOTBE) 1636C 1637 KOFF15 = KCMO + ILRHSI(ISYMAL) 1638 KOFF16 = IAODIS(ISYMAL,ISYMBE) + 1 1639C 1640 NTOTAL = MAX(NBAS(ISYMAL),1) 1641 NTOTBE = MAX(NBAS(ISYMBE),1) 1642C 1643 CALL DGEMM('N','T',NBAS(ISYMAL),NBAS(ISYMBE),NRHFFR(ISYMI2), 1644 * ONE,WORK(KOFF15),NTOTAL,WORK(KSCR2),NTOTBE,ONE, 1645 * AODEN(KOFF16),NTOTAL) 1646C 1647 110 CONTINUE 1648C 1649C---------------------------------------- 1650C Add terms from frozen-frozen block. 1651C---------------------------------------- 1652C 1653 IF (ICON .EQ. 1) THEN 1654 MODEL = 'DUMMY' 1655 CALL CC_FCD1AO(AODEN,WORK(KEND1),LWRK1,MODEL) 1656 ENDIF 1657C 1658 RETURN 1659 END 1660C /* Deck cc_gtofro */ 1661 SUBROUTINE CC_GTOFRO(XINT,DSFRO,CMO,WORK,LWORK,ISYDIS) 1662C 1663C Written by Asger Halkier 28/5 - 1998. 1664C 1665C To calculate the integral batch (al be | fro del). 1666C 1667#include "implicit.h" 1668 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0) 1669 DIMENSION XINT(*), DSFRO(*), CMO(*), WORK(LWORK) 1670#include "priunit.h" 1671#include "maxorb.h" 1672#include "ccsdinp.h" 1673#include "ccorb.h" 1674#include "ccsdsym.h" 1675#include "ccfro.h" 1676C 1677 DO 100 ISYMG = 1,NSYM 1678C 1679 ISYMI = ISYMG 1680 ISALBE = MULD2H(ISYMG,ISYDIS) 1681C 1682 KOFF1 = IDSAOG(ISYMG,ISYDIS) + 1 1683 KOFF2 = ILRHSI(ISYMG) + 1 1684 KOFF3 = IDSFRO(ISALBE,ISYMI) + 1 1685C 1686 NTOTAB = MAX(NNBST(ISALBE),1) 1687 NTOTG = MAX(NBAS(ISYMG),1) 1688C 1689 CALL DGEMM('N','N',NNBST(ISALBE),NRHFFR(ISYMI),NBAS(ISYMG), 1690 * ONE,XINT(KOFF1),NTOTAB,CMO(KOFF2),NTOTG,ZERO, 1691 * DSFRO(KOFF3),NTOTAB) 1692C 1693 100 CONTINUE 1694C 1695 RETURN 1696 END 1697C /* Deck cc_ofroin */ 1698 SUBROUTINE CC_OFROIN(DSRHF,DSRES,CMO,WORK,LWORK,ISYDIS) 1699C 1700C Written by Asger Halkier 28/5 - 1998. 1701C 1702C To calculate the integral batch (cor fro | cor del). 1703C 1704#include "implicit.h" 1705 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0) 1706 DIMENSION DSRHF(*), DSRES(*), CMO(*), WORK(LWORK) 1707#include "priunit.h" 1708#include "maxorb.h" 1709#include "ccsdinp.h" 1710#include "ccorb.h" 1711#include "ccsdsym.h" 1712#include "ccfro.h" 1713C 1714C------------------------------------------------------- 1715C Outer symmetry loop and work space allocation one. 1716C------------------------------------------------------- 1717C 1718 DO 100 ISYML = 1,NSYM 1719C 1720 ISALBE = MULD2H(ISYDIS,ISYML) 1721 ISYMKI = MULD2H(ISYDIS,ISYML) 1722C 1723C---------------------------------- 1724C Work space allocation one. 1725C---------------------------------- 1726C 1727 KAOINT = 1 1728 KEND1 = KAOINT + N2BST(ISALBE) 1729 LWRK1 = LWORK - KEND1 1730C 1731 IF (LWRK1 .LT. 0) THEN 1732 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 1733 CALL QUIT('Insufficient memory for allocation '// 1734 & 'in CC_OFROIN') 1735 ENDIF 1736C 1737 CALL DZERO(WORK(KAOINT),N2BST(ISALBE)) 1738C 1739 DO 110 L = 1,NRHF(ISYML) 1740C 1741C---------------------------------------- 1742C Unpack integral distribution. 1743C---------------------------------------- 1744C 1745 KOFF1 = IDSRHF(ISALBE,ISYML) + NNBST(ISALBE)*(L - 1) + 1 1746C 1747 CALL CCSD_SYMSQ(DSRHF(KOFF1),ISALBE,WORK(KAOINT)) 1748C 1749C------------------------------------------------------------- 1750C Inner symmetry loop and work space allocation two. 1751C------------------------------------------------------------- 1752C 1753 DO 120 ISYMAL = 1,NSYM 1754C 1755 ISYMBE = MULD2H(ISYMAL,ISALBE) 1756 ISYMK = ISYMAL 1757 ISYMI = ISYMBE 1758C 1759 KSCR1 = KEND1 1760 KEND2 = KSCR1 + NBAS(ISYMAL)*NRHFFR(ISYMI) 1761 LWRK2 = LWORK - KEND2 1762C 1763 IF (LWRK2 .LT. 0) THEN 1764 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND2 1765 CALL QUIT('Insufficient space for allocation '// 1766 & 'in CC_OFROIN') 1767 ENDIF 1768C 1769 CALL DZERO(WORK(KSCR1),NBAS(ISYMAL)*NRHFFR(ISYMI)) 1770C 1771C---------------------------------------------------------------- 1772C Perform the contractions for obtaining the result. 1773C---------------------------------------------------------------- 1774C 1775 KOFF2 = KAOINT + IAODIS(ISYMAL,ISYMBE) 1776 KOFF3 = ILRHSI(ISYMBE) + 1 1777C 1778 NTOTAL = MAX(NBAS(ISYMAL),1) 1779 NTOTBE = MAX(NBAS(ISYMBE),1) 1780C 1781 CALL DGEMM('N','N',NBAS(ISYMAL),NRHFFR(ISYMI), 1782 * NBAS(ISYMBE),ONE,WORK(KOFF2),NTOTAL, 1783 * CMO(KOFF3),NTOTBE,ZERO,WORK(KSCR1),NTOTAL) 1784C 1785 KOFF4 = ILRHSI(ISYMAL) + NBAS(ISYMAL)*NRHFFR(ISYMK) + 1 1786 KOFF5 = IOFROO(ISYMKI,ISYML) + NCOFRO(ISYMKI)*(L - 1) 1787 * + ICOFRO(ISYMK,ISYMI) + 1 1788C 1789 NTOTAL = MAX(NBAS(ISYMAL),1) 1790 NTOTK = MAX(NRHF(ISYMK),1) 1791C 1792 CALL DGEMM('T','N',NRHF(ISYMK),NRHFFR(ISYMI), 1793 * NBAS(ISYMAL),ONE,CMO(KOFF4),NTOTAL, 1794 * WORK(KSCR1),NTOTAL,ZERO,DSRES(KOFF5),NTOTK) 1795C 1796 120 CONTINUE 1797 110 CONTINUE 1798 100 CONTINUE 1799C 1800 RETURN 1801 END 1802C /* Deck mp2_etfrdi */ 1803 SUBROUTINE MP2_ETFRD(ETAFRO,DSINT,TINT,ISYDIS) 1804C 1805C Written by Asger Halkier 28/5 - 1998. 1806C 1807C To calculate the direct contribution to the frozen part of 1808C eta in MP2 calculations. 1809C 1810#include "implicit.h" 1811 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0) 1812 DIMENSION ETAFRO(*), DSINT(*), TINT(*) 1813#include "priunit.h" 1814#include "maxorb.h" 1815#include "ccsdinp.h" 1816#include "ccorb.h" 1817#include "ccsdsym.h" 1818#include "ccfro.h" 1819C 1820C---------------------------- 1821C Calculate contribution. 1822C---------------------------- 1823C 1824 DO 100 ISYML = 1,NSYM 1825C 1826 ISYMKI = MULD2H(ISYDIS,ISYML) 1827 ISYMAK = MULD2H(ISYDIS,ISYML) 1828C 1829 DO 110 L = 1,NRHF(ISYML) 1830 DO 120 ISYMA = 1,NSYM 1831C 1832 ISYMK = MULD2H(ISYMAK,ISYMA) 1833 ISYMI = ISYMA 1834C 1835 KOFF1 = IT2BCD(ISYMAK,ISYML) + NT1AM(ISYMAK)*(L - 1) 1836 * + IT1AM(ISYMA,ISYMK) + 1 1837 KOFF2 = IOFROO(ISYMKI,ISYML) + NCOFRO(ISYMKI)*(L - 1) 1838 * + ICOFRO(ISYMK,ISYMI) + 1 1839 KOFF3 = IT1FRO(ISYMA,ISYMI) + 1 1840C 1841 NTOTA = MAX(NVIR(ISYMA),1) 1842 NTOTK = MAX(NRHF(ISYMK),1) 1843C 1844 CALL DGEMM('N','N',NVIR(ISYMA),NRHFFR(ISYMI), 1845 * NRHF(ISYMK),-ONE,TINT(KOFF1),NTOTA, 1846 * DSINT(KOFF2),NTOTK,ONE,ETAFRO(KOFF3),NTOTA) 1847C 1848 120 CONTINUE 1849 110 CONTINUE 1850 100 CONTINUE 1851C 1852 RETURN 1853 END 1854C /* Deck mp2_eidv1 */ 1855 SUBROUTINE MP2_EIDV1(ETAFRO,XINTFR,ZKAB,CMO,WORK,LWORK,IDEL, 1856 * ISYDEL) 1857C 1858C Written by Asger Halkier 29/5 - 1998. 1859C 1860C To calculate the coulomb part of the indirect virtual 1861C contribution to the frozen part of eta in MP2 calculations. 1862C 1863#include "implicit.h" 1864 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, FOUR = 4.0D0) 1865 DIMENSION ETAFRO(*), XINTFR(*), ZKAB(*), CMO(*), WORK(LWORK) 1866#include "priunit.h" 1867#include "maxorb.h" 1868#include "ccsdinp.h" 1869#include "ccorb.h" 1870#include "ccsdsym.h" 1871#include "ccfro.h" 1872C 1873C-------------------------------------------------- 1874C Initial symmetries and work space allocation. 1875C-------------------------------------------------- 1876C 1877 ISYMA = ISYDEL 1878 ISYMI = ISYMA 1879 ISYMBC = 1 1880 ISALBE = ISYMBC 1881C 1882 IF (NRHFFR(ISYMI) .EQ. 0) RETURN 1883C 1884 KAVEC = 1 1885 KAOINT = KAVEC + NVIR(ISYMA) 1886 KMOINT = KAOINT + N2BST(ISALBE) 1887 KEND1 = KMOINT + NMATAB(ISYMBC) 1888 LWRK1 = LWORK - KEND1 1889C 1890 IF (LWRK1 .LT. 0) THEN 1891 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 1892 CALL QUIT('Insufficient work space for allocation '// 1893 & 'in MP2_EIDV1') 1894 ENDIF 1895C 1896 CALL DZERO(WORK(KAVEC),NVIR(ISYMA)) 1897 CALL DZERO(WORK(KAOINT),N2BST(ISALBE)) 1898C 1899C----------------------------------- 1900C Copy vector out of CMO matrix. 1901C----------------------------------- 1902C 1903 KOFF1 = ILVISI(ISYDEL) + IDEL - IBAS(ISYDEL) 1904C 1905 CALL DCOPY(NVIR(ISYMA),CMO(KOFF1),NBAS(ISYDEL),WORK(KAVEC),1) 1906C 1907 DO 100 I = 1,NRHFFR(ISYMI) 1908C 1909 CALL DZERO(WORK(KMOINT),NMATAB(ISYMBC)) 1910C 1911C-------------------------------------------------------- 1912C Square up integral distribution (al be | I del). 1913C-------------------------------------------------------- 1914C 1915 KOFF2 = IDSFRO(ISALBE,ISYMI) + NNBST(ISALBE)*(I - 1) + 1 1916C 1917 CALL CCSD_SYMSQ(XINTFR(KOFF2),ISALBE,WORK(KAOINT)) 1918C 1919C------------------------------------------------------ 1920C Inner symmetry loop and final work allocation. 1921C------------------------------------------------------ 1922C 1923 DO 110 ISYMB = 1,NSYM 1924C 1925 ISYMC = MULD2H(ISYMBC,ISYMB) 1926 ISYMAL = ISYMB 1927 ISYMBE = ISYMC 1928C 1929 KSCR1 = KEND1 1930 KEND2 = KSCR1 + NBAS(ISYMAL)*NVIR(ISYMC) 1931 LWRK2 = LWORK - KEND2 1932C 1933 IF (LWRK2 .LT. 0) THEN 1934 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND2 1935 CALL QUIT('Insufficient memory for allocation '// 1936 & 'in MP2_EIDV1') 1937 ENDIF 1938C 1939 CALL DZERO(WORK(KSCR1),NBAS(ISYMAL)*NVIR(ISYMC)) 1940C 1941C--------------------------------------------- 1942C Calculate integrals (b c | I del). 1943C--------------------------------------------- 1944C 1945 KOFF3 = KAOINT + IAODIS(ISYMAL,ISYMBE) 1946 KOFF4 = ILVISI(ISYMBE) + 1 1947C 1948 NTOTAL = MAX(NBAS(ISYMAL),1) 1949 NTOTBE = MAX(NBAS(ISYMBE),1) 1950C 1951 CALL DGEMM('N','N',NBAS(ISYMAL),NVIR(ISYMC),NBAS(ISYMBE), 1952 * ONE,WORK(KOFF3),NTOTAL,CMO(KOFF4),NTOTBE,ZERO, 1953 * WORK(KSCR1),NTOTAL) 1954C 1955 KOFF5 = ILVISI(ISYMAL) + 1 1956 KOFF6 = KMOINT + IMATAB(ISYMB,ISYMC) 1957C 1958 NTOTAL = MAX(NBAS(ISYMAL),1) 1959 NTOTB = MAX(NVIR(ISYMB),1) 1960C 1961 CALL DGEMM('T','N',NVIR(ISYMB),NVIR(ISYMC),NBAS(ISYMAL), 1962 * ONE,CMO(KOFF5),NTOTAL,WORK(KSCR1),NTOTAL,ZERO, 1963 * WORK(KOFF6),NTOTB) 1964C 1965 110 CONTINUE 1966C 1967C----------------------------------- 1968C Calculate the contribution. 1969C----------------------------------- 1970C 1971 FACT = FOUR*DDOT(NMATAB(ISYMBC),WORK(KMOINT),1,ZKAB(1),1) 1972 KOFF7 = IT1FRO(ISYMA,ISYMI) + NVIR(ISYMA)*(I - 1) + 1 1973 CALL DAXPY(NVIR(ISYMA),FACT,WORK(KAVEC),1,ETAFRO(KOFF7),1) 1974C 1975 100 CONTINUE 1976C 1977 RETURN 1978 END 1979C /* Deck mp2_eidv2 */ 1980 SUBROUTINE MP2_EIDV2(ETAFRO,XINTFR,ZKAB,CMO,WORK,LWORK,IDEL, 1981 * ISYDEL) 1982C 1983C Written by Asger Halkier 29/5 - 1998. 1984C 1985C To calculate the exchange part of the indirect virtual 1986C contribution to the frozen part of eta in MP2 calculations. 1987C 1988#include "implicit.h" 1989 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 1990 DIMENSION ETAFRO(*), XINTFR(*), ZKAB(*), CMO(*), WORK(LWORK) 1991#include "priunit.h" 1992#include "maxorb.h" 1993#include "ccsdinp.h" 1994#include "ccorb.h" 1995#include "ccsdsym.h" 1996#include "ccfro.h" 1997C 1998C------------------------------------------------ 1999C Initial symmetries & work space allocation. 2000C------------------------------------------------ 2001C 2002 ISYMAI = 1 2003 ISYMBC = 1 2004 ISYMC = ISYDEL 2005 ISYMB = MULD2H(ISYMBC,ISYMC) 2006C 2007 KCVEC = 1 2008 KBVEC = KCVEC + NVIR(ISYMC) 2009 KEND1 = KBVEC + NVIR(ISYMB) 2010 LWRK1 = LWORK - KEND1 2011C 2012 IF (LWRK1 .LT. 0) THEN 2013 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 2014 CALL QUIT('Insufficient work space for allocation '// 2015 & 'in MP2_EIDV2') 2016 ENDIF 2017C 2018 CALL DZERO(WORK(KCVEC),NVIR(ISYMC)) 2019 CALL DZERO(WORK(KBVEC),NVIR(ISYMB)) 2020C 2021C----------------------------------- 2022C Copy vector out of CMO matrix. 2023C----------------------------------- 2024C 2025 KOFFC = ILVISI(ISYDEL) + IDEL - IBAS(ISYDEL) 2026C 2027 CALL DCOPY(NVIR(ISYMC),CMO(KOFFC),NBAS(ISYDEL),WORK(KCVEC),1) 2028C 2029C----------------------------------- 2030C Contract ZKAB with CMO vector. 2031C----------------------------------- 2032C 2033 KOFF1 = IMATAB(ISYMB,ISYMC) + 1 2034C 2035 NTOTB = MAX(NVIR(ISYMB),1) 2036C 2037 CALL DGEMV('N',NVIR(ISYMB),NVIR(ISYMC),ONE,ZKAB(KOFF1),NTOTB, 2038 * WORK(KCVEC),1,ZERO,WORK(KBVEC),1) 2039C 2040C---------------------------------------------------- 2041C Symmetry loop and second work space allocation. 2042C---------------------------------------------------- 2043C 2044 DO 100 ISYMI = 1,NSYM 2045C 2046 IF (NRHFFR(ISYMI) .EQ. 0) GOTO 100 2047C 2048 ISYMA = MULD2H(ISYMAI,ISYMI) 2049 ISYMAL = ISYMA 2050 ISYMBE = ISYMB 2051 ISALBE = MULD2H(ISYMAL,ISYMBE) 2052C 2053 KAOINT = KEND1 2054 KSCR1 = KAOINT + N2BST(ISALBE) 2055 KSCR2 = KSCR1 + NBAS(ISYMAL)*NVIR(ISYMB) 2056 KEND2 = KSCR2 + NVIR(ISYMA)*NVIR(ISYMB) 2057 LWRK2 = LWORK - KEND2 2058C 2059 IF (LWRK2 .LT. 0) THEN 2060 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND2 2061 CALL QUIT('Insufficient work space for allocation '// 2062 & 'in MP2_EIDV2') 2063 ENDIF 2064C 2065 CALL DZERO(WORK(KAOINT),N2BST(ISALBE)) 2066 CALL DZERO(WORK(KSCR1),NBAS(ISYMAL)*NVIR(ISYMB)) 2067 CALL DZERO(WORK(KSCR2),NVIR(ISYMA)*NVIR(ISYMB)) 2068C 2069 DO 110 I = 1,NRHFFR(ISYMI) 2070C 2071C---------------------------------------- 2072C Unpack integral distribution. 2073C---------------------------------------- 2074C 2075 KOFF2 = IDSFRO(ISALBE,ISYMI) + NNBST(ISALBE)*(I - 1) + 1 2076C 2077 CALL CCSD_SYMSQ(XINTFR(KOFF2),ISALBE,WORK(KAOINT)) 2078C 2079C------------------------------------------- 2080C Transform integrals to MO basis. 2081C------------------------------------------- 2082C 2083 KOFF3 = KAOINT + IAODIS(ISYMAL,ISYMBE) 2084 KOFF4 = ILVISI(ISYMBE) + 1 2085C 2086 NTOTAL = MAX(NBAS(ISYMAL),1) 2087 NTOTBE = MAX(NBAS(ISYMBE),1) 2088C 2089 CALL DGEMM('N','N',NBAS(ISYMAL),NVIR(ISYMB),NBAS(ISYMBE), 2090 * ONE,WORK(KOFF3),NTOTAL,CMO(KOFF4),NTOTBE,ZERO, 2091 * WORK(KSCR1),NTOTAL) 2092C 2093 KOFF5 = ILVISI(ISYMAL) + 1 2094C 2095 NTOTAL = MAX(NBAS(ISYMAL),1) 2096 NTOTA = MAX(NVIR(ISYMA),1) 2097C 2098 CALL DGEMM('T','N',NVIR(ISYMA),NVIR(ISYMB),NBAS(ISYMAL), 2099 * ONE,CMO(KOFF5),NTOTAL,WORK(KSCR1),NTOTAL,ZERO, 2100 * WORK(KSCR2),NTOTA) 2101C 2102C-------------------------------------------------- 2103C Final matrix vector product for result. 2104C-------------------------------------------------- 2105C 2106 KOFF6 = IT1FRO(ISYMA,ISYMI) + NVIR(ISYMA)*(I - 1) + 1 2107C 2108 NTOTA = MAX(NVIR(ISYMA),1) 2109C 2110 CALL DGEMV('N',NVIR(ISYMA),NVIR(ISYMB),-TWO,WORK(KSCR2), 2111 * NTOTA,WORK(KBVEC),1,ONE,ETAFRO(KOFF6),1) 2112C 2113 110 CONTINUE 2114 100 CONTINUE 2115C 2116 RETURN 2117 END 2118C /* Deck mp2_eidc1 */ 2119 SUBROUTINE MP2_EIDC1(ETAFRO,XINTFR,ZKIJ,CMO,WORK,LWORK,IDEL, 2120 * ISYDEL) 2121C 2122C Written by Asger Halkier 30/5 - 1998. 2123C 2124C To calculate the coulomb part of the indirect correlated 2125C contribution to the frozen part of eta in MP2 calculations. 2126C 2127#include "implicit.h" 2128 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, FOUR = 4.0D0) 2129 DIMENSION ETAFRO(*), XINTFR(*), ZKIJ(*), CMO(*), WORK(LWORK) 2130#include "priunit.h" 2131#include "maxorb.h" 2132#include "ccsdinp.h" 2133#include "ccorb.h" 2134#include "ccsdsym.h" 2135#include "ccfro.h" 2136C 2137C-------------------------------------------------- 2138C Initial symmetries and work space allocation. 2139C-------------------------------------------------- 2140C 2141 ISYMA = ISYDEL 2142 ISYMI = ISYMA 2143 ISYMJK = 1 2144 ISALBE = ISYMJK 2145C 2146 IF (NRHFFR(ISYMI) .EQ. 0) RETURN 2147C 2148 KAVEC = 1 2149 KAOINT = KAVEC + NVIR(ISYMA) 2150 KMOINT = KAOINT + N2BST(ISALBE) 2151 KEND1 = KMOINT + NMATIJ(ISYMJK) 2152 LWRK1 = LWORK - KEND1 2153C 2154 IF (LWRK1 .LT. 0) THEN 2155 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 2156 CALL QUIT('Insufficient work space for allocation '// 2157 & 'in MP2_EIDC1') 2158 ENDIF 2159C 2160 CALL DZERO(WORK(KAVEC),NVIR(ISYMA)) 2161 CALL DZERO(WORK(KAOINT),N2BST(ISALBE)) 2162C 2163C----------------------------------- 2164C Copy vector out of CMO matrix. 2165C----------------------------------- 2166C 2167 KOFF1 = ILVISI(ISYDEL) + IDEL - IBAS(ISYDEL) 2168C 2169 CALL DCOPY(NVIR(ISYMA),CMO(KOFF1),NBAS(ISYDEL),WORK(KAVEC),1) 2170C 2171 DO 100 I = 1,NRHFFR(ISYMI) 2172C 2173 CALL DZERO(WORK(KMOINT),NMATIJ(ISYMJK)) 2174C 2175C-------------------------------------------------------- 2176C Square up integral distribution (al be | I del). 2177C-------------------------------------------------------- 2178C 2179 KOFF2 = IDSFRO(ISALBE,ISYMI) + NNBST(ISALBE)*(I - 1) + 1 2180C 2181 CALL CCSD_SYMSQ(XINTFR(KOFF2),ISALBE,WORK(KAOINT)) 2182C 2183C------------------------------------------------------ 2184C Inner symmetry loop and final work allocation. 2185C------------------------------------------------------ 2186C 2187 DO 110 ISYMJ = 1,NSYM 2188C 2189 ISYMK = MULD2H(ISYMJK,ISYMJ) 2190 ISYMAL = ISYMJ 2191 ISYMBE = ISYMK 2192C 2193 KSCR1 = KEND1 2194 KEND2 = KSCR1 + NBAS(ISYMAL)*NRHF(ISYMK) 2195 LWRK2 = LWORK - KEND2 2196C 2197 IF (LWRK2 .LT. 0) THEN 2198 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND2 2199 CALL QUIT('Insufficient memory for allocation '// 2200 & 'in MP2_EIDC1') 2201 ENDIF 2202C 2203 CALL DZERO(WORK(KSCR1),NBAS(ISYMAL)*NRHF(ISYMK)) 2204C 2205C--------------------------------------------- 2206C Calculate integrals (j k | I del). 2207C--------------------------------------------- 2208C 2209 KOFF3 = KAOINT + IAODIS(ISYMAL,ISYMBE) 2210 KOFF4 = ILRHSI(ISYMBE) + NBAS(ISYMBE)*NRHFFR(ISYMK) + 1 2211C 2212 NTOTAL = MAX(NBAS(ISYMAL),1) 2213 NTOTBE = MAX(NBAS(ISYMBE),1) 2214C 2215 CALL DGEMM('N','N',NBAS(ISYMAL),NRHF(ISYMK),NBAS(ISYMBE), 2216 * ONE,WORK(KOFF3),NTOTAL,CMO(KOFF4),NTOTBE,ZERO, 2217 * WORK(KSCR1),NTOTAL) 2218C 2219 KOFF5 = ILRHSI(ISYMAL) + NBAS(ISYMAL)*NRHFFR(ISYMJ) + 1 2220 KOFF6 = KMOINT + IMATIJ(ISYMJ,ISYMK) 2221C 2222 NTOTAL = MAX(NBAS(ISYMAL),1) 2223 NTOTJ = MAX(NRHF(ISYMJ),1) 2224C 2225 CALL DGEMM('T','N',NRHF(ISYMJ),NRHF(ISYMK),NBAS(ISYMAL), 2226 * ONE,CMO(KOFF5),NTOTAL,WORK(KSCR1),NTOTAL,ZERO, 2227 * WORK(KOFF6),NTOTJ) 2228C 2229 110 CONTINUE 2230C 2231C----------------------------------- 2232C Calculate the contribution. 2233C----------------------------------- 2234C 2235 FACT = FOUR*DDOT(NMATIJ(ISYMJK),WORK(KMOINT),1,ZKIJ(1),1) 2236 KOFF7 = IT1FRO(ISYMA,ISYMI) + NVIR(ISYMA)*(I - 1) + 1 2237 CALL DAXPY(NVIR(ISYMA),FACT,WORK(KAVEC),1,ETAFRO(KOFF7),1) 2238C 2239 100 CONTINUE 2240C 2241 RETURN 2242 END 2243C /* Deck mp2_eidc2 */ 2244 SUBROUTINE MP2_EIDC2(ETAFRO,XINOFO,ZKIJ,CMO,WORK,LWORK,IDEL, 2245 * ISYDEL) 2246C 2247C Written by Asger Halkier 30/5 - 1998. 2248C 2249C To calculate the exchange part of the indirect correlated 2250C contribution to the frozen part of eta in MP2 calculations. 2251C 2252#include "implicit.h" 2253 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 2254 DIMENSION ETAFRO(*), XINOFO(*), ZKIJ(*), CMO(*), WORK(LWORK) 2255#include "priunit.h" 2256#include "maxorb.h" 2257#include "ccsdinp.h" 2258#include "ccorb.h" 2259#include "ccsdsym.h" 2260#include "ccfro.h" 2261C 2262C-------------------------------------------------- 2263C Initial symmetries and work space allocation. 2264C-------------------------------------------------- 2265C 2266 ISYMA = ISYDEL 2267 ISYMI = ISYMA 2268 ISYMJK = 1 2269C 2270 IF (NRHFFR(ISYMI) .EQ. 0) RETURN 2271C 2272 KAVEC = 1 2273 KIVEC = KAVEC + NVIR(ISYMA) 2274 KEND1 = KIVEC + NRHFFR(ISYMI) 2275 LWRK1 = LWORK - KEND1 2276C 2277 IF (LWRK1 .LT. 0) THEN 2278 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 2279 CALL QUIT('Insufficient work space for allocation '// 2280 & 'in MP2_EIDC2') 2281 ENDIF 2282C 2283 CALL DZERO(WORK(KAVEC),NVIR(ISYMA)) 2284 CALL DZERO(WORK(KIVEC),NRHFFR(ISYMI)) 2285C 2286C----------------------------------- 2287C Copy vector out of CMO matrix. 2288C----------------------------------- 2289C 2290 KOFF1 = ILVISI(ISYDEL) + IDEL - IBAS(ISYDEL) 2291C 2292 CALL DCOPY(NVIR(ISYMA),CMO(KOFF1),NBAS(ISYDEL),WORK(KAVEC),1) 2293C 2294C------------------------------------------------- 2295C Outer loops over third integral batch index. 2296C------------------------------------------------- 2297C 2298 DO 100 ISYMK = 1,NSYM 2299C 2300 ISYMJ = MULD2H(ISYMJK,ISYMK) 2301 ISYMJI = MULD2H(ISYMJ,ISYMI) 2302C 2303 DO 110 K = 1,NRHF(ISYMK) 2304C 2305C------------------------------------------------------- 2306C Inner loop over frozen index and contraction 2307C of integrals and ZKIJ-vector. 2308C------------------------------------------------------- 2309C 2310 DO 120 I = 1,NRHFFR(ISYMI) 2311C 2312 NJK = IMATIJ(ISYMJ,ISYMK) + NRHF(ISYMJ)*(K - 1) + 1 2313 NJIK = IOFROO(ISYMJI,ISYMK) + NCOFRO(ISYMJI)*(K - 1) 2314 * + ICOFRO(ISYMJ,ISYMI) + NRHF(ISYMJ)*(I - 1) + 1 2315C 2316 FACT = DDOT(NRHF(ISYMJ),ZKIJ(NJK),1,XINOFO(NJIK),1) 2317 WORK(KIVEC+I-1) = WORK(KIVEC+I-1) + FACT 2318C 2319 120 CONTINUE 2320 110 CONTINUE 2321 100 CONTINUE 2322C 2323C----------------------------- 2324C Final storage in result. 2325C----------------------------- 2326C 2327 DO 130 I = 1,NRHFFR(ISYMI) 2328C 2329 KOFF2 = KIVEC + I - 1 2330 KOFF3 = IT1FRO(ISYMA,ISYMI) + NVIR(ISYMA)*(I - 1) + 1 2331C 2332 CALL DAXPY(NVIR(ISYMA),-TWO*WORK(KOFF2),WORK(KAVEC),1, 2333 * ETAFRO(KOFF3),1) 2334C 2335 130 CONTINUE 2336C 2337 RETURN 2338 END 2339C /* Deck mp2_eidf1 */ 2340 SUBROUTINE MP2_EIDF1(ETAAI,XINOFO,ZKJK,CMO,WORK,LWORK,IDEL, 2341 * ISYDEL) 2342C 2343C Written by Asger Halkier 30/5 - 1998. 2344C 2345C To calculate the coulomb part of the indirect frozen 2346C contribution to the correlated part of eta in MP2 calculations. 2347C 2348#include "implicit.h" 2349 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, EIGHT = 8.0D0) 2350 DIMENSION ETAAI(*), XINOFO(*), ZKJK(*), CMO(*), WORK(LWORK) 2351#include "priunit.h" 2352#include "maxorb.h" 2353#include "ccsdinp.h" 2354#include "ccorb.h" 2355#include "ccsdsym.h" 2356#include "ccfro.h" 2357C 2358C-------------------------------------------------- 2359C Initial symmetries and work space allocation. 2360C-------------------------------------------------- 2361C 2362 ISYMA = ISYDEL 2363 ISYMI = ISYMA 2364 ISYMJK = 1 2365C 2366 IF (NRHF(ISYMI) .EQ. 0) RETURN 2367C 2368 KAVEC = 1 2369 KEND1 = KAVEC + NVIR(ISYMA) 2370 LWRK1 = LWORK - KEND1 2371C 2372 IF (LWRK1 .LT. 0) THEN 2373 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 2374 CALL QUIT('Insufficient work space for allocation '// 2375 & 'in MP2_EIDF1') 2376 ENDIF 2377C 2378 CALL DZERO(WORK(KAVEC),NVIR(ISYMA)) 2379C 2380C----------------------------------- 2381C Copy vector out of CMO matrix. 2382C----------------------------------- 2383C 2384 KOFF1 = ILVISI(ISYDEL) + IDEL - IBAS(ISYDEL) 2385C 2386 CALL DCOPY(NVIR(ISYMA),CMO(KOFF1),NBAS(ISYDEL),WORK(KAVEC),1) 2387C 2388C---------------------- 2389C Calculate result. 2390C---------------------- 2391C 2392 DO 100 I = 1,NRHF(ISYMI) 2393C 2394 KOFF1 = IOFROO(ISYMJK,ISYMI) + NCOFRO(ISYMJK)*(I - 1) + 1 2395 KOFF2 = IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I - 1) + 1 2396 FACT = EIGHT*DDOT(NCOFRO(ISYMJK),ZKJK(1),1,XINOFO(KOFF1),1) 2397 CALL DAXPY(NVIR(ISYMA),FACT,WORK(KAVEC),1,ETAAI(KOFF2),1) 2398C 2399 100 CONTINUE 2400C 2401 RETURN 2402 END 2403C /* Deck mp2_eidf2 */ 2404 SUBROUTINE MP2_EIDF2(ETAFRO,XINTFR,ZKJK,CMO,WORK,LWORK,IDEL, 2405 * ISYDEL) 2406C 2407C Written by Asger Halkier 30/5 - 1998. 2408C 2409C To calculate the coulomb part of the indirect frozen 2410C contribution to the frozen part of eta in MP2 calculations. 2411C 2412#include "implicit.h" 2413 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, EIGHT = 8.0D0) 2414 DIMENSION ETAFRO(*), XINTFR(*), ZKJK(*), CMO(*), WORK(LWORK) 2415#include "priunit.h" 2416#include "maxorb.h" 2417#include "ccsdinp.h" 2418#include "ccorb.h" 2419#include "ccsdsym.h" 2420#include "ccfro.h" 2421C 2422C-------------------------------------------------- 2423C Initial symmetries and work space allocation. 2424C-------------------------------------------------- 2425C 2426 ISYMA = ISYDEL 2427 ISYMI = ISYMA 2428 ISYMJK = 1 2429 ISALBE = ISYMJK 2430C 2431 IF (NRHFFR(ISYMI) .EQ. 0) RETURN 2432C 2433 KAVEC = 1 2434 KAOINT = KAVEC + NVIR(ISYMA) 2435 KMOINT = KAOINT + N2BST(ISALBE) 2436 KEND1 = KMOINT + NCOFRO(ISYMJK) 2437 LWRK1 = LWORK - KEND1 2438C 2439 IF (LWRK1 .LT. 0) THEN 2440 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 2441 CALL QUIT('Insufficient work space for allocation '// 2442 & 'in MP2_EIDF2') 2443 ENDIF 2444C 2445 CALL DZERO(WORK(KAVEC),NVIR(ISYMA)) 2446 CALL DZERO(WORK(KAOINT),N2BST(ISALBE)) 2447C 2448C----------------------------------- 2449C Copy vector out of CMO matrix. 2450C----------------------------------- 2451C 2452 KOFF1 = ILVISI(ISYDEL) + IDEL - IBAS(ISYDEL) 2453C 2454 CALL DCOPY(NVIR(ISYMA),CMO(KOFF1),NBAS(ISYDEL),WORK(KAVEC),1) 2455C 2456 DO 100 I = 1,NRHFFR(ISYMI) 2457C 2458 CALL DZERO(WORK(KMOINT),NCOFRO(ISYMJK)) 2459C 2460C-------------------------------------------------------- 2461C Square up integral distribution (al be | I del). 2462C-------------------------------------------------------- 2463C 2464 KOFF2 = IDSFRO(ISALBE,ISYMI) + NNBST(ISALBE)*(I - 1) + 1 2465C 2466 CALL CCSD_SYMSQ(XINTFR(KOFF2),ISALBE,WORK(KAOINT)) 2467C 2468C------------------------------------------------------ 2469C Inner symmetry loop and final work allocation. 2470C------------------------------------------------------ 2471C 2472 DO 110 ISYMJ = 1,NSYM 2473C 2474 ISYMK = MULD2H(ISYMJK,ISYMJ) 2475 ISYMAL = ISYMJ 2476 ISYMBE = ISYMK 2477C 2478 KSCR1 = KEND1 2479 KEND2 = KSCR1 + NBAS(ISYMAL)*NRHFFR(ISYMK) 2480 LWRK2 = LWORK - KEND2 2481C 2482 IF (LWRK2 .LT. 0) THEN 2483 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND2 2484 CALL QUIT('Insufficient memory for allocation '// 2485 & 'in MP2_EIDF2') 2486 ENDIF 2487C 2488 CALL DZERO(WORK(KSCR1),NBAS(ISYMAL)*NRHFFR(ISYMK)) 2489C 2490C--------------------------------------------- 2491C Calculate integrals (j K | I del). 2492C--------------------------------------------- 2493C 2494 KOFF3 = KAOINT + IAODIS(ISYMAL,ISYMBE) 2495 KOFF4 = ILRHSI(ISYMBE) + 1 2496C 2497 NTOTAL = MAX(NBAS(ISYMAL),1) 2498 NTOTBE = MAX(NBAS(ISYMBE),1) 2499C 2500 CALL DGEMM('N','N',NBAS(ISYMAL),NRHFFR(ISYMK),NBAS(ISYMBE), 2501 * ONE,WORK(KOFF3),NTOTAL,CMO(KOFF4),NTOTBE,ZERO, 2502 * WORK(KSCR1),NTOTAL) 2503C 2504 KOFF5 = ILRHSI(ISYMAL) + NBAS(ISYMAL)*NRHFFR(ISYMJ) + 1 2505 KOFF6 = KMOINT + ICOFRO(ISYMJ,ISYMK) 2506C 2507 NTOTAL = MAX(NBAS(ISYMAL),1) 2508 NTOTJ = MAX(NRHF(ISYMJ),1) 2509C 2510 CALL DGEMM('T','N',NRHF(ISYMJ),NRHFFR(ISYMK),NBAS(ISYMAL), 2511 * ONE,CMO(KOFF5),NTOTAL,WORK(KSCR1),NTOTAL,ZERO, 2512 * WORK(KOFF6),NTOTJ) 2513C 2514 110 CONTINUE 2515C 2516C----------------------------------- 2517C Calculate the contribution. 2518C----------------------------------- 2519C 2520 FACT = EIGHT*DDOT(NCOFRO(ISYMJK),WORK(KMOINT),1,ZKJK(1),1) 2521 KOFF7 = IT1FRO(ISYMA,ISYMI) + NVIR(ISYMA)*(I - 1) + 1 2522 CALL DAXPY(NVIR(ISYMA),FACT,WORK(KAVEC),1,ETAFRO(KOFF7),1) 2523C 2524 100 CONTINUE 2525C 2526 RETURN 2527 END 2528C /* Deck mp2_eidf3 */ 2529 SUBROUTINE MP2_EIDF3(ETAAI,XINOFO,ZKJK,CMO,WORK,LWORK,IDEL, 2530 * ISYDEL) 2531C 2532C Written by Asger Halkier 30/5 - 1998. 2533C 2534C To calculate the first exchange part of the indirect frozen 2535C contribution to the correlated part of eta in MP2 calculations. 2536C 2537#include "implicit.h" 2538 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 2539 DIMENSION ETAAI(*), XINOFO(*), ZKJK(*), CMO(*), WORK(LWORK) 2540#include "priunit.h" 2541#include "maxorb.h" 2542#include "ccsdinp.h" 2543#include "ccorb.h" 2544#include "ccsdsym.h" 2545#include "ccfro.h" 2546C 2547C-------------------------------------------------- 2548C Initial symmetries and work space allocation. 2549C-------------------------------------------------- 2550C 2551 ISYMA = ISYDEL 2552 ISYMI = ISYMA 2553 ISYMJK = 1 2554C 2555 IF (NRHF(ISYMI) .EQ. 0) RETURN 2556C 2557 KAVEC = 1 2558 KIVEC = KAVEC + NVIR(ISYMA) 2559 KEND1 = KIVEC + NRHF(ISYMI) 2560 LWRK1 = LWORK - KEND1 2561C 2562 IF (LWRK1 .LT. 0) THEN 2563 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 2564 CALL QUIT('Insufficient work space for allocation '// 2565 & 'in MP2_EIDF3') 2566 ENDIF 2567C 2568 CALL DZERO(WORK(KAVEC),NVIR(ISYMA)) 2569 CALL DZERO(WORK(KIVEC),NRHF(ISYMI)) 2570C 2571C----------------------------------- 2572C Copy vector out of CMO matrix. 2573C----------------------------------- 2574C 2575 KOFF1 = ILVISI(ISYDEL) + IDEL - IBAS(ISYDEL) 2576C 2577 CALL DCOPY(NVIR(ISYMA),CMO(KOFF1),NBAS(ISYDEL),WORK(KAVEC),1) 2578C 2579C---------------------------------------------- 2580C Loops over third index in integral batch. 2581C---------------------------------------------- 2582C 2583 DO 100 ISYMJ = 1,NSYM 2584C 2585 ISYMK = MULD2H(ISYMJK,ISYMJ) 2586 ISYMIK = MULD2H(ISYMI,ISYMK) 2587C 2588 DO 110 J = 1,NRHF(ISYMJ) 2589C 2590 KOFF2 = IOFROO(ISYMIK,ISYMJ) + NCOFRO(ISYMIK)*(J - 1) 2591 * + ICOFRO(ISYMI,ISYMK) + 1 2592 KOFF3 = ICOFRO(ISYMJ,ISYMK) + J 2593C 2594 NTOTI = MAX(NRHF(ISYMI),1) 2595C 2596 CALL DGEMV('N',NRHF(ISYMI),NRHFFR(ISYMK),ONE, 2597 * XINOFO(KOFF2),NTOTI,ZKJK(KOFF3),NRHF(ISYMJ), 2598 * ONE,WORK(KIVEC),1) 2599C 2600 110 CONTINUE 2601 100 CONTINUE 2602C 2603C----------------------------- 2604C Final storage in result. 2605C----------------------------- 2606C 2607 DO 120 I = 1,NRHF(ISYMI) 2608C 2609 KOFF4 = KIVEC + I - 1 2610 KOFF5 = IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I - 1) + 1 2611C 2612 CALL DAXPY(NVIR(ISYMA),-TWO*WORK(KOFF4),WORK(KAVEC),1, 2613 * ETAAI(KOFF5),1) 2614C 2615 120 CONTINUE 2616C 2617 RETURN 2618 END 2619C /* Deck mp2_eidf4 */ 2620 SUBROUTINE MP2_EIDF4(ETAAI,DSFRO,ZKJK,CMO,WORK,LWORK,IDEL, 2621 * ISYDEL) 2622C 2623C Written by Asger Halkier 30/5 - 1998. 2624C 2625C To calculate the second exchange part of the indirect frozen 2626C contribution to the correlated part of eta in MP2 calculations. 2627C 2628#include "implicit.h" 2629 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 2630 DIMENSION ETAAI(*), DSFRO(*), ZKJK(*), CMO(*), WORK(LWORK) 2631#include "priunit.h" 2632#include "maxorb.h" 2633#include "ccsdinp.h" 2634#include "ccorb.h" 2635#include "ccsdsym.h" 2636#include "ccfro.h" 2637C 2638C-------------------------------------------------- 2639C Initial symmetries and work space allocation. 2640C-------------------------------------------------- 2641C 2642 ISYMA = ISYDEL 2643 ISYMI = ISYMA 2644 ISYMJK = 1 2645C 2646 IF (NRHF(ISYMI) .EQ. 0) RETURN 2647C 2648 KAVEC = 1 2649 KIVEC = KAVEC + NVIR(ISYMA) 2650 KEND1 = KIVEC + NRHF(ISYMI) 2651 LWRK1 = LWORK - KEND1 2652C 2653 IF (LWRK1 .LT. 0) THEN 2654 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 2655 CALL QUIT('Insufficient work space for allocation '// 2656 & 'in MP2_EIDF4') 2657 ENDIF 2658C 2659 CALL DZERO(WORK(KAVEC),NVIR(ISYMA)) 2660 CALL DZERO(WORK(KIVEC),NRHF(ISYMI)) 2661C 2662C----------------------------------- 2663C Copy vector out of CMO matrix. 2664C----------------------------------- 2665C 2666 KOFF1 = ILVISI(ISYDEL) + IDEL - IBAS(ISYDEL) 2667C 2668 CALL DCOPY(NVIR(ISYMA),CMO(KOFF1),NBAS(ISYDEL),WORK(KAVEC),1) 2669C 2670C--------------------------------------------------------------------------- 2671C Loops over third index in integral batch & final wok space allocation. 2672C--------------------------------------------------------------------------- 2673C 2674 DO 100 ISYMK = 1,NSYM 2675C 2676 ISYMJ = MULD2H(ISYMJK,ISYMK) 2677 ISYMIJ = MULD2H(ISYMI,ISYMJ) 2678 ISALBE = ISYMIJ 2679 ISYMAL = ISYMI 2680 ISYMBE = ISYMJ 2681C 2682 KAOINT = KEND1 2683 KSCR1 = KAOINT + N2BST(ISALBE) 2684 KMOINT = KSCR1 + NBAS(ISYMAL)*NRHF(ISYMJ) 2685 KEND2 = KMOINT + NRHF(ISYMI)*NRHF(ISYMJ) 2686 LWRK2 = LWORK - KEND2 2687C 2688 IF (LWRK2 .LT. 0) THEN 2689 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND2 2690 CALL QUIT('Insufficient work space for allocation '// 2691 & 'in MP2_EIDF4') 2692 ENDIF 2693C 2694 CALL DZERO(WORK(KAOINT),N2BST(ISALBE)) 2695C 2696 DO 110 K = 1,NRHFFR(ISYMK) 2697C 2698C----------------------------------------------------------- 2699C Square up integral distribution (al be | K del). 2700C----------------------------------------------------------- 2701C 2702 KOFF2 = IDSFRO(ISALBE,ISYMK) + NNBST(ISALBE)*(K - 1) + 1 2703C 2704 CALL CCSD_SYMSQ(DSFRO(KOFF2),ISALBE,WORK(KAOINT)) 2705C 2706C--------------------------------------------- 2707C Calculate integrals (i j | K del). 2708C--------------------------------------------- 2709C 2710 KOFF3 = KAOINT + IAODIS(ISYMAL,ISYMBE) 2711 KOFF4 = ILRHSI(ISYMBE) + NBAS(ISYMBE)*NRHFFR(ISYMJ) + 1 2712C 2713 NTOTAL = MAX(NBAS(ISYMAL),1) 2714 NTOTBE = MAX(NBAS(ISYMBE),1) 2715C 2716 CALL DGEMM('N','N',NBAS(ISYMAL),NRHF(ISYMJ),NBAS(ISYMBE), 2717 * ONE,WORK(KOFF3),NTOTAL,CMO(KOFF4),NTOTBE,ZERO, 2718 * WORK(KSCR1),NTOTAL) 2719C 2720 KOFF5 = ILRHSI(ISYMAL) + NBAS(ISYMAL)*NRHFFR(ISYMI) + 1 2721C 2722 NTOTAL = MAX(NBAS(ISYMAL),1) 2723 NTOTI = MAX(NRHF(ISYMI),1) 2724C 2725 CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),NBAS(ISYMAL), 2726 * ONE,CMO(KOFF5),NTOTAL,WORK(KSCR1),NTOTAL,ZERO, 2727 * WORK(KMOINT),NTOTI) 2728C 2729C------------------------------------------- 2730C Contract MO integrals with ZKJK. 2731C------------------------------------------- 2732C 2733 KOFF6 = ICOFRO(ISYMJ,ISYMK) + NRHF(ISYMJ)*(K - 1) + 1 2734C 2735 NTOTI = MAX(NRHF(ISYMI),1) 2736C 2737 CALL DGEMV('N',NRHF(ISYMI),NRHF(ISYMJ),ONE,WORK(KMOINT), 2738 * NTOTI,ZKJK(KOFF6),1,ONE,WORK(KIVEC),1) 2739C 2740 110 CONTINUE 2741 100 CONTINUE 2742C 2743C----------------------------- 2744C Final storage in result. 2745C----------------------------- 2746C 2747 DO 120 I = 1,NRHF(ISYMI) 2748C 2749 KOFF7 = KIVEC + I - 1 2750 KOFF8 = IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I - 1) + 1 2751C 2752 CALL DAXPY(NVIR(ISYMA),-TWO*WORK(KOFF7),WORK(KAVEC),1, 2753 * ETAAI(KOFF8),1) 2754C 2755 120 CONTINUE 2756C 2757 RETURN 2758 END 2759C /* Deck mp2_eidf6 */ 2760 SUBROUTINE MP2_EIDF6(ETAFRO,DSFRO,ZKJK,CMO,WORK,LWORK,IDEL, 2761 * ISYDEL) 2762C 2763C Written by Asger Halkier 31/5 - 1998. 2764C 2765C To calculate the second exchange part of the indirect frozen 2766C contribution to the frozen part of eta in MP2 calculations. 2767C 2768#include "implicit.h" 2769 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 2770 DIMENSION ETAFRO(*), DSFRO(*), ZKJK(*), CMO(*), WORK(LWORK) 2771#include "priunit.h" 2772#include "maxorb.h" 2773#include "ccsdinp.h" 2774#include "ccorb.h" 2775#include "ccsdsym.h" 2776#include "ccfro.h" 2777C 2778C-------------------------------------------------- 2779C Initial symmetries and work space allocation. 2780C-------------------------------------------------- 2781C 2782 ISYMA = ISYDEL 2783 ISYMI = ISYMA 2784 ISYMJK = 1 2785C 2786 IF (NRHFFR(ISYMI) .EQ. 0) RETURN 2787C 2788 KAVEC = 1 2789 KIVEC = KAVEC + NVIR(ISYMA) 2790 KEND1 = KIVEC + NRHFFR(ISYMI) 2791 LWRK1 = LWORK - KEND1 2792C 2793 IF (LWRK1 .LT. 0) THEN 2794 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 2795 CALL QUIT('Insufficient work space for allocation '// 2796 & 'in MP2_EIDF6') 2797 ENDIF 2798C 2799 CALL DZERO(WORK(KAVEC),NVIR(ISYMA)) 2800 CALL DZERO(WORK(KIVEC),NRHFFR(ISYMI)) 2801C 2802C----------------------------------- 2803C Copy vector out of CMO matrix. 2804C----------------------------------- 2805C 2806 KOFF1 = ILVISI(ISYDEL) + IDEL - IBAS(ISYDEL) 2807C 2808 CALL DCOPY(NVIR(ISYMA),CMO(KOFF1),NBAS(ISYDEL),WORK(KAVEC),1) 2809C 2810C--------------------------------------------------------------------------- 2811C Loops over third index in integral batch & final wok space allocation. 2812C--------------------------------------------------------------------------- 2813C 2814 DO 100 ISYMK = 1,NSYM 2815C 2816 ISYMJ = MULD2H(ISYMJK,ISYMK) 2817 ISYMIJ = MULD2H(ISYMI,ISYMJ) 2818 ISALBE = ISYMIJ 2819 ISYMAL = ISYMI 2820 ISYMBE = ISYMJ 2821C 2822 KAOINT = KEND1 2823 KSCR1 = KAOINT + N2BST(ISALBE) 2824 KMOINT = KSCR1 + NBAS(ISYMAL)*NRHF(ISYMJ) 2825 KEND2 = KMOINT + NRHFFR(ISYMI)*NRHF(ISYMJ) 2826 LWRK2 = LWORK - KEND2 2827C 2828 IF (LWRK2 .LT. 0) THEN 2829 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND2 2830 CALL QUIT('Insufficient work space for allocation '// 2831 & 'in MP2_EIDF6') 2832 ENDIF 2833C 2834 CALL DZERO(WORK(KAOINT),N2BST(ISALBE)) 2835C 2836 DO 110 K = 1,NRHFFR(ISYMK) 2837C 2838C----------------------------------------------------------- 2839C Square up integral distribution (al be | K del). 2840C----------------------------------------------------------- 2841C 2842 KOFF2 = IDSFRO(ISALBE,ISYMK) + NNBST(ISALBE)*(K - 1) + 1 2843C 2844 CALL CCSD_SYMSQ(DSFRO(KOFF2),ISALBE,WORK(KAOINT)) 2845C 2846C--------------------------------------------- 2847C Calculate integrals (I j | K del). 2848C--------------------------------------------- 2849C 2850 KOFF3 = KAOINT + IAODIS(ISYMAL,ISYMBE) 2851 KOFF4 = ILRHSI(ISYMBE) + NBAS(ISYMBE)*NRHFFR(ISYMJ) + 1 2852C 2853 NTOTAL = MAX(NBAS(ISYMAL),1) 2854 NTOTBE = MAX(NBAS(ISYMBE),1) 2855C 2856 CALL DGEMM('N','N',NBAS(ISYMAL),NRHF(ISYMJ),NBAS(ISYMBE), 2857 * ONE,WORK(KOFF3),NTOTAL,CMO(KOFF4),NTOTBE,ZERO, 2858 * WORK(KSCR1),NTOTAL) 2859C 2860 KOFF5 = ILRHSI(ISYMAL) + 1 2861C 2862 NTOTAL = MAX(NBAS(ISYMAL),1) 2863 NTOTI = MAX(NRHFFR(ISYMI),1) 2864C 2865 CALL DGEMM('T','N',NRHFFR(ISYMI),NRHF(ISYMJ),NBAS(ISYMAL), 2866 * ONE,CMO(KOFF5),NTOTAL,WORK(KSCR1),NTOTAL,ZERO, 2867 * WORK(KMOINT),NTOTI) 2868C 2869C------------------------------------------- 2870C Contract MO integrals with ZKJK. 2871C------------------------------------------- 2872C 2873 KOFF6 = ICOFRO(ISYMJ,ISYMK) + NRHF(ISYMJ)*(K - 1) + 1 2874C 2875 NTOTI = MAX(NRHFFR(ISYMI),1) 2876C 2877 CALL DGEMV('N',NRHFFR(ISYMI),NRHF(ISYMJ),ONE,WORK(KMOINT), 2878 * NTOTI,ZKJK(KOFF6),1,ONE,WORK(KIVEC),1) 2879C 2880 110 CONTINUE 2881 100 CONTINUE 2882C 2883C----------------------------- 2884C Final storage in result. 2885C----------------------------- 2886C 2887 DO 120 I = 1,NRHFFR(ISYMI) 2888C 2889 KOFF7 = KIVEC + I - 1 2890 KOFF8 = IT1FRO(ISYMA,ISYMI) + NVIR(ISYMA)*(I - 1) + 1 2891C 2892 CALL DAXPY(NVIR(ISYMA),-TWO*WORK(KOFF7),WORK(KAVEC),1, 2893 * ETAFRO(KOFF8),1) 2894C 2895 120 CONTINUE 2896C 2897 RETURN 2898 END 2899C /* Deck mp2_eidf5 */ 2900 SUBROUTINE MP2_EIDF5(ETAFRO,DSRHF,ZKJK,CMO,WORK,LWORK,IDEL, 2901 * ISYDEL) 2902C 2903C Written by Asger Halkier 31/5 - 1998. 2904C 2905C To calculate the first exchange part of the indirect frozen 2906C contribution to the frozen part of eta in MP2 calculations. 2907C 2908#include "implicit.h" 2909 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 2910 DIMENSION ETAFRO(*), DSRHF(*), ZKJK(*), CMO(*), WORK(LWORK) 2911#include "priunit.h" 2912#include "maxorb.h" 2913#include "ccsdinp.h" 2914#include "ccorb.h" 2915#include "ccsdsym.h" 2916#include "ccfro.h" 2917C 2918C-------------------------------------------------- 2919C Initial symmetries and work space allocation. 2920C-------------------------------------------------- 2921C 2922 ISYMA = ISYDEL 2923 ISYMI = ISYMA 2924 ISYMJK = 1 2925C 2926 IF (NRHFFR(ISYMI) .EQ. 0) RETURN 2927C 2928 KAVEC = 1 2929 KIVEC = KAVEC + NVIR(ISYMA) 2930 KEND1 = KIVEC + NRHFFR(ISYMI) 2931 LWRK1 = LWORK - KEND1 2932C 2933 IF (LWRK1 .LT. 0) THEN 2934 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 2935 CALL QUIT('Insufficient work space for allocation '// 2936 & 'in MP2_EIDF5') 2937 ENDIF 2938C 2939 CALL DZERO(WORK(KAVEC),NVIR(ISYMA)) 2940 CALL DZERO(WORK(KIVEC),NRHFFR(ISYMI)) 2941C 2942C----------------------------------- 2943C Copy vector out of CMO matrix. 2944C----------------------------------- 2945C 2946 KOFF1 = ILVISI(ISYDEL) + IDEL - IBAS(ISYDEL) 2947C 2948 CALL DCOPY(NVIR(ISYMA),CMO(KOFF1),NBAS(ISYDEL),WORK(KAVEC),1) 2949C 2950C--------------------------------------------------------------------------- 2951C Loops over third index in integral batch & final wok space allocation. 2952C--------------------------------------------------------------------------- 2953C 2954 DO 100 ISYMJ = 1,NSYM 2955C 2956 ISYMK = MULD2H(ISYMJK,ISYMJ) 2957 ISYMIK = MULD2H(ISYMI,ISYMK) 2958 ISALBE = ISYMIK 2959 ISYMAL = ISYMI 2960 ISYMBE = ISYMK 2961C 2962 IF (NRHFFR(ISYMK) .EQ. 0) GOTO 100 2963C 2964 KAOINT = KEND1 2965 KSCR1 = KAOINT + N2BST(ISALBE) 2966 KMOINT = KSCR1 + NBAS(ISYMAL)*NRHFFR(ISYMK) 2967 KEND2 = KMOINT + NRHFFR(ISYMI)*NRHFFR(ISYMK) 2968 LWRK2 = LWORK - KEND2 2969C 2970 IF (LWRK2 .LT. 0) THEN 2971 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND2 2972 CALL QUIT('Insufficient work space for allocation '// 2973 & 'in MP2_EIDF5') 2974 ENDIF 2975C 2976 CALL DZERO(WORK(KAOINT),KEND2-KEND1) 2977C 2978 DO 110 J = 1,NRHF(ISYMJ) 2979C 2980C----------------------------------------------------------- 2981C Square up integral distribution (al be | j del). 2982C----------------------------------------------------------- 2983C 2984 KOFF2 = IDSRHF(ISALBE,ISYMJ) + NNBST(ISALBE)*(J - 1) + 1 2985C 2986 CALL CCSD_SYMSQ(DSRHF(KOFF2),ISALBE,WORK(KAOINT)) 2987C 2988C--------------------------------------------- 2989C Calculate integrals (I K | j del). 2990C--------------------------------------------- 2991C 2992 KOFF3 = KAOINT + IAODIS(ISYMAL,ISYMBE) 2993 KOFF4 = ILRHSI(ISYMBE) + 1 2994C 2995 NTOTAL = MAX(NBAS(ISYMAL),1) 2996 NTOTBE = MAX(NBAS(ISYMBE),1) 2997C 2998 CALL DGEMM('N','N',NBAS(ISYMAL),NRHFFR(ISYMK),NBAS(ISYMBE), 2999 * ONE,WORK(KOFF3),NTOTAL,CMO(KOFF4),NTOTBE,ZERO, 3000 * WORK(KSCR1),NTOTAL) 3001C 3002 KOFF5 = ILRHSI(ISYMAL) + 1 3003C 3004 NTOTAL = MAX(NBAS(ISYMAL),1) 3005 NTOTI = MAX(NRHFFR(ISYMI),1) 3006C 3007 CALL DGEMM('T','N',NRHFFR(ISYMI),NRHFFR(ISYMK), 3008 * NBAS(ISYMAL),ONE,CMO(KOFF5),NTOTAL,WORK(KSCR1), 3009 * NTOTAL,ZERO,WORK(KMOINT),NTOTI) 3010C 3011C------------------------------------------- 3012C Contract MO integrals with ZKJK. 3013C------------------------------------------- 3014C 3015 KOFF6 = ICOFRO(ISYMJ,ISYMK) + J 3016C 3017 NTOTI = MAX(NRHFFR(ISYMI),1) 3018C 3019 CALL DGEMV('N',NRHFFR(ISYMI),NRHFFR(ISYMK),ONE, 3020 * WORK(KMOINT),NTOTI,ZKJK(KOFF6),NRHF(ISYMJ), 3021 * ONE,WORK(KIVEC),1) 3022C 3023 110 CONTINUE 3024 100 CONTINUE 3025C 3026C----------------------------- 3027C Final storage in result. 3028C----------------------------- 3029C 3030 DO 120 I = 1,NRHFFR(ISYMI) 3031C 3032 KOFF7 = KIVEC + I - 1 3033 KOFF8 = IT1FRO(ISYMA,ISYMI) + NVIR(ISYMA)*(I - 1) + 1 3034C 3035 CALL DAXPY(NVIR(ISYMA),-TWO*WORK(KOFF7),WORK(KAVEC),1, 3036 * ETAFRO(KOFF8),1) 3037C 3038 120 CONTINUE 3039C 3040 RETURN 3041 END 3042C /* Deck ccifromo */ 3043 SUBROUTINE CCIFROMO(XFIJ,XFJI,XFAI,XFIA,XFII,XAOINT, 3044 * XLAMDP,XLAMDH,CMO,WORK,LWORK,ISYM) 3045C 3046C Written by Asger Halkier 13/8 - 1998 3047C 3048C Version: 1.0 3049C 3050C Purpose: To transform the AO integrals in (XAOINT) to 3051C frozen MO-basis (stored in XINTIJ through XINTAI)! 3052C Note that the AO integrals either can be the one 3053C electron integrals or a submatrix alfa,beta of the two 3054C electron integrals (alfa beta|gamma delta)! 3055C ISYM is the integral symmetry! 3056C 3057#include "implicit.h" 3058 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 3059 DIMENSION XFIJ(*), XFJI(*), XFAI(*), XFIA(*), XFII(*) 3060 DIMENSION XAOINT(*), XLAMDP(*), XLAMDH(*), CMO(*), WORK(LWORK) 3061#include "priunit.h" 3062#include "maxorb.h" 3063#include "ccsdinp.h" 3064#include "ccorb.h" 3065#include "ccsdsym.h" 3066#include "ccfro.h" 3067C 3068C------------------------------ 3069C Initialize result arrays. 3070C------------------------------ 3071C 3072 CALL DZERO(XFIJ,NCOFRO(ISYM)) 3073 CALL DZERO(XFJI,NCOFRO(ISYM)) 3074 CALL DZERO(XFAI,NT1FRO(ISYM)) 3075 CALL DZERO(XFIA,NT1FRO(ISYM)) 3076 CALL DZERO(XFII,NFROFR(ISYM)) 3077C 3078 DO 100 ISYMAL = 1,NSYM 3079C 3080 ISYMBE = MULD2H(ISYMAL,ISYM) 3081 ISYMJ = ISYMBE 3082 ISYMI = ISYMAL 3083C 3084C------------------------------ 3085C Work space allocation. 3086C------------------------------ 3087C 3088 LSCR = MAX(NBAS(ISYMAL)*NVIR(ISYMBE), 3089 * NBAS(ISYMAL)*NRHFFR(ISYMBE)) 3090C 3091 KSCR1 = 1 3092 KEND1 = KSCR1 + LSCR 3093 LWRK1 = LWORK - KEND1 3094C 3095 IF (LWRK1 .LT. 0) THEN 3096 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 3097 CALL QUIT('Insufficient memory in CCIFROMO') 3098 ENDIF 3099C 3100 CALL DZERO(WORK(KSCR1),LSCR) 3101C 3102C-------------------------------------------------- 3103C Transform second integral index to frozen. 3104C-------------------------------------------------- 3105C 3106 KOFF1 = IAODIS(ISYMAL,ISYMBE) + 1 3107 KOFF2 = ILRHSI(ISYMBE) + 1 3108C 3109 NTOTAL = MAX(NBAS(ISYMAL),1) 3110 NTOTBE = MAX(NBAS(ISYMBE),1) 3111C 3112 CALL DGEMM('N','N',NBAS(ISYMAL),NRHFFR(ISYMJ),NBAS(ISYMBE), 3113 * ONE,XAOINT(KOFF1),NTOTAL,CMO(KOFF2),NTOTBE,ZERO, 3114 * WORK(KSCR1),NTOTAL) 3115C 3116C---------------------------------------------- 3117C Calculate the frozen-frozen integrals. 3118C---------------------------------------------- 3119C 3120 KOFF3 = ILRHSI(ISYMAL) + 1 3121 KOFF4 = IFROFR(ISYMI,ISYMJ) + 1 3122C 3123 NTOTAL = MAX(NBAS(ISYMAL),1) 3124 NTOTI = MAX(NRHFFR(ISYMI),1) 3125C 3126 CALL DGEMM('T','N',NRHFFR(ISYMI),NRHFFR(ISYMJ),NBAS(ISYMAL), 3127 * ONE,CMO(KOFF3),NTOTAL,WORK(KSCR1),NTOTAL,ZERO, 3128 * XFII(KOFF4),NTOTI) 3129C 3130C-------------------------------------------------- 3131C Calculate the correlated-frozen integrals. 3132C-------------------------------------------------- 3133C 3134 KOFF5 = ILMRHF(ISYMAL) + 1 3135 KOFF6 = ICOFRO(ISYMI,ISYMJ) + 1 3136C 3137 NTOTAL = MAX(NBAS(ISYMAL),1) 3138 NTOTI = MAX(NRHF(ISYMI),1) 3139C 3140 CALL DGEMM('T','N',NRHF(ISYMI),NRHFFR(ISYMJ),NBAS(ISYMAL), 3141 * ONE,XLAMDP(KOFF5),NTOTAL,WORK(KSCR1),NTOTAL, 3142 * ZERO,XFIJ(KOFF6),NTOTI) 3143C 3144C----------------------------------------------- 3145C Calculate the virtual-frozen integrals. 3146C----------------------------------------------- 3147C 3148 ISYMA = ISYMAL 3149C 3150 KOFF7 = ILMVIR(ISYMAL) + 1 3151 KOFF8 = IT1FRO(ISYMA,ISYMJ) + 1 3152C 3153 NTOTAL = MAX(NBAS(ISYMAL),1) 3154 NTOTA = MAX(NVIR(ISYMA),1) 3155C 3156 CALL DGEMM('T','N',NVIR(ISYMA),NRHFFR(ISYMJ),NBAS(ISYMAL), 3157 * ONE,XLAMDP(KOFF7),NTOTAL,WORK(KSCR1),NTOTAL, 3158 * ZERO,XFAI(KOFF8),NTOTA) 3159C 3160C------------------------------------------------------ 3161C Transform second integral index to correlated. 3162C------------------------------------------------------ 3163C 3164 CALL DZERO(WORK(KSCR1),LSCR) 3165C 3166 KOFF9 = IAODIS(ISYMAL,ISYMBE) + 1 3167 KOFF10 = ILMRHF(ISYMBE) + 1 3168C 3169 NTOTAL = MAX(NBAS(ISYMAL),1) 3170 NTOTBE = MAX(NBAS(ISYMBE),1) 3171C 3172 CALL DGEMM('N','N',NBAS(ISYMAL),NRHF(ISYMJ),NBAS(ISYMBE), 3173 * ONE,XAOINT(KOFF9),NTOTAL,XLAMDH(KOFF10),NTOTBE, 3174 * ZERO,WORK(KSCR1),NTOTAL) 3175C 3176C-------------------------------------------------- 3177C Calculate the frozen-correlated integrals. 3178C-------------------------------------------------- 3179C 3180 KOFF11 = ILRHSI(ISYMAL) + 1 3181 KOFF12 = ICOFRO(ISYMJ,ISYMI) + 1 3182C 3183 NTOTAL = MAX(NBAS(ISYMAL),1) 3184 NTOTJ = MAX(NRHF(ISYMJ),1) 3185C 3186 CALL DGEMM('T','N',NRHF(ISYMJ),NRHFFR(ISYMI),NBAS(ISYMAL), 3187 * ONE,WORK(KSCR1),NTOTAL,CMO(KOFF11),NTOTAL,ZERO, 3188 * XFJI(KOFF12),NTOTJ) 3189C 3190C--------------------------------------------------- 3191C Transform second integral index to virtual. 3192C--------------------------------------------------- 3193C 3194 CALL DZERO(WORK(KSCR1),LSCR) 3195C 3196 ISYMA = ISYMBE 3197C 3198 KOFF13 = IAODIS(ISYMAL,ISYMBE) + 1 3199 KOFF14 = ILMVIR(ISYMBE) + 1 3200C 3201 NTOTAL = MAX(NBAS(ISYMAL),1) 3202 NTOTBE = MAX(NBAS(ISYMBE),1) 3203C 3204 CALL DGEMM('N','N',NBAS(ISYMAL),NVIR(ISYMA),NBAS(ISYMBE), 3205 * ONE,XAOINT(KOFF13),NTOTAL,XLAMDH(KOFF14),NTOTBE, 3206 * ZERO,WORK(KSCR1),NTOTAL) 3207C 3208C----------------------------------------------- 3209C Calculate the frozen-virtual integrals. 3210C----------------------------------------------- 3211C 3212 KOFF13 = ILRHSI(ISYMAL) + 1 3213 KOFF14 = IT1FRO(ISYMA,ISYMI) + 1 3214C 3215 NTOTAL = MAX(NBAS(ISYMAL),1) 3216 NTOTA = MAX(NVIR(ISYMA),1) 3217C 3218 CALL DGEMM('T','N',NVIR(ISYMA),NRHFFR(ISYMI),NBAS(ISYMAL), 3219 * ONE,WORK(KSCR1),NTOTAL,CMO(KOFF13),NTOTAL,ZERO, 3220 * XFIA(KOFF14),NTOTA) 3221C 3222 100 CONTINUE 3223C 3224 RETURN 3225 END 3226C /* Deck ccfretai */ 3227 SUBROUTINE CCFRETAI(ETAAI,XFIJ,XFJI,XFAI,XFIA,DFIJ,DFJI, 3228 * DFAI,DFIA,T1AM,WORK,LWORK,ISYM) 3229C 3230C Written by Asger Halkier 13/8 - 1998 3231C 3232C Version: 1.0 3233C 3234C Purpose: To calculate the contributions to eta-ai-0 from 3235C intermediate looping over frozen. 3236C 3237#include "implicit.h" 3238 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 3239 DIMENSION ETAAI(*), XFIJ(*), XFJI(*), XFAI(*), XFIA(*), DFIJ(*) 3240 DIMENSION DFJI(*), DFAI(*), DFIA(*), T1AM(*), WORK(LWORK) 3241#include "priunit.h" 3242#include "maxorb.h" 3243#include "ccsdinp.h" 3244#include "ccorb.h" 3245#include "ccsdsym.h" 3246#include "ccfro.h" 3247C 3248C------------------------------ 3249C Do the four direct terms. 3250C------------------------------ 3251C 3252 DO 100 ISYMA = 1,NSYM 3253C 3254 ISYMI = ISYMA 3255 ISYMJ = MULD2H(ISYMA,ISYM) 3256C 3257 IF (NRHFFR(ISYMJ) .EQ. 0) GOTO 100 3258C 3259 KOFFR = IT1AM(ISYMA,ISYMI) + 1 3260 KOFF1 = IT1FRO(ISYMA,ISYMJ) + 1 3261 KOFF2 = ICOFRO(ISYMI,ISYMJ) + 1 3262C 3263 NTOTA = MAX(NVIR(ISYMA),1) 3264 NTOTI = MAX(NRHF(ISYMI),1) 3265C 3266 CALL DGEMM('N','T',NVIR(ISYMA),NRHF(ISYMI),NRHFFR(ISYMJ), 3267 * ONE,XFIA(KOFF1),NTOTA,DFJI(KOFF2),NTOTI,ONE, 3268 * ETAAI(KOFFR),NTOTA) 3269C 3270 CALL DGEMM('N','T',NVIR(ISYMA),NRHF(ISYMI),NRHFFR(ISYMJ), 3271 * -ONE,DFAI(KOFF1),NTOTA,XFIJ(KOFF2),NTOTI,ONE, 3272 * ETAAI(KOFFR),NTOTA) 3273C 3274 CALL DGEMM('N','T',NVIR(ISYMA),NRHF(ISYMI),NRHFFR(ISYMJ), 3275 * -ONE,DFIA(KOFF1),NTOTA,XFJI(KOFF2),NTOTI,ONE, 3276 * ETAAI(KOFFR),NTOTA) 3277C 3278 CALL DGEMM('N','T',NVIR(ISYMA),NRHF(ISYMI),NRHFFR(ISYMJ), 3279 * ONE,XFAI(KOFF1),NTOTA,DFIJ(KOFF2),NTOTI,ONE, 3280 * ETAAI(KOFFR),NTOTA) 3281C 3282 100 CONTINUE 3283C 3284C----------------------------------------------------------- 3285C Do the terms involving one T1AM and loop over virtual. 3286C----------------------------------------------------------- 3287C 3288 DO 110 ISYMA = 1,NSYM 3289C 3290 ISYMI = ISYMA 3291 ISYMC = ISYMI 3292 ISYMJ = MULD2H(ISYMA,ISYM) 3293C 3294 LSCR = NVIR(ISYMA)*NVIR(ISYMC) 3295C 3296 KSCR1 = 1 3297 KEND1 = KSCR1 + LSCR 3298 LWRK1 = LWORK - KEND1 3299C 3300 IF (LWRK1 .LT. 0) THEN 3301 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 3302 CALL QUIT('Insufficient memory in CCFRETAI') 3303 ENDIF 3304C 3305 CALL DZERO(WORK(KSCR1),LSCR) 3306C 3307 KOFF3 = IT1FRO(ISYMA,ISYMJ) + 1 3308 KOFF4 = IT1FRO(ISYMC,ISYMJ) + 1 3309 KOFFT = IT1AM(ISYMC,ISYMI) + 1 3310 KOFFR = IT1AM(ISYMA,ISYMI) + 1 3311C 3312 NTOTA = MAX(NVIR(ISYMA),1) 3313 NTOTC = MAX(NVIR(ISYMC),1) 3314C 3315 CALL DGEMM('N','T',NVIR(ISYMA),NVIR(ISYMC),NRHFFR(ISYMJ), 3316 * ONE,DFIA(KOFF3),NTOTA,XFIA(KOFF4),NTOTC,ZERO, 3317 * WORK(KSCR1),NTOTA) 3318C 3319 CALL DGEMM('N','T',NVIR(ISYMA),NVIR(ISYMC),NRHFFR(ISYMJ), 3320 * -ONE,XFAI(KOFF3),NTOTA,DFAI(KOFF4),NTOTC,ONE, 3321 * WORK(KSCR1),NTOTA) 3322C 3323 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI),NVIR(ISYMC), 3324 * ONE,WORK(KSCR1),NTOTA,T1AM(KOFFT),NTOTC,ONE, 3325 * ETAAI(KOFFR),NTOTA) 3326C 3327 110 CONTINUE 3328C 3329C-------------------------------------------------------------- 3330C Do the terms involving one T1AM and loop over correlated. 3331C-------------------------------------------------------------- 3332C 3333 DO 120 ISYMA = 1,NSYM 3334C 3335 ISYMI = ISYMA 3336 ISYMK = ISYMA 3337 ISYMJ = MULD2H(ISYMI,ISYM) 3338C 3339 LSCR = NRHF(ISYMK)*NRHF(ISYMI) 3340C 3341 KSCR2 = 1 3342 KEND2 = KSCR2 + LSCR 3343 LWRK2 = LWORK - KEND2 3344C 3345 IF (LWRK2 .LT. 0) THEN 3346 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND2 3347 CALL QUIT('Insufficient memory in CCFRETAI') 3348 ENDIF 3349C 3350 CALL DZERO(WORK(KSCR2),LSCR) 3351C 3352 KOFF5 = ICOFRO(ISYMK,ISYMJ) + 1 3353 KOFF6 = ICOFRO(ISYMI,ISYMJ) + 1 3354 KOFFT = IT1AM(ISYMA,ISYMK) + 1 3355 KOFFR = IT1AM(ISYMA,ISYMI) + 1 3356C 3357 NTOTK = MAX(NRHF(ISYMK),1) 3358 NTOTI = MAX(NRHF(ISYMI),1) 3359 NTOTA = MAX(NVIR(ISYMA),1) 3360C 3361 CALL DGEMM('N','T',NRHF(ISYMK),NRHF(ISYMI),NRHFFR(ISYMJ), 3362 * ONE,DFJI(KOFF5),NTOTK,XFJI(KOFF6),NTOTI,ZERO, 3363 * WORK(KSCR2),NTOTK) 3364C 3365 CALL DGEMM('N','T',NRHF(ISYMK),NRHF(ISYMI),NRHFFR(ISYMJ), 3366 * -ONE,XFIJ(KOFF5),NTOTK,DFIJ(KOFF6),NTOTI,ONE, 3367 * WORK(KSCR2),NTOTK) 3368C 3369 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI),NRHF(ISYMK), 3370 * -ONE,T1AM(KOFFT),NTOTA,WORK(KSCR2),NTOTK,ONE, 3371 * ETAAI(KOFFR),NTOTA) 3372C 3373 120 CONTINUE 3374C 3375C-------------------------------------- 3376C Do the terms involving two T1AMs. 3377C-------------------------------------- 3378C 3379 DO 130 ISYMA = 1,NSYM 3380C 3381 ISYMI = ISYMA 3382 ISYMK = ISYMA 3383 ISYMC = ISYMI 3384 ISYMJ = MULD2H(ISYMK,ISYM) 3385C 3386 LSCR3 = NRHF(ISYMK)*NVIR(ISYMC) 3387 LSCR4 = NRHF(ISYMK)*NRHF(ISYMI) 3388C 3389 KSCR3 = 1 3390 KSCR4 = KSCR3 + LSCR3 3391 KEND3 = KSCR4 + LSCR4 3392 LWRK3 = LWORK - KEND3 3393C 3394 IF (LWRK3 .LT. 0) THEN 3395 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND3 3396 CALL QUIT('Insufficient memory in CCFRETAI') 3397 ENDIF 3398C 3399 CALL DZERO(WORK(KSCR3),KEND3) 3400C 3401 KOFF7 = ICOFRO(ISYMK,ISYMJ) + 1 3402 KOFF8 = IT1FRO(ISYMC,ISYMJ) + 1 3403 KOFT1 = IT1AM(ISYMC,ISYMI) + 1 3404 KOFT2 = IT1AM(ISYMA,ISYMK) + 1 3405 KOFFR = IT1AM(ISYMA,ISYMI) + 1 3406C 3407 NTOTK = MAX(NRHF(ISYMK),1) 3408 NTOTC = MAX(NVIR(ISYMC),1) 3409 NTOTA = MAX(NVIR(ISYMA),1) 3410C 3411 CALL DGEMM('N','T',NRHF(ISYMK),NVIR(ISYMC),NRHFFR(ISYMJ), 3412 * ONE,DFJI(KOFF7),NTOTK,XFIA(KOFF8),NTOTC,ZERO, 3413 * WORK(KSCR3),NTOTK) 3414C 3415 CALL DGEMM('N','T',NRHF(ISYMK),NVIR(ISYMC),NRHFFR(ISYMJ), 3416 * -ONE,XFIJ(KOFF7),NTOTK,DFAI(KOFF8),NTOTC,ONE, 3417 * WORK(KSCR3),NTOTK) 3418C 3419 CALL DGEMM('N','N',NRHF(ISYMK),NRHF(ISYMI),NVIR(ISYMC), 3420 * ONE,WORK(KSCR3),NTOTK,T1AM(KOFT1),NTOTC,ZERO, 3421 * WORK(KSCR4),NTOTK) 3422C 3423 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI),NRHF(ISYMK), 3424 * ONE,T1AM(KOFT2),NTOTA,WORK(KSCR4),NTOTK,ONE, 3425 * ETAAI(KOFFR),NTOTA) 3426C 3427 130 CONTINUE 3428C 3429 RETURN 3430 END 3431C /* Deck ccfretab */ 3432 SUBROUTINE CCFRETAB(ETAAB,XFIJ,XFJI,XFAI,XFIA,DFIJ,DFJI, 3433 * DFAI,DFIA,T1AM,WORK,LWORK,ISYM) 3434C 3435C Written by Asger Halkier 14/8 - 1998 3436C 3437C Version: 1.0 3438C 3439C Purpose: To calculate the contributions to eta-ab-0 from 3440C intermediate looping over frozen. 3441C 3442#include "implicit.h" 3443 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 3444 DIMENSION ETAAB(*), XFIJ(*), XFJI(*), XFAI(*), XFIA(*), DFIJ(*) 3445 DIMENSION DFJI(*), DFAI(*), DFIA(*), T1AM(*), WORK(LWORK) 3446#include "priunit.h" 3447#include "maxorb.h" 3448#include "ccsdinp.h" 3449#include "ccorb.h" 3450#include "ccsdsym.h" 3451#include "ccfro.h" 3452C 3453C------------------------------ 3454C Do the four direct terms. 3455C------------------------------ 3456C 3457 DO 100 ISYMA = 1,NSYM 3458C 3459 ISYMB = ISYMA 3460 ISYMJ = MULD2H(ISYMA,ISYM) 3461C 3462 IF (NRHFFR(ISYMJ) .EQ. 0) GOTO 100 3463C 3464 KOFF1 = IT1FRO(ISYMA,ISYMJ) + 1 3465 KOFF2 = IT1FRO(ISYMB,ISYMJ) + 1 3466 KOFFR = NMATIJ(1) + IMATAB(ISYMA,ISYMB) + 1 3467C 3468 NTOTA = MAX(NVIR(ISYMA),1) 3469 NTOTB = MAX(NVIR(ISYMB),1) 3470C 3471 CALL DGEMM('N','T',NVIR(ISYMA),NVIR(ISYMB),NRHFFR(ISYMJ), 3472 * ONE,XFIA(KOFF1),NTOTA,DFIA(KOFF2),NTOTB,ONE, 3473 * ETAAB(KOFFR),NTOTA) 3474C 3475 CALL DGEMM('N','T',NVIR(ISYMA),NVIR(ISYMB),NRHFFR(ISYMJ), 3476 * -ONE,DFAI(KOFF1),NTOTA,XFAI(KOFF2),NTOTB,ONE, 3477 * ETAAB(KOFFR),NTOTA) 3478C 3479 CALL DGEMM('N','T',NVIR(ISYMA),NVIR(ISYMB),NRHFFR(ISYMJ), 3480 * -ONE,DFIA(KOFF1),NTOTA,XFIA(KOFF2),NTOTB,ONE, 3481 * ETAAB(KOFFR),NTOTA) 3482C 3483 CALL DGEMM('N','T',NVIR(ISYMA),NVIR(ISYMB),NRHFFR(ISYMJ), 3484 * ONE,XFAI(KOFF1),NTOTA,DFAI(KOFF2),NTOTB,ONE, 3485 * ETAAB(KOFFR),NTOTA) 3486C 3487 100 CONTINUE 3488C 3489C-------------------------------------- 3490C Do the two first terms with T1AM. 3491C-------------------------------------- 3492C 3493 DO 110 ISYMA = 1,NSYM 3494C 3495 ISYMB = ISYMA 3496 ISYMK = ISYMB 3497 ISYMJ = MULD2H(ISYMA,ISYM) 3498C 3499 LSCR = NVIR(ISYMA)*NRHF(ISYMK) 3500C 3501 KSCR1 = 1 3502 KEND1 = KSCR1 + LSCR 3503 LWRK1 = LWORK - KEND1 3504C 3505 IF (LWRK1 .LT. 0) THEN 3506 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 3507 CALL QUIT('Insufficient memory in CCFRETAB') 3508 ENDIF 3509C 3510 CALL DZERO(WORK(KSCR1),LSCR) 3511C 3512 KOFF3 = IT1FRO(ISYMA,ISYMJ) + 1 3513 KOFF4 = ICOFRO(ISYMK,ISYMJ) + 1 3514 KOFFT = IT1AM(ISYMB,ISYMK) + 1 3515 KOFFR = NMATIJ(1) + IMATAB(ISYMA,ISYMB) + 1 3516C 3517 NTOTA = MAX(NVIR(ISYMA),1) 3518 NTOTK = MAX(NRHF(ISYMK),1) 3519 NTOTB = MAX(NVIR(ISYMB),1) 3520C 3521 CALL DGEMM('N','T',NVIR(ISYMA),NRHF(ISYMK),NRHFFR(ISYMJ), 3522 * ONE,XFIA(KOFF3),NTOTA,DFJI(KOFF4),NTOTK,ZERO, 3523 * WORK(KSCR1),NTOTA) 3524C 3525 CALL DGEMM('N','T',NVIR(ISYMA),NRHF(ISYMK),NRHFFR(ISYMJ), 3526 * -ONE,DFAI(KOFF3),NTOTA,XFIJ(KOFF4),NTOTK,ONE, 3527 * WORK(KSCR1),NTOTA) 3528C 3529 CALL DGEMM('N','T',NVIR(ISYMA),NVIR(ISYMB),NRHF(ISYMK), 3530 * ONE,WORK(KSCR1),NTOTA,T1AM(KOFFT),NTOTB,ONE, 3531 * ETAAB(KOFFR),NTOTA) 3532C 3533 110 CONTINUE 3534C 3535C------------------------------------- 3536C Do the last two terms with T1AM. 3537C------------------------------------- 3538C 3539 DO 120 ISYMA = 1,NSYM 3540C 3541 ISYMB = ISYMA 3542 ISYMK = ISYMA 3543 ISYMJ = MULD2H(ISYMA,ISYM) 3544C 3545 LSCR = NVIR(ISYMB)*NRHF(ISYMK) 3546C 3547 KSCR2 = 1 3548 KEND2 = KSCR2 + LSCR 3549 LWRK2 = LWORK - KEND2 3550C 3551 IF (LWRK2 .LT. 0) THEN 3552 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND2 3553 CALL QUIT('Insufficient memory in CCFRETAB') 3554 ENDIF 3555C 3556 CALL DZERO(WORK(KSCR2),LSCR) 3557C 3558 KOFF5 = ICOFRO(ISYMK,ISYMJ) + 1 3559 KOFF6 = IT1FRO(ISYMB,ISYMJ) + 1 3560 KOFFT = IT1AM(ISYMA,ISYMK) + 1 3561 KOFFR = NMATIJ(1) + IMATAB(ISYMA,ISYMB) + 1 3562C 3563 NTOTK = MAX(NRHF(ISYMK),1) 3564 NTOTB = MAX(NVIR(ISYMB),1) 3565 NTOTA = MAX(NVIR(ISYMA),1) 3566C 3567 CALL DGEMM('N','T',NRHF(ISYMK),NVIR(ISYMB),NRHFFR(ISYMJ), 3568 * ONE,DFJI(KOFF5),NTOTK,XFIA(KOFF6),NTOTB,ZERO, 3569 * WORK(KSCR2),NTOTK) 3570C 3571 CALL DGEMM('N','T',NRHF(ISYMK),NVIR(ISYMB),NRHFFR(ISYMJ), 3572 * -ONE,XFIJ(KOFF5),NTOTK,DFAI(KOFF6),NTOTB,ONE, 3573 * WORK(KSCR2),NTOTK) 3574C 3575 CALL DGEMM('N','N',NVIR(ISYMA),NVIR(ISYMB),NRHF(ISYMK), 3576 * -ONE,T1AM(KOFFT),NTOTA,WORK(KSCR2),NTOTK,ONE, 3577 * ETAAB(KOFFR),NTOTA) 3578C 3579 120 CONTINUE 3580C 3581 RETURN 3582 END 3583C /* Deck ccfretij */ 3584 SUBROUTINE CCFRETIJ(ETAIJ,XFIJ,XFJI,XFAI,XFIA,DFIJ,DFJI, 3585 * DFAI,DFIA,T1AM,WORK,LWORK,ISYM) 3586C 3587C Written by Asger Halkier 14/8 - 1998 3588C 3589C Version: 1.0 3590C 3591C Purpose: To calculate the contributions to eta-ij-0 from 3592C intermediate looping over frozen. 3593C 3594#include "implicit.h" 3595 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 3596 DIMENSION ETAIJ(*), XFIJ(*), XFJI(*), XFAI(*), XFIA(*), DFIJ(*) 3597 DIMENSION DFJI(*), DFAI(*), DFIA(*), T1AM(*), WORK(LWORK) 3598#include "priunit.h" 3599#include "maxorb.h" 3600#include "ccsdinp.h" 3601#include "ccorb.h" 3602#include "ccsdsym.h" 3603#include "ccfro.h" 3604C 3605C------------------------------ 3606C Do the four direct terms. 3607C------------------------------ 3608C 3609 DO 100 ISYMI = 1,NSYM 3610C 3611 ISYMJ = ISYMI 3612 ISYMK = MULD2H(ISYMI,ISYM) 3613C 3614 IF (NRHFFR(ISYMK) .EQ. 0) GOTO 100 3615C 3616 KOFF1 = ICOFRO(ISYMI,ISYMK) + 1 3617 KOFF2 = ICOFRO(ISYMJ,ISYMK) + 1 3618 KOFFR = IMATIJ(ISYMI,ISYMJ) + 1 3619C 3620 NTOTI = MAX(NRHF(ISYMI),1) 3621 NTOTJ = MAX(NRHF(ISYMJ),1) 3622C 3623 CALL DGEMM('N','T',NRHF(ISYMI),NRHF(ISYMJ),NRHFFR(ISYMK), 3624 * ONE,XFJI(KOFF1),NTOTI,DFJI(KOFF2),NTOTJ,ONE, 3625 * ETAIJ(KOFFR),NTOTI) 3626C 3627 CALL DGEMM('N','T',NRHF(ISYMI),NRHF(ISYMJ),NRHFFR(ISYMK), 3628 * -ONE,DFIJ(KOFF1),NTOTI,XFIJ(KOFF2),NTOTJ,ONE, 3629 * ETAIJ(KOFFR),NTOTI) 3630C 3631 CALL DGEMM('N','T',NRHF(ISYMI),NRHF(ISYMJ),NRHFFR(ISYMK), 3632 * -ONE,DFJI(KOFF1),NTOTI,XFJI(KOFF2),NTOTJ,ONE, 3633 * ETAIJ(KOFFR),NTOTI) 3634C 3635 CALL DGEMM('N','T',NRHF(ISYMI),NRHF(ISYMJ),NRHFFR(ISYMK), 3636 * ONE,XFIJ(KOFF1),NTOTI,DFIJ(KOFF2),NTOTJ,ONE, 3637 * ETAIJ(KOFFR),NTOTI) 3638C 3639 100 CONTINUE 3640C 3641C-------------------------------------- 3642C Do the two first terms with T1AM. 3643C-------------------------------------- 3644C 3645 DO 110 ISYMI = 1,NSYM 3646C 3647 ISYMJ = ISYMI 3648 ISYMC = ISYMJ 3649 ISYMK = MULD2H(ISYMI,ISYM) 3650C 3651 LSCR = NRHF(ISYMI)*NVIR(ISYMC) 3652C 3653 KSCR1 = 1 3654 KEND1 = KSCR1 + LSCR 3655 LWRK1 = LWORK - KEND1 3656C 3657 IF (LWRK1 .LT. 0) THEN 3658 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 3659 CALL QUIT('Insufficient memory in CCFRETIJ') 3660 ENDIF 3661C 3662 CALL DZERO(WORK(KSCR1),LSCR) 3663C 3664 KOFF3 = ICOFRO(ISYMI,ISYMK) + 1 3665 KOFF4 = IT1FRO(ISYMC,ISYMK) + 1 3666 KOFFT = IT1AM(ISYMC,ISYMJ) + 1 3667 KOFFR = IMATIJ(ISYMI,ISYMJ) + 1 3668C 3669 NTOTI = MAX(NRHF(ISYMI),1) 3670 NTOTC = MAX(NVIR(ISYMC),1) 3671C 3672 CALL DGEMM('N','T',NRHF(ISYMI),NVIR(ISYMC),NRHFFR(ISYMK), 3673 * ONE,DFJI(KOFF3),NTOTI,XFIA(KOFF4),NTOTC,ZERO, 3674 * WORK(KSCR1),NTOTI) 3675C 3676 CALL DGEMM('N','T',NRHF(ISYMI),NVIR(ISYMC),NRHFFR(ISYMK), 3677 * -ONE,XFIJ(KOFF3),NTOTI,DFAI(KOFF4),NTOTC,ONE, 3678 * WORK(KSCR1),NTOTI) 3679C 3680 CALL DGEMM('N','N',NRHF(ISYMI),NRHF(ISYMJ),NVIR(ISYMC), 3681 * ONE,WORK(KSCR1),NTOTI,T1AM(KOFFT),NTOTC,ONE, 3682 * ETAIJ(KOFFR),NTOTI) 3683C 3684 110 CONTINUE 3685C 3686C------------------------------------- 3687C Do the last two terms with T1AM. 3688C------------------------------------- 3689C 3690 DO 120 ISYMI = 1,NSYM 3691C 3692 ISYMJ = ISYMI 3693 ISYMC = ISYMI 3694 ISYMK = MULD2H(ISYMJ,ISYM) 3695C 3696 LSCR = NRHF(ISYMJ)*NVIR(ISYMC) 3697C 3698 KSCR2 = 1 3699 KEND2 = KSCR2 + LSCR 3700 LWRK2 = LWORK - KEND2 3701C 3702 IF (LWRK2 .LT. 0) THEN 3703 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND2 3704 CALL QUIT('Insufficient memory in CCFRETIJ') 3705 ENDIF 3706C 3707 CALL DZERO(WORK(KSCR2),LSCR) 3708C 3709 KOFF5 = IT1FRO(ISYMC,ISYMK) + 1 3710 KOFF6 = ICOFRO(ISYMJ,ISYMK) + 1 3711 KOFFT = IT1AM(ISYMC,ISYMI) + 1 3712 KOFFR = IMATIJ(ISYMI,ISYMJ) + 1 3713C 3714 NTOTC = MAX(NVIR(ISYMC),1) 3715 NTOTJ = MAX(NRHF(ISYMJ),1) 3716 NTOTI = MAX(NRHF(ISYMI),1) 3717C 3718 CALL DGEMM('N','T',NVIR(ISYMC),NRHF(ISYMJ),NRHFFR(ISYMK), 3719 * ONE,XFIA(KOFF5),NTOTC,DFJI(KOFF6),NTOTJ,ZERO, 3720 * WORK(KSCR2),NTOTC) 3721C 3722 CALL DGEMM('N','T',NVIR(ISYMC),NRHF(ISYMJ),NRHFFR(ISYMK), 3723 * -ONE,DFAI(KOFF5),NTOTC,XFIJ(KOFF6),NTOTJ,ONE, 3724 * WORK(KSCR2),NTOTC) 3725C 3726 CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),NVIR(ISYMC), 3727 * -ONE,T1AM(KOFFT),NTOTC,WORK(KSCR2),NTOTC,ONE, 3728 * ETAIJ(KOFFR),NTOTI) 3729C 3730 120 CONTINUE 3731C 3732 RETURN 3733 END 3734C /* Deck ccfd1ii */ 3735 SUBROUTINE CCFD1II(DFII) 3736C 3737C Written by Asger Halkier 14/8 - 1998 3738C 3739C Version: 1.0 3740C 3741C Purpose: To calculate the simple frozen-frozen 1 electron density. 3742C 3743#include "implicit.h" 3744 PARAMETER(TWO = 2.0D0) 3745 DIMENSION DFII(*) 3746#include "maxorb.h" 3747#include "ccsdinp.h" 3748#include "ccorb.h" 3749#include "ccsdsym.h" 3750#include "ccfro.h" 3751C 3752 CALL DZERO(DFII,NFROFR(1)) 3753C 3754 DO 100 ISYM = 1,NSYM 3755 DO 110 I = 1,NRHFFR(ISYM) 3756 NII = IFROFR(ISYM,ISYM) + NRHFFR(ISYM)*(I - 1) + I 3757 DFII(NII) = TWO 3758 110 CONTINUE 3759 100 CONTINUE 3760C 3761 RETURN 3762 END 3763C /* Deck cc_etijf */ 3764 SUBROUTINE CC_ETIJF(ETAIJ,DIJ,DAB,DAI,DIA,DIJT,DFII,DFIJ,DFJI, 3765 * DFAI,DFIA,XIJ,XAI,XIA,XIJT,XABT,XFIJ,XFJI, 3766 * XFAI,XFIA,XFII,T1AM,WORK,LWORK,ISYM,IOPT) 3767C 3768C Written by Asger Halkier 14/8 - 1998 3769C 3770C Version: 1.0 3771C 3772C Purpose: To calculate eta-iJ-0. ISYM is the symmetry of the 3773C integrals and densities and IOPT controls the one 3774C and two electron contributions. 3775C 3776#include "implicit.h" 3777 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 3778 DIMENSION ETAIJ(*), DIJ(*), DAB(*), DAI(*), DIA(*), DIJT(*) 3779 DIMENSION DFII(*), DFIJ(*), DFJI(*), DFAI(*), DFIA(*), XIJ(*) 3780 DIMENSION XAI(*), XIA(*), XIJT(*), XABT(*), XFIJ(*), XFJI(*) 3781 DIMENSION XFAI(*), XFIA(*), XFII(*), T1AM(*), WORK(LWORK) 3782#include "priunit.h" 3783#include "maxorb.h" 3784#include "ccsdinp.h" 3785#include "ccorb.h" 3786#include "ccsdsym.h" 3787#include "ccfro.h" 3788C 3789C------------------------------- 3790C Do the direct terms first. 3791C------------------------------- 3792C 3793 DO 100 ISYMI = 1,NSYM 3794C 3795 ISYMJ = ISYMI 3796 ISYMC = MULD2H(ISYMI,ISYM) 3797 ISYMK = ISYMC 3798C 3799 IF (NRHFFR(ISYMJ) .EQ. 0 ) GOTO 100 3800C 3801 KOFFR = ICOFRO(ISYMI,ISYMJ) + 1 3802 KOFF1 = IT1AM(ISYMC,ISYMI) + 1 3803 KOFF2 = IT1FRO(ISYMC,ISYMJ) + 1 3804C 3805 NTOTI = MAX(NRHF(ISYMI),1) 3806 NTOTC = MAX(NVIR(ISYMC),1) 3807C 3808 CALL DGEMM('T','N',NRHF(ISYMI),NRHFFR(ISYMJ),NVIR(ISYMC), 3809 * -ONE,DIA(KOFF1),NTOTC,XFIA(KOFF2),NTOTC,ONE, 3810 * ETAIJ(KOFFR),NTOTI) 3811C 3812 CALL DGEMM('T','N',NRHF(ISYMI),NRHFFR(ISYMJ),NVIR(ISYMC), 3813 * -ONE,DAI(KOFF1),NTOTC,XFAI(KOFF2),NTOTC,ONE, 3814 * ETAIJ(KOFFR),NTOTI) 3815C 3816 KOFF3 = IMATIJ(ISYMI,ISYMK) + 1 3817 KOFF4 = ICOFRO(ISYMK,ISYMJ) + 1 3818C 3819 NTOTK = MAX(NRHF(ISYMK),1) 3820C 3821 CALL DGEMM('N','N',NRHF(ISYMI),NRHFFR(ISYMJ),NRHF(ISYMK), 3822 * -ONE,DIJ(KOFF3),NTOTI,XFJI(KOFF4),NTOTK,ONE, 3823 * ETAIJ(KOFFR),NTOTI) 3824C 3825 CALL DGEMM('N','N',NRHF(ISYMI),NRHFFR(ISYMJ),NRHF(ISYMK), 3826 * -ONE,DIJT(KOFF3),NTOTI,XFIJ(KOFF4),NTOTK,ONE, 3827 * ETAIJ(KOFFR),NTOTI) 3828C 3829 KOFF5 = ICOFRO(ISYMI,ISYMK) + 1 3830 KOFF6 = IFROFR(ISYMK,ISYMJ) + 1 3831 KOFF7 = IFROFR(ISYMJ,ISYMK) + 1 3832C 3833 NTOTK = MAX(NRHFFR(ISYMK),1) 3834 NTOTJ = MAX(NRHFFR(ISYMJ),1) 3835C 3836 CALL DGEMM('N','N',NRHF(ISYMI),NRHFFR(ISYMJ),NRHFFR(ISYMK), 3837 * ONE,XFJI(KOFF5),NTOTI,DFII(KOFF6),NTOTK,ONE, 3838 * ETAIJ(KOFFR),NTOTI) 3839C 3840 CALL DGEMM('N','T',NRHF(ISYMI),NRHFFR(ISYMJ),NRHFFR(ISYMK), 3841 * ONE,XFIJ(KOFF5),NTOTI,DFII(KOFF7),NTOTJ,ONE, 3842 * ETAIJ(KOFFR),NTOTI) 3843C 3844C------------------------------------------------------------ 3845C The terms that only has contributions from 2e- part. 3846C------------------------------------------------------------ 3847C 3848 IF (IOPT .EQ. 2) THEN 3849C 3850 CALL DGEMM('T','N',NRHF(ISYMI),NRHFFR(ISYMJ),NVIR(ISYMC), 3851 * ONE,XAI(KOFF1),NTOTC,DFAI(KOFF2),NTOTC,ONE, 3852 * ETAIJ(KOFFR),NTOTI) 3853C 3854 CALL DGEMM('T','N',NRHF(ISYMI),NRHFFR(ISYMJ),NVIR(ISYMC), 3855 * ONE,XIA(KOFF1),NTOTC,DFIA(KOFF2),NTOTC,ONE, 3856 * ETAIJ(KOFFR),NTOTI) 3857C 3858 NTOTK = MAX(NRHF(ISYMK),1) 3859C 3860 CALL DGEMM('N','N',NRHF(ISYMI),NRHFFR(ISYMJ),NRHF(ISYMK), 3861 * ONE,XIJT(KOFF3),NTOTI,DFIJ(KOFF4),NTOTK,ONE, 3862 * ETAIJ(KOFFR),NTOTI) 3863C 3864 CALL DGEMM('N','N',NRHF(ISYMI),NRHFFR(ISYMJ),NRHF(ISYMK), 3865 * ONE,XIJ(KOFF3),NTOTI,DFJI(KOFF4),NTOTK,ONE, 3866 * ETAIJ(KOFFR),NTOTI) 3867C 3868 NTOTK = MAX(NRHFFR(ISYMK),1) 3869C 3870 CALL DGEMM('N','T',NRHF(ISYMI),NRHFFR(ISYMJ),NRHFFR(ISYMK), 3871 * -ONE,DFIJ(KOFF5),NTOTI,XFII(KOFF7),NTOTJ,ONE, 3872 * ETAIJ(KOFFR),NTOTI) 3873C 3874 CALL DGEMM('N','N',NRHF(ISYMI),NRHFFR(ISYMJ),NRHFFR(ISYMK), 3875 * -ONE,DFJI(KOFF5),NTOTI,XFII(KOFF6),NTOTK,ONE, 3876 * ETAIJ(KOFFR),NTOTI) 3877C 3878 ENDIF 3879 100 CONTINUE 3880C 3881C--------------------------------- 3882C Do the terms involving T1AM. 3883C--------------------------------- 3884C 3885 DO 110 ISYMI = 1,NSYM 3886C 3887 ISYMJ = ISYMI 3888 ISYMC = ISYMI 3889 ISYMD = MULD2H(ISYMJ,ISYM) 3890 ISYMK = MULD2H(ISYMJ,ISYM) 3891C 3892 IF (NRHFFR(ISYMJ) .EQ. 0 ) GOTO 110 3893C 3894 LSCR = NRHFFR(ISYMJ)*NVIR(ISYMC) 3895C 3896 KSCR1 = 1 3897 KEND1 = KSCR1 + LSCR 3898 LWRK1 = LWORK - KEND1 3899C 3900 IF (LWRK1 .LT. 0) THEN 3901 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 3902 CALL QUIT('Insufficient memory in CC_ETIJF') 3903 ENDIF 3904C 3905 CALL DZERO(WORK(KSCR1),LSCR) 3906C 3907 KOFFT = IT1AM(ISYMC,ISYMI) + 1 3908 KOFFR = ICOFRO(ISYMI,ISYMJ) + 1 3909C 3910 NTOTC = MAX(NVIR(ISYMC),1) 3911 NTOTI = MAX(NRHF(ISYMI),1) 3912 NTOTD = MAX(NVIR(ISYMD),1) 3913 NTOTK = MAX(NRHF(ISYMK),1) 3914 NTOKF = MAX(NRHFFR(ISYMK),1) 3915C 3916 KOFF8 = IMATAB(ISYMC,ISYMD) + 1 3917 KOFF9 = IT1FRO(ISYMD,ISYMJ) + 1 3918 KOFF10 = IT1AM(ISYMC,ISYMK) + 1 3919 KOFF11 = ICOFRO(ISYMK,ISYMJ) + 1 3920 KOFF12 = IT1FRO(ISYMC,ISYMK) + 1 3921 KOFF13 = IFROFR(ISYMK,ISYMJ) + 1 3922C 3923 CALL DGEMM('N','N',NVIR(ISYMC),NRHFFR(ISYMJ),NVIR(ISYMD), 3924 * -ONE,DAB(KOFF8),NTOTC,XFIA(KOFF9),NTOTD,ZERO, 3925 * WORK(KSCR1),NTOTC) 3926C 3927 CALL DGEMM('N','N',NVIR(ISYMC),NRHFFR(ISYMJ),NRHF(ISYMK), 3928 * -ONE,DAI(KOFF10),NTOTC,XFJI(KOFF11),NTOTK,ONE, 3929 * WORK(KSCR1),NTOTC) 3930C 3931 CALL DGEMM('N','N',NVIR(ISYMC),NRHFFR(ISYMJ),NRHFFR(ISYMK), 3932 * ONE,XFIA(KOFF12),NTOTC,DFII(KOFF13),NTOKF,ONE, 3933 * WORK(KSCR1),NTOTC) 3934C 3935C------------------------------------------------------------ 3936C The terms that only has contributions from 2e- part. 3937C------------------------------------------------------------ 3938C 3939 IF (IOPT .EQ. 2) THEN 3940C 3941 CALL DGEMM('N','N',NVIR(ISYMC),NRHFFR(ISYMJ),NVIR(ISYMD), 3942 * ONE,XABT(KOFF8),NTOTC,DFAI(KOFF9),NTOTD,ONE, 3943 * WORK(KSCR1),NTOTC) 3944C 3945 CALL DGEMM('N','N',NVIR(ISYMC),NRHFFR(ISYMJ),NRHF(ISYMK), 3946 * ONE,XIA(KOFF10),NTOTC,DFIJ(KOFF11),NTOTK,ONE, 3947 * WORK(KSCR1),NTOTC) 3948C 3949 KOFF14 = IFROFR(ISYMJ,ISYMK) + 1 3950 NTOTJ = MAX(NRHFFR(ISYMJ),1) 3951C 3952 CALL DGEMM('N','T',NVIR(ISYMC),NRHFFR(ISYMJ),NRHFFR(ISYMK), 3953 * -ONE,DFAI(KOFF12),NTOTC,XFII(KOFF14),NTOTJ,ONE, 3954 * WORK(KSCR1),NTOTC) 3955C 3956 ENDIF 3957C 3958 CALL DGEMM('T','N',NRHF(ISYMI),NRHFFR(ISYMJ),NVIR(ISYMC), 3959 * -ONE,T1AM(KOFFT),NTOTC,WORK(KSCR1),NTOTC,ONE, 3960 * ETAIJ(KOFFR),NTOTI) 3961C 3962 110 CONTINUE 3963C 3964 RETURN 3965 END 3966C /* Deck cc_etifjf */ 3967 SUBROUTINE CC_ETIFJF(ETAIJ,DFII,DFIJ,DFJI,DFAI,DFIA,XFIJ,XFJI, 3968 * XFAI,XFIA,XFII,ISYM,IOPT) 3969C 3970C Written by Asger Halkier 15/8 - 1998 3971C 3972C Version: 1.0 3973C 3974C Purpose: To calculate eta-IJ-0. ISYM is the symmetry of the 3975C integrals and densities and IOPT controls the one 3976C and two electron contributions. 3977C 3978#include "implicit.h" 3979 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 3980 DIMENSION ETAIJ(*), DFII(*), DFIJ(*), DFJI(*), DFAI(*), DFIA(*) 3981 DIMENSION XFIJ(*), XFJI(*), XFAI(*), XFIA(*), XFII(*) 3982#include "priunit.h" 3983#include "maxorb.h" 3984#include "ccsdinp.h" 3985#include "ccorb.h" 3986#include "ccsdsym.h" 3987#include "ccfro.h" 3988C 3989 DO 100 ISYMI = 1,NSYM 3990C 3991 ISYMJ = ISYMI 3992 ISYMK = MULD2H(ISYMI,ISYM) 3993 ISYMC = MULD2H(ISYMI,ISYM) 3994C 3995 IF ((NRHFFR(ISYMI) .EQ. 0).OR.(NRHFFR(ISYMJ) .EQ. 0)) GOTO 100 3996C 3997 KOFFR = IFROFR(ISYMI,ISYMJ) + 1 3998 KOFF1 = IFROFR(ISYMK,ISYMI) + 1 3999 KOFF2 = IFROFR(ISYMK,ISYMJ) + 1 4000 KOFF3 = IFROFR(ISYMI,ISYMK) + 1 4001 KOFF4 = IFROFR(ISYMJ,ISYMK) + 1 4002C 4003 NTOTI = MAX(NRHFFR(ISYMI),1) 4004 NTOTK = MAX(NRHFFR(ISYMK),1) 4005 NTOTJ = MAX(NRHFFR(ISYMJ),1) 4006C 4007 CALL DGEMM('T','N',NRHFFR(ISYMI),NRHFFR(ISYMJ),NRHFFR(ISYMK), 4008 * ONE,XFII(KOFF1),NTOTK,DFII(KOFF2),NTOTK,ONE, 4009 * ETAIJ(KOFFR),NTOTI) 4010C 4011 CALL DGEMM('N','T',NRHFFR(ISYMI),NRHFFR(ISYMJ),NRHFFR(ISYMK), 4012 * -ONE,DFII(KOFF3),NTOTI,XFII(KOFF4),NTOTJ,ONE, 4013 * ETAIJ(KOFFR),NTOTI) 4014C 4015 CALL DGEMM('T','N',NRHFFR(ISYMI),NRHFFR(ISYMJ),NRHFFR(ISYMK), 4016 * -ONE,DFII(KOFF1),NTOTK,XFII(KOFF2),NTOTK,ONE, 4017 * ETAIJ(KOFFR),NTOTI) 4018C 4019 CALL DGEMM('N','T',NRHFFR(ISYMI),NRHFFR(ISYMJ),NRHFFR(ISYMK), 4020 * ONE,XFII(KOFF3),NTOTI,DFII(KOFF4),NTOTJ,ONE, 4021 * ETAIJ(KOFFR),NTOTI) 4022C 4023C------------------------------------------------------------ 4024C The terms that only has contributions from 2e- part. 4025C------------------------------------------------------------ 4026C 4027 IF (IOPT .EQ. 2) THEN 4028C 4029 KOFF5 = ICOFRO(ISYMK,ISYMI) + 1 4030 KOFF6 = ICOFRO(ISYMK,ISYMJ) + 1 4031C 4032 NTOTK = MAX(NRHF(ISYMK),1) 4033C 4034 CALL DGEMM('T','N',NRHFFR(ISYMI),NRHFFR(ISYMJ),NRHF(ISYMK), 4035 * ONE,XFIJ(KOFF5),NTOTK,DFIJ(KOFF6),NTOTK,ONE, 4036 * ETAIJ(KOFFR),NTOTI) 4037C 4038 CALL DGEMM('T','N',NRHFFR(ISYMI),NRHFFR(ISYMJ),NRHF(ISYMK), 4039 * -ONE,DFJI(KOFF5),NTOTK,XFJI(KOFF6),NTOTK,ONE, 4040 * ETAIJ(KOFFR),NTOTI) 4041C 4042 CALL DGEMM('T','N',NRHFFR(ISYMI),NRHFFR(ISYMJ),NRHF(ISYMK), 4043 * -ONE,DFIJ(KOFF5),NTOTK,XFIJ(KOFF6),NTOTK,ONE, 4044 * ETAIJ(KOFFR),NTOTI) 4045C 4046 CALL DGEMM('T','N',NRHFFR(ISYMI),NRHFFR(ISYMJ),NRHF(ISYMK), 4047 * ONE,XFJI(KOFF5),NTOTK,DFJI(KOFF6),NTOTK,ONE, 4048 * ETAIJ(KOFFR),NTOTI) 4049C 4050 KOFF7 = IT1FRO(ISYMC,ISYMI) + 1 4051 KOFF8 = IT1FRO(ISYMC,ISYMJ) + 1 4052C 4053 NTOTC = MAX(NVIR(ISYMC),1) 4054C 4055 CALL DGEMM('T','N',NRHFFR(ISYMI),NRHFFR(ISYMJ),NVIR(ISYMC), 4056 * ONE,XFAI(KOFF7),NTOTC,DFAI(KOFF8),NTOTC,ONE, 4057 * ETAIJ(KOFFR),NTOTI) 4058C 4059 CALL DGEMM('T','N',NRHFFR(ISYMI),NRHFFR(ISYMJ),NVIR(ISYMC), 4060 * -ONE,DFIA(KOFF7),NTOTC,XFIA(KOFF8),NTOTC,ONE, 4061 * ETAIJ(KOFFR),NTOTI) 4062C 4063 CALL DGEMM('T','N',NRHFFR(ISYMI),NRHFFR(ISYMJ),NVIR(ISYMC), 4064 * -ONE,DFAI(KOFF7),NTOTC,XFAI(KOFF8),NTOTC,ONE, 4065 * ETAIJ(KOFFR),NTOTI) 4066C 4067 CALL DGEMM('T','N',NRHFFR(ISYMI),NRHFFR(ISYMJ),NVIR(ISYMC), 4068 * ONE,XFIA(KOFF7),NTOTC,DFIA(KOFF8),NTOTC,ONE, 4069 * ETAIJ(KOFFR),NTOTI) 4070C 4071 ENDIF 4072 100 CONTINUE 4073C 4074 RETURN 4075 END 4076C /* Deck cc_zkfcb */ 4077 SUBROUTINE CC_ZKFCB(ZKDIA,WORK,LWORK) 4078C 4079C Written by Asger Halkier 15/8 - 1998. 4080C 4081C To calculate the frozen-occupied block (iJ) of kappa-bar-0. 4082C 4083#include "implicit.h" 4084C 4085 PARAMETER (ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0) 4086 DIMENSION ZKDIA(*), WORK(LWORK) 4087C 4088#include "priunit.h" 4089#include "maxorb.h" 4090#include "ccsdinp.h" 4091#include "ccorb.h" 4092#include "ccsdsym.h" 4093#include "ccfro.h" 4094C 4095C--------------------------- 4096C Work space allocation. 4097C--------------------------- 4098C 4099 KETAIJ = 1 4100 KFOCKD = KETAIJ + NCOFRO(1) 4101 KEND1 = KFOCKD + NORBTS 4102 LWRK1 = LWORK - KEND1 4103C 4104 IF (LWRK1 .LT. 0) THEN 4105 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 4106 CALL QUIT('Insufficient work space for allocation in CC_ZKFCB') 4107 ENDIF 4108C 4109 CALL DZERO(WORK(KETAIJ),NCOFRO(1)) 4110 CALL DZERO(WORK(KFOCKD),NORBTS) 4111 CALL DCOPY(NCOFRO(1),ZKDIA(1),1,WORK(KETAIJ),1) 4112 CALL DZERO(ZKDIA(1),NCOFRO(1)) 4113C 4114C-------------------------- 4115C Get orbital energies. 4116C-------------------------- 4117C 4118 CALL FOCK_ALL(WORK(KFOCKD),WORK(KEND1),LWRK1) 4119C 4120C--------------------------- 4121C Calculate the results. 4122C--------------------------- 4123C 4124 DO 130 ISYMI = 1,NSYM 4125 ISYMJ = ISYMI 4126 DO 140 J = 1,NRHFFR(ISYMJ) 4127 DO 150 I = 1,NRHF(ISYMI) 4128 NIJ = ICOFRO(ISYMI,ISYMJ) + NRHF(ISYMI)*(J - 1) + I 4129 KOI = KFOCKD + IRHFS(ISYMI) + NRHFFR(ISYMI) + I - 1 4130 KOJ = KFOCKD + IRHFS(ISYMJ) + J - 1 4131C 4132 DENOM = WORK(KOJ) - WORK(KOI) 4133 ZKDIA(NIJ) = HALF*WORK(KETAIJ+NIJ-1)/DENOM 4134C 4135 150 CONTINUE 4136 140 CONTINUE 4137 130 CONTINUE 4138C 4139 CALL DCOPY(NCOFRO(1),ZKDIA(1),1,ZKDIA(1+NCOFRO(1)),1) 4140C 4141 RETURN 4142 END 4143C /* Deck cc_etaif */ 4144 SUBROUTINE CC_ETAIF(ETAAI,DAB,DAI,DIA,DIJT,DABT,DFII,DFIJ,DFJI, 4145 * DFAI,DFIA,XIJ,XAB,XAI,XIA,XABT,XFIJ,XFJI, 4146 * XFAI,XFIA,XFII,T1AM,WORK,LWORK,ISYM,IOPT) 4147C 4148C Written by Asger Halkier 16/8 - 1998 4149C 4150C Version: 1.0 4151C 4152C Purpose: To calculate eta-aI-0. ISYM is the symmetry of the 4153C integrals and densities and IOPT controls the one 4154C and two electron contributions. 4155C 4156#include "implicit.h" 4157 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 4158 DIMENSION ETAAI(*), DAB(*), DAI(*), DIA(*), DIJT(*), DABT(*) 4159 DIMENSION DFII(*), DFIJ(*), DFJI(*), DFAI(*), DFIA(*), XIJ(*) 4160 DIMENSION XAB(*), XAI(*), XIA(*), XABT(*), XFIJ(*), XFJI(*) 4161 DIMENSION XFAI(*), XFIA(*), XFII(*), T1AM(*), WORK(LWORK) 4162#include "priunit.h" 4163#include "maxorb.h" 4164#include "ccsdinp.h" 4165#include "ccorb.h" 4166#include "ccsdsym.h" 4167#include "ccfro.h" 4168C 4169C------------------------------- 4170C Do the direct terms first. 4171C------------------------------- 4172C 4173 DO 100 ISYMA = 1,NSYM 4174C 4175 ISYMI = ISYMA 4176 ISYMC = MULD2H(ISYMA,ISYM) 4177 ISYMK = MULD2H(ISYMA,ISYM) 4178 ISYMJ = MULD2H(ISYMA,ISYM) 4179C 4180 IF (NRHFFR(ISYMI) .EQ. 0) GOTO 100 4181C 4182 KOFFR = IT1FRO(ISYMA,ISYMI) + 1 4183 KOFF1 = IMATAB(ISYMA,ISYMC) + 1 4184 KOFF2 = IT1FRO(ISYMC,ISYMI) + 1 4185 KOFF3 = IT1AM(ISYMA,ISYMK) + 1 4186 KOFF4 = ICOFRO(ISYMK,ISYMI) + 1 4187 KOFF5 = IT1FRO(ISYMA,ISYMJ) + 1 4188 KOFF6 = IFROFR(ISYMJ,ISYMI) + 1 4189 KOFF7 = IFROFR(ISYMI,ISYMJ) + 1 4190C 4191 NTOTA = MAX(NVIR(ISYMA),1) 4192 NTOTC = MAX(NVIR(ISYMC),1) 4193 NTOTK = MAX(NRHF(ISYMK),1) 4194 NTOTJ = MAX(NRHFFR(ISYMJ),1) 4195 NTOTI = MAX(NRHFFR(ISYMI),1) 4196C 4197 CALL DGEMM('N','N',NVIR(ISYMA),NRHFFR(ISYMI),NVIR(ISYMC), 4198 * -ONE,DAB(KOFF1),NTOTA,XFIA(KOFF2),NTOTC,ONE, 4199 * ETAAI(KOFFR),NTOTA) 4200C 4201 CALL DGEMM('N','N',NVIR(ISYMA),NRHFFR(ISYMI),NVIR(ISYMC), 4202 * -ONE,DABT(KOFF1),NTOTA,XFAI(KOFF2),NTOTC,ONE, 4203 * ETAAI(KOFFR),NTOTA) 4204C 4205 CALL DGEMM('N','N',NVIR(ISYMA),NRHFFR(ISYMI),NRHF(ISYMK), 4206 * -ONE,DAI(KOFF3),NTOTA,XFJI(KOFF4),NTOTK,ONE, 4207 * ETAAI(KOFFR),NTOTA) 4208C 4209 CALL DGEMM('N','N',NVIR(ISYMA),NRHFFR(ISYMI),NRHF(ISYMK), 4210 * -ONE,DIA(KOFF3),NTOTA,XFIJ(KOFF4),NTOTK,ONE, 4211 * ETAAI(KOFFR),NTOTA) 4212C 4213 CALL DGEMM('N','N',NVIR(ISYMA),NRHFFR(ISYMI),NRHFFR(ISYMJ), 4214 * ONE,XFIA(KOFF5),NTOTA,DFII(KOFF6),NTOTJ,ONE, 4215 * ETAAI(KOFFR),NTOTA) 4216C 4217 CALL DGEMM('N','T',NVIR(ISYMA),NRHFFR(ISYMI),NRHFFR(ISYMJ), 4218 * ONE,XFAI(KOFF5),NTOTA,DFII(KOFF7),NTOTI,ONE, 4219 * ETAAI(KOFFR),NTOTA) 4220C 4221C------------------------------------------------------------ 4222C The terms that only has contributions from 2e- part. 4223C------------------------------------------------------------ 4224C 4225 IF (IOPT .EQ. 2) THEN 4226C 4227 CALL DGEMM('N','N',NVIR(ISYMA),NRHFFR(ISYMI),NVIR(ISYMC), 4228 * ONE,XABT(KOFF1),NTOTA,DFAI(KOFF2),NTOTC,ONE, 4229 * ETAAI(KOFFR),NTOTA) 4230C 4231 CALL DGEMM('N','N',NVIR(ISYMA),NRHFFR(ISYMI),NVIR(ISYMC), 4232 * ONE,XAB(KOFF1),NTOTA,DFIA(KOFF2),NTOTC,ONE, 4233 * ETAAI(KOFFR),NTOTA) 4234C 4235 CALL DGEMM('N','N',NVIR(ISYMA),NRHFFR(ISYMI),NRHF(ISYMK), 4236 * ONE,XIA(KOFF3),NTOTA,DFIJ(KOFF4),NTOTK,ONE, 4237 * ETAAI(KOFFR),NTOTA) 4238C 4239 CALL DGEMM('N','N',NVIR(ISYMA),NRHFFR(ISYMI),NRHF(ISYMK), 4240 * ONE,XAI(KOFF3),NTOTA,DFJI(KOFF4),NTOTK,ONE, 4241 * ETAAI(KOFFR),NTOTA) 4242C 4243 CALL DGEMM('N','T',NVIR(ISYMA),NRHFFR(ISYMI),NRHFFR(ISYMJ), 4244 * -ONE,DFAI(KOFF5),NTOTA,XFII(KOFF7),NTOTI,ONE, 4245 * ETAAI(KOFFR),NTOTA) 4246C 4247 CALL DGEMM('N','N',NVIR(ISYMA),NRHFFR(ISYMI),NRHFFR(ISYMJ), 4248 * -ONE,DFIA(KOFF5),NTOTA,XFII(KOFF6),NTOTJ,ONE, 4249 * ETAAI(KOFFR),NTOTA) 4250C 4251 ENDIF 4252 100 CONTINUE 4253C 4254C--------------------------------- 4255C Do the terms involving T1AM. 4256C--------------------------------- 4257C 4258 DO 110 ISYMA = 1,NSYM 4259C 4260 ISYMI = ISYMA 4261 ISYMK = ISYMA 4262 ISYMC = MULD2H(ISYMI,ISYM) 4263 ISYML = MULD2H(ISYMI,ISYM) 4264 ISYMJ = MULD2H(ISYMI,ISYM) 4265C 4266 IF (NRHFFR(ISYMI) .EQ. 0) GOTO 110 4267C 4268 LSCR = NRHF(ISYMK)*NRHFFR(ISYMI) 4269C 4270 KSCR1 = 1 4271 KEND1 = KSCR1 + LSCR 4272 LWRK1 = LWORK - KEND1 4273C 4274 IF (LWRK1 .LT. 0) THEN 4275 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 4276 CALL QUIT('Insufficient memory in CC_ETAIF') 4277 ENDIF 4278C 4279 CALL DZERO(WORK(KSCR1),LSCR) 4280C 4281 KOFFR = IT1FRO(ISYMA,ISYMI) + 1 4282 KOFFT = IT1AM(ISYMA,ISYMK) + 1 4283 KOFF8 = IT1AM(ISYMC,ISYMK) + 1 4284 KOFF9 = IT1FRO(ISYMC,ISYMI) + 1 4285 KOFF10 = IMATIJ(ISYMK,ISYML) + 1 4286 KOFF11 = ICOFRO(ISYML,ISYMI) + 1 4287 KOFF12 = ICOFRO(ISYMK,ISYMJ) + 1 4288 KOFF13 = IFROFR(ISYMI,ISYMJ) + 1 4289C 4290 NTOTA = MAX(NVIR(ISYMA),1) 4291 NTOTC = MAX(NVIR(ISYMC),1) 4292 NTOTK = MAX(NRHF(ISYMK),1) 4293 NTOTL = MAX(NRHF(ISYML),1) 4294 NTOTI = MAX(NRHFFR(ISYMI),1) 4295C 4296 CALL DGEMM('T','N',NRHF(ISYMK),NRHFFR(ISYMI),NVIR(ISYMC), 4297 * ONE,DAI(KOFF8),NTOTC,XFAI(KOFF9),NTOTC,ZERO, 4298 * WORK(KSCR1),NTOTK) 4299C 4300 CALL DGEMM('N','N',NRHF(ISYMK),NRHFFR(ISYMI),NRHF(ISYML), 4301 * ONE,DIJT(KOFF10),NTOTK,XFIJ(KOFF11),NTOTL,ONE, 4302 * WORK(KSCR1),NTOTK) 4303C 4304 CALL DGEMM('N','T',NRHF(ISYMK),NRHFFR(ISYMI),NRHFFR(ISYMJ), 4305 * -ONE,XFIJ(KOFF12),NTOTK,DFII(KOFF13),NTOTI,ONE, 4306 * WORK(KSCR1),NTOTK) 4307C 4308C------------------------------------------------------------ 4309C The terms that only has contributions from 2e- part. 4310C------------------------------------------------------------ 4311C 4312 IF (IOPT .EQ. 2) THEN 4313C 4314 CALL DGEMM('T','N',NRHF(ISYMK),NRHFFR(ISYMI),NVIR(ISYMC), 4315 * -ONE,XIA(KOFF8),NTOTC,DFIA(KOFF9),NTOTC,ONE, 4316 * WORK(KSCR1),NTOTK) 4317C 4318 CALL DGEMM('N','N',NRHF(ISYMK),NRHFFR(ISYMI),NRHF(ISYML), 4319 * -ONE,XIJ(KOFF10),NTOTK,DFJI(KOFF11),NTOTL,ONE, 4320 * WORK(KSCR1),NTOTK) 4321C 4322 KOFF14 = IFROFR(ISYMJ,ISYMI) + 1 4323 NTOTJ = MAX(NRHFFR(ISYMJ),1) 4324C 4325 CALL DGEMM('N','N',NRHF(ISYMK),NRHFFR(ISYMI),NRHFFR(ISYMJ), 4326 * ONE,DFJI(KOFF12),NTOTK,XFII(KOFF14),NTOTJ,ONE, 4327 * WORK(KSCR1),NTOTK) 4328C 4329 ENDIF 4330C 4331 CALL DGEMM('N','N',NVIR(ISYMA),NRHFFR(ISYMI),NRHF(ISYMK), 4332 * -ONE,T1AM(KOFFT),NTOTA,WORK(KSCR1),NTOTK,ONE, 4333 * ETAAI(KOFFR),NTOTA) 4334C 4335 110 CONTINUE 4336C 4337 RETURN 4338 END 4339C /* Deck cc_frcoc1 */ 4340 SUBROUTINE CC_FRCOC1(XINTAI,XINOFO,ISYDIS) 4341C 4342C Written by Asger Halkier 17/8 - 1998. 4343C 4344C To calculate the coulomb part and the first exchange part of the 4345C intermediate needed for the frozen-correlated correction to eta-ai. 4346C 4347#include "implicit.h" 4348 PARAMETER (ZERO = 0.0D0, TWO = 2.0D0, EIGHT = 8.0D0) 4349 DIMENSION XINTAI(*), XINOFO(*) 4350#include "priunit.h" 4351#include "maxorb.h" 4352#include "ccsdinp.h" 4353#include "ccorb.h" 4354#include "ccsdsym.h" 4355#include "ccfro.h" 4356C 4357C-------------------------------------------------- 4358C Initial symmetries and work space allocation. 4359C-------------------------------------------------- 4360C 4361 ISYMA = ISYDIS 4362 ISYMI = ISYMA 4363 ISYMJK = 1 4364C 4365 IF (NRHF(ISYMI) .EQ. 0) RETURN 4366C 4367C------------------------------- 4368C Copy the coulomb part out. 4369C------------------------------- 4370C 4371 KOFF1 = IOFROO(ISYMJK,ISYMI) + 1 4372 KOFF2 = IOFROO(ISYMJK,ISYMI) + 1 4373 LEN = NCOFRO(ISYMJK)*NRHF(ISYMI) 4374 CALL DAXPY(LEN,EIGHT,XINOFO(KOFF1),1,XINTAI(KOFF2),1) 4375C 4376C-------------------------------- 4377C Do the first exchange part. 4378C-------------------------------- 4379C 4380 DO 100 ISYMJ = 1,NSYM 4381 ISYMK = MULD2H(ISYMJ,ISYMJK) 4382 ISYMIK = MULD2H(ISYMI,ISYMK) 4383 DO 110 J = 1,NRHF(ISYMJ) 4384 DO 120 K = 1,NRHFFR(ISYMK) 4385C 4386 KOFF2 = IOFROO(ISYMIK,ISYMJ) + NCOFRO(ISYMIK)*(J - 1) 4387 * + ICOFRO(ISYMI,ISYMK) + NRHF(ISYMI)*(K - 1) + 1 4388 KOFF3 = IOFROO(ISYMJK,ISYMI) + ICOFRO(ISYMJ,ISYMK) 4389 * + NRHF(ISYMJ)*(K - 1) + J 4390C 4391 CALL DAXPY(NRHF(ISYMI),-TWO,XINOFO(KOFF2),1, 4392 * XINTAI(KOFF3),NCOFRO(ISYMJK)) 4393C 4394 120 CONTINUE 4395 110 CONTINUE 4396 100 CONTINUE 4397C 4398 RETURN 4399 END 4400C /* Deck cc_frcoc1 */ 4401 SUBROUTINE CC_FRINDU(XINTAI,XINAIF,IDEL,ISYMD,LUN1,LUN2) 4402C 4403C Written by Asger Halkier 17/8 - 1998. 4404C 4405C To dump intermediates needed for the frozen-correlated corrections 4406C to eta-ai & eta-aI to disc. 4407C 4408#include "implicit.h" 4409 PARAMETER (ZERO = 0.0D0, TWO = 2.0D0, EIGHT = 8.0D0) 4410 DIMENSION XINTAI(*), XINAIF(*) 4411#include "priunit.h" 4412#include "maxorb.h" 4413#include "ccsdinp.h" 4414#include "ccorb.h" 4415#include "ccsdsym.h" 4416#include "ccfro.h" 4417C 4418 CHARACTER NAME1*8 4419 CHARACTER NAME2*8 4420C 4421 NAME1 = 'CCFRO1IN' 4422 NAME2 = 'CCFRO2IN' 4423C 4424C------------------------------------------ 4425C Write integral intermediates to disc. 4426C------------------------------------------ 4427C 4428 D = IDEL - IBAS(ISYMD) 4429 NTOT = NOFROO(ISYMD) 4430 IOFF = IOFOAO(ISYMD,ISYMD) + NTOT*(D - 1) + 1 4431C 4432 IF (NTOT .GT. 0) THEN 4433 CALL PUTWA2(LUN1,NAME1,XINTAI,IOFF,NTOT) 4434 ENDIF 4435C 4436 D = IDEL - IBAS(ISYMD) 4437 NTOT = NCOFRF(ISYMD) 4438 IOFF = IOFFAO(ISYMD,ISYMD) + NTOT*(D - 1) + 1 4439C 4440 IF (NTOT .GT. 0) THEN 4441 CALL PUTWA2(LUN2,NAME2,XINAIF,IOFF,NTOT) 4442 ENDIF 4443C 4444 RETURN 4445 END 4446C /* Deck cc_etacor */ 4447 SUBROUTINE CC_ETACOR(ETAAI,ETAAIF,ZETAJK,CMO,LUN1,LUN2,WORK,LWORK) 4448C 4449C Written by Asger Halkier 17/8 - 1998. 4450C 4451C To calculate the frozen-correlated corrections to eta-ai & 4452C eta-aI from intermediates on disc. 4453C 4454#include "implicit.h" 4455 PARAMETER (ZERO = 0.0D0, TWO = 2.0D0, ONE = 1.0D0) 4456 DIMENSION ETAAI(*), ETAAIF(*), ZETAJK(*), CMO(*), WORK(LWORK) 4457#include "priunit.h" 4458#include "maxorb.h" 4459#include "ccsdinp.h" 4460#include "ccorb.h" 4461#include "ccsdsym.h" 4462#include "ccfro.h" 4463C 4464 CHARACTER NAME1*8 4465 CHARACTER NAME2*8 4466C 4467 NAME1 = 'CCFRO1IN' 4468 NAME2 = 'CCFRO2IN' 4469C 4470C----------------------------- 4471C First we do the ai part. 4472C----------------------------- 4473C 4474 DO 100 ISYMD = 1,NSYM 4475C 4476 IF (NBAS(ISYMD) .EQ. 0) GOTO 100 4477C 4478 ISYJKI = ISYMD 4479 ISYMA = ISYMD 4480 ISYMI = ISYMA 4481 ISYMJK = 1 4482C 4483C---------------------------------- 4484C Work space allocation one. 4485C---------------------------------- 4486C 4487 KSCR1 = 1 4488 KINIM = KSCR1 + NBAS(ISYMD)*NOFROO(ISYJKI) 4489 KEND1 = KINIM + NVIR(ISYMA)*NOFROO(ISYJKI) 4490 LWRK1 = LWORK - KEND1 4491C 4492 IF (LWRK1 .LT. 0) THEN 4493 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 4494 CALL QUIT('Insufficient memory in CC_ETACOR') 4495 ENDIF 4496C 4497 CALL DZERO(WORK(KSCR1),KEND1) 4498C 4499C--------------------------------- 4500C Read integrals from disc. 4501C--------------------------------- 4502C 4503 NTOT = NOFROO(ISYJKI)*NBAS(ISYMD) 4504 IOFF = IOFOAO(ISYJKI,ISYMD) + 1 4505C 4506 IF (NTOT .GT. 0) THEN 4507 CALL GETWA2(LUN1,NAME1,WORK(KSCR1),IOFF,NTOT) 4508 END IF 4509C 4510C-------------------------------------- 4511C Transform AO index to virtual. 4512C-------------------------------------- 4513C 4514 KOFF1 = ILVISI(ISYMD) + 1 4515C 4516 NTOJKI = MAX(NOFROO(ISYJKI),1) 4517 NTOTD = MAX(NBAS(ISYMD),1) 4518C 4519 CALL DGEMM('N','N',NOFROO(ISYJKI),NVIR(ISYMA),NBAS(ISYMD), 4520 * ONE,WORK(KSCR1),NTOJKI,CMO(KOFF1),NTOTD,ZERO, 4521 * WORK(KINIM),NTOJKI) 4522C 4523C----------------------------------- 4524C Calculate the contribution. 4525C----------------------------------- 4526C 4527 DO 110 A = 1,NVIR(ISYMA) 4528 DO 120 I = 1,NRHF(ISYMI) 4529 KOFF2 = KINIM + NOFROO(ISYJKI)*(A - 1) 4530 * + IOFROO(ISYMJK,ISYMI) + NCOFRO(ISYMJK)*(I - 1) 4531 KOFF3 = IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I - 1) + A 4532 FACT = DDOT(NCOFRO(1),WORK(KOFF2),1,ZETAJK(1),1) 4533 ETAAI(KOFF3) = ETAAI(KOFF3) + FACT 4534C 4535 120 CONTINUE 4536 110 CONTINUE 4537 100 CONTINUE 4538C 4539C---------------------------- 4540C Then we do the aI part. 4541C---------------------------- 4542C 4543 DO 130 ISYMD = 1,NSYM 4544C 4545 IF (NBAS(ISYMD) .EQ. 0) GOTO 130 4546C 4547 ISYJKI = ISYMD 4548 ISYMA = ISYMD 4549 ISYMI = ISYMA 4550 ISYMJK = 1 4551C 4552C---------------------------------- 4553C Work space allocation two. 4554C---------------------------------- 4555C 4556 KSCR2 = 1 4557 KINIM = KSCR2 + NBAS(ISYMD)*NCOFRF(ISYJKI) 4558 KEND2 = KINIM + NVIR(ISYMA)*NCOFRF(ISYJKI) 4559 LWRK2 = LWORK - KEND2 4560C 4561 IF (LWRK2 .LT. 0) THEN 4562 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND2 4563 CALL QUIT('Insufficient memory in CC_ETACOR') 4564 ENDIF 4565C 4566 CALL DZERO(WORK(KSCR2),KEND2) 4567C 4568C--------------------------------- 4569C Read integrals from disc. 4570C--------------------------------- 4571C 4572 NTOT = NCOFRF(ISYJKI)*NBAS(ISYMD) 4573 IOFF = IOFFAO(ISYJKI,ISYMD) + 1 4574C 4575 IF (NTOT .GT. 0) THEN 4576 CALL GETWA2(LUN2,NAME2,WORK(KSCR2),IOFF,NTOT) 4577 END IF 4578C 4579C-------------------------------------- 4580C Transform AO index to virtual. 4581C-------------------------------------- 4582C 4583 KOFF4 = ILVISI(ISYMD) + 1 4584C 4585 NTOJKI = MAX(NCOFRF(ISYJKI),1) 4586 NTOTD = MAX(NBAS(ISYMD),1) 4587C 4588 CALL DGEMM('N','N',NCOFRF(ISYJKI),NVIR(ISYMA),NBAS(ISYMD), 4589 * ONE,WORK(KSCR2),NTOJKI,CMO(KOFF4),NTOTD,ZERO, 4590 * WORK(KINIM),NTOJKI) 4591C 4592C----------------------------------- 4593C Calculate the contribution. 4594C----------------------------------- 4595C 4596 DO 140 A = 1,NVIR(ISYMA) 4597 DO 150 I = 1,NRHFFR(ISYMI) 4598 KOFF5 = KINIM + NCOFRF(ISYJKI)*(A - 1) 4599 * + ICOFRF(ISYMJK,ISYMI) + NCOFRO(ISYMJK)*(I - 1) 4600 KOFF6 = IT1FRO(ISYMA,ISYMI) + NVIR(ISYMA)*(I - 1) + A 4601 FACT = DDOT(NCOFRO(1),WORK(KOFF5),1,ZETAJK(1),1) 4602 ETAAIF(KOFF6) = ETAAIF(KOFF6) + FACT 4603C 4604 150 CONTINUE 4605 140 CONTINUE 4606 130 CONTINUE 4607C 4608 RETURN 4609 END 4610C /* Deck cc_frcomi */ 4611 SUBROUTINE CC_FRCOMI(XINTAI,XINAIF,DSFRO,CMO,WORK,LWORK, 4612 * ISYDEL) 4613C 4614C Written by Asger Halkier 17/8 - 1998. 4615C 4616C To calculate the second exchange part of the intermediates needed 4617C for the frozen-correlated correction to eta-ai and eta-aI. 4618C 4619#include "implicit.h" 4620 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 4621 DIMENSION XINTAI(*), XINAIF(*), DSFRO(*), CMO(*), WORK(LWORK) 4622#include "priunit.h" 4623#include "maxorb.h" 4624#include "ccsdinp.h" 4625#include "ccorb.h" 4626#include "ccsdsym.h" 4627#include "ccfro.h" 4628C 4629C-------------------------------------------------- 4630C Initial symmetries and work space allocation. 4631C-------------------------------------------------- 4632C 4633 ISYMA = ISYDEL 4634 ISYMI = ISYMA 4635 ISYMJK = 1 4636C 4637C--------------------------------------------------------------------------- 4638C Loops over third index in integral batch & final wok space allocation. 4639C--------------------------------------------------------------------------- 4640C 4641 DO 100 ISYMK = 1,NSYM 4642C 4643 ISYMJ = MULD2H(ISYMJK,ISYMK) 4644 ISYMIJ = MULD2H(ISYMI,ISYMJ) 4645 ISALBE = ISYMIJ 4646 ISYMAL = ISYMI 4647 ISYMBE = ISYMJ 4648C 4649 KAOINT = 1 4650 KSCR1 = KAOINT + N2BST(ISALBE) 4651 KMOIN1 = KSCR1 + NBAS(ISYMAL)*NRHF(ISYMJ) 4652 KMOIN2 = KMOIN1 + NRHF(ISYMI)*NRHF(ISYMJ) 4653 KEND1 = KMOIN2 + NRHFFR(ISYMI)*NRHF(ISYMJ) 4654 LWRK1 = LWORK - KEND1 4655C 4656 IF (LWRK1 .LT. 0) THEN 4657 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 4658 CALL QUIT('Insufficient work space for allocation '// 4659 & 'in CC_FRCOMI') 4660 ENDIF 4661C 4662 CALL DZERO(WORK(KAOINT),N2BST(ISALBE)) 4663 CALL DZERO(WORK(KMOIN1),NRHF(ISYMI)*NRHF(ISYMJ)) 4664 CALL DZERO(WORK(KMOIN2),NRHFFR(ISYMI)*NRHF(ISYMJ)) 4665C 4666 DO 110 K = 1,NRHFFR(ISYMK) 4667C 4668C----------------------------------------------------------- 4669C Square up integral distribution (al be | K del). 4670C----------------------------------------------------------- 4671C 4672 KOFF1 = IDSFRO(ISALBE,ISYMK) + NNBST(ISALBE)*(K - 1) + 1 4673C 4674 CALL CCSD_SYMSQ(DSFRO(KOFF1),ISALBE,WORK(KAOINT)) 4675C 4676C--------------------------------------------- 4677C Calculate integrals (i j | K del). 4678C--------------------------------------------- 4679C 4680 KOFF2 = KAOINT + IAODIS(ISYMAL,ISYMBE) 4681 KOFF3 = ILRHSI(ISYMBE) + NBAS(ISYMBE)*NRHFFR(ISYMJ) + 1 4682C 4683 NTOTAL = MAX(NBAS(ISYMAL),1) 4684 NTOTBE = MAX(NBAS(ISYMBE),1) 4685C 4686 CALL DGEMM('N','N',NBAS(ISYMAL),NRHF(ISYMJ),NBAS(ISYMBE), 4687 * ONE,WORK(KOFF2),NTOTAL,CMO(KOFF3),NTOTBE,ZERO, 4688 * WORK(KSCR1),NTOTAL) 4689C 4690 KOFF4 = ILRHSI(ISYMAL) + NBAS(ISYMAL)*NRHFFR(ISYMI) + 1 4691C 4692 NTOTAL = MAX(NBAS(ISYMAL),1) 4693 NTOTI = MAX(NRHF(ISYMI),1) 4694C 4695 CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),NBAS(ISYMAL), 4696 * ONE,CMO(KOFF4),NTOTAL,WORK(KSCR1),NTOTAL,ZERO, 4697 * WORK(KMOIN1),NTOTI) 4698C 4699C----------------------------------------------- 4700C Add contribution to ai-intermediate. 4701C----------------------------------------------- 4702C 4703 DO 120 J = 1,NRHF(ISYMJ) 4704C 4705 KOFF5 = KMOIN1 + NRHF(ISYMI)*(J - 1) 4706 KOFF6 = IOFROO(ISYMJK,ISYMI) + ICOFRO(ISYMJ,ISYMK) 4707 * + NRHF(ISYMJ)*(K - 1) + J 4708 CALL DAXPY(NRHF(ISYMI),-TWO,WORK(KOFF5),1, 4709 * XINTAI(KOFF6),NCOFRO(ISYMJK)) 4710C 4711 120 CONTINUE 4712C 4713C--------------------------------------------- 4714C Calculate integrals (I j | K del). 4715C--------------------------------------------- 4716C 4717 KOFF7 = ILRHSI(ISYMAL) + 1 4718C 4719 NTOTAL = MAX(NBAS(ISYMAL),1) 4720 NTOTI = MAX(NRHFFR(ISYMI),1) 4721C 4722 CALL DGEMM('T','N',NRHFFR(ISYMI),NRHF(ISYMJ),NBAS(ISYMAL), 4723 * ONE,CMO(KOFF7),NTOTAL,WORK(KSCR1),NTOTAL,ZERO, 4724 * WORK(KMOIN2),NTOTI) 4725C 4726C----------------------------------------------- 4727C Add contribution to aI-intermediate. 4728C----------------------------------------------- 4729C 4730 DO 130 J = 1,NRHF(ISYMJ) 4731C 4732 KOFF8 = KMOIN2 + NRHFFR(ISYMI)*(J - 1) 4733 KOFF9 = ICOFRF(ISYMJK,ISYMI) + ICOFRO(ISYMJ,ISYMK) 4734 * + NRHF(ISYMJ)*(K - 1) + J 4735 CALL DAXPY(NRHFFR(ISYMI),-TWO,WORK(KOFF8),1, 4736 * XINAIF(KOFF9),NCOFRO(ISYMJK)) 4737C 4738 130 CONTINUE 4739 110 CONTINUE 4740 100 CONTINUE 4741C 4742 RETURN 4743 END 4744C /* Deck cc_frcof1 */ 4745 SUBROUTINE CC_FRCOF1(XINAIF,DSFRO,CMO,WORK,LWORK,ISYDEL) 4746C 4747C Written by Asger Halkier 17/8 - 1998. 4748C 4749C To calculate the coulomb part of the intermediates needed 4750C for the frozen-correlated correction to eta-aI. 4751C 4752#include "implicit.h" 4753 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, EIGHT = 8.0D0) 4754 DIMENSION XINAIF(*), DSFRO(*), CMO(*), WORK(LWORK) 4755#include "priunit.h" 4756#include "maxorb.h" 4757#include "ccsdinp.h" 4758#include "ccorb.h" 4759#include "ccsdsym.h" 4760#include "ccfro.h" 4761C 4762C-------------------------------------------------- 4763C Initial symmetries and work space allocation. 4764C-------------------------------------------------- 4765C 4766 ISYMA = ISYDEL 4767 ISYMI = ISYMA 4768 ISYMJK = 1 4769 ISALBE = ISYMJK 4770C 4771 IF (NRHFFR(ISYMI) .EQ. 0) RETURN 4772C 4773 KAOINT = 1 4774 KMOINT = KAOINT + N2BST(ISALBE) 4775 KEND1 = KMOINT + NCOFRO(ISYMJK) 4776 LWRK1 = LWORK - KEND1 4777C 4778 IF (LWRK1 .LT. 0) THEN 4779 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 4780 CALL QUIT('Insufficient work space for allocation '// 4781 & 'in CC_FRCOF1') 4782 ENDIF 4783C 4784 CALL DZERO(WORK(KAOINT),N2BST(ISALBE)) 4785C 4786 DO 100 I = 1,NRHFFR(ISYMI) 4787C 4788 CALL DZERO(WORK(KMOINT),NCOFRO(ISYMJK)) 4789C 4790C-------------------------------------------------------- 4791C Square up integral distribution (al be | I del). 4792C-------------------------------------------------------- 4793C 4794 KOFF1 = IDSFRO(ISALBE,ISYMI) + NNBST(ISALBE)*(I - 1) + 1 4795C 4796 CALL CCSD_SYMSQ(DSFRO(KOFF1),ISALBE,WORK(KAOINT)) 4797C 4798C------------------------------------------------------ 4799C Inner symmetry loop and final work allocation. 4800C------------------------------------------------------ 4801C 4802 DO 110 ISYMJ = 1,NSYM 4803C 4804 ISYMK = MULD2H(ISYMJK,ISYMJ) 4805 ISYMAL = ISYMJ 4806 ISYMBE = ISYMK 4807C 4808 KSCR1 = KEND1 4809 KEND2 = KSCR1 + NBAS(ISYMAL)*NRHFFR(ISYMK) 4810 LWRK2 = LWORK - KEND2 4811C 4812 IF (LWRK2 .LT. 0) THEN 4813 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND2 4814 CALL QUIT('Insufficient memory for allocation '// 4815 & 'in CC_FRCOF1') 4816 ENDIF 4817C 4818 CALL DZERO(WORK(KSCR1),NBAS(ISYMAL)*NRHFFR(ISYMK)) 4819C 4820C--------------------------------------------- 4821C Calculate integrals (j K | I del). 4822C--------------------------------------------- 4823C 4824 KOFF2 = KAOINT + IAODIS(ISYMAL,ISYMBE) 4825 KOFF3 = ILRHSI(ISYMBE) + 1 4826C 4827 NTOTAL = MAX(NBAS(ISYMAL),1) 4828 NTOTBE = MAX(NBAS(ISYMBE),1) 4829C 4830 CALL DGEMM('N','N',NBAS(ISYMAL),NRHFFR(ISYMK),NBAS(ISYMBE), 4831 * ONE,WORK(KOFF2),NTOTAL,CMO(KOFF3),NTOTBE,ZERO, 4832 * WORK(KSCR1),NTOTAL) 4833C 4834 KOFF4 = ILRHSI(ISYMAL) + NBAS(ISYMAL)*NRHFFR(ISYMJ) + 1 4835 KOFF5 = KMOINT + ICOFRO(ISYMJ,ISYMK) 4836C 4837 NTOTAL = MAX(NBAS(ISYMAL),1) 4838 NTOTJ = MAX(NRHF(ISYMJ),1) 4839C 4840 CALL DGEMM('T','N',NRHF(ISYMJ),NRHFFR(ISYMK),NBAS(ISYMAL), 4841 * ONE,CMO(KOFF4),NTOTAL,WORK(KSCR1),NTOTAL,ZERO, 4842 * WORK(KOFF5),NTOTJ) 4843C 4844 110 CONTINUE 4845C 4846C------------------------------------------------------- 4847C Calculate the contribution to the intermediate. 4848C------------------------------------------------------- 4849C 4850 KOFF6 = ICOFRF(ISYMJK,ISYMI) + NCOFRO(ISYMJK)*(I - 1) + 1 4851 CALL DAXPY(NCOFRO(ISYMJK),EIGHT,WORK(KMOINT),1,XINAIF(KOFF6),1) 4852C 4853 100 CONTINUE 4854C 4855 RETURN 4856 END 4857C /* Deck cc_frcof2 */ 4858 SUBROUTINE CC_FRCOF2(XINAIF,DSRHF,CMO,WORK,LWORK,ISYDEL) 4859C 4860C Written by Asger Halkier 17/8 - 1998. 4861C 4862C To calculate the first exchange part of the intermediates 4863C needed for the frozen-correlated correction to eta-aI. 4864C 4865#include "implicit.h" 4866 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 4867 DIMENSION XINAIF(*), DSRHF(*), CMO(*), WORK(LWORK) 4868#include "priunit.h" 4869#include "maxorb.h" 4870#include "ccsdinp.h" 4871#include "ccorb.h" 4872#include "ccsdsym.h" 4873#include "ccfro.h" 4874C 4875C-------------------------------------------------- 4876C Initial symmetries and work space allocation. 4877C-------------------------------------------------- 4878C 4879 ISYMA = ISYDEL 4880 ISYMI = ISYMA 4881 ISYMJK = 1 4882C 4883 IF (NRHFFR(ISYMI) .EQ. 0) RETURN 4884C 4885C--------------------------------------------------------------------- 4886C Loops over third index in integral batch & work space allocation. 4887C--------------------------------------------------------------------- 4888C 4889 DO 100 ISYMJ = 1,NSYM 4890C 4891 ISYMK = MULD2H(ISYMJK,ISYMJ) 4892 ISYMIK = MULD2H(ISYMI,ISYMK) 4893 ISALBE = ISYMIK 4894 ISYMAL = ISYMI 4895 ISYMBE = ISYMK 4896C 4897 IF (NRHFFR(ISYMK) .EQ. 0) GOTO 100 4898C 4899 KAOINT = 1 4900 KSCR1 = KAOINT + N2BST(ISALBE) 4901 KMOINT = KSCR1 + NBAS(ISYMAL)*NRHFFR(ISYMK) 4902 KEND1 = KMOINT + NRHFFR(ISYMI)*NRHFFR(ISYMK) 4903 LWRK1 = LWORK - KEND1 4904C 4905 IF (LWRK1 .LT. 0) THEN 4906 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 4907 CALL QUIT('Insufficient work space for allocation '// 4908 & 'in CC_FRCOF2') 4909 ENDIF 4910C 4911 CALL DZERO(WORK(KAOINT),KEND1) 4912C 4913 DO 110 J = 1,NRHF(ISYMJ) 4914C 4915C----------------------------------------------------------- 4916C Square up integral distribution (al be | j del). 4917C----------------------------------------------------------- 4918C 4919 KOFF1 = IDSRHF(ISALBE,ISYMJ) + NNBST(ISALBE)*(J - 1) + 1 4920C 4921 CALL CCSD_SYMSQ(DSRHF(KOFF1),ISALBE,WORK(KAOINT)) 4922C 4923C--------------------------------------------- 4924C Calculate integrals (I K | j del). 4925C--------------------------------------------- 4926C 4927 KOFF2 = KAOINT + IAODIS(ISYMAL,ISYMBE) 4928 KOFF3 = ILRHSI(ISYMBE) + 1 4929C 4930 NTOTAL = MAX(NBAS(ISYMAL),1) 4931 NTOTBE = MAX(NBAS(ISYMBE),1) 4932C 4933 CALL DGEMM('N','N',NBAS(ISYMAL),NRHFFR(ISYMK),NBAS(ISYMBE), 4934 * ONE,WORK(KOFF2),NTOTAL,CMO(KOFF3),NTOTBE,ZERO, 4935 * WORK(KSCR1),NTOTAL) 4936C 4937 KOFF4 = ILRHSI(ISYMAL) + 1 4938C 4939 NTOTAL = MAX(NBAS(ISYMAL),1) 4940 NTOTI = MAX(NRHFFR(ISYMI),1) 4941C 4942 CALL DGEMM('T','N',NRHFFR(ISYMI),NRHFFR(ISYMK), 4943 * NBAS(ISYMAL),ONE,CMO(KOFF4),NTOTAL,WORK(KSCR1), 4944 * NTOTAL,ZERO,WORK(KMOINT),NTOTI) 4945C 4946C---------------------------------------------------------- 4947C Calculate the contribution to the intermediate. 4948C---------------------------------------------------------- 4949C 4950 DO 120 K = 1,NRHFFR(ISYMK) 4951C 4952 KOFF5 = KMOINT + NRHFFR(ISYMI)*(K - 1) 4953 KOFF6 = ICOFRF(ISYMJK,ISYMI) + ICOFRO(ISYMJ,ISYMK) 4954 * + NRHF(ISYMJ)*(K - 1) + J 4955 CALL DAXPY(NRHFFR(ISYMI),-TWO,WORK(KOFF5),1, 4956 * XINAIF(KOFF6),NCOFRO(ISYMJK)) 4957C 4958 120 CONTINUE 4959 110 CONTINUE 4960 100 CONTINUE 4961C 4962 RETURN 4963 END 4964