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_d1ao */ 20 SUBROUTINE CC_D1AO(IPDD,R12PRP,AODEN,ZKAM,WORK,LWORK,MODEL,LLIST, 21 & ILLNR,NO,FNDPTIA, FNDPTIA2, FNDPTAB, FNDPTIJ) 22C 23C Written by Asger Halkier January 1996. 24C 25C Version: 1.0 26C 27C Purpose: To calculate the one electron Coupled Cluster 28C density matrix and returning it backtransformed 29C to AO-basis in AODEN! 30C 31C Current models: 32C MP2, CCD, CCSD, CCSD(T), CC3, DRCCD (DRPA), RCCD (singlet RPA) 33C 34C NO = .true. compute natural occupation numbers 35C 36C OC: 15-4-1997: general llist. 37C 38C May 1998: MP2 Frozen core density by Asger Halkier. 39C 40C Spring 2002 : CC3 section by K. Hald. 41C Summer 2010 : CCD + RCCD/DRCCD (RPA) session by Sonia Coriani 42C 43#include "implicit.h" 44#include "priunit.h" 45#include "dummy.h" 46#include "maxash.h" 47#include "maxorb.h" 48#include "mxcent.h" 49#include "aovec.h" 50#include "iratdef.h" 51 logical r12prp 52 PARAMETER (ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, TWO = 2.0D0) 53 DIMENSION INDEXA(MXCORB_CC) 54 DIMENSION AODEN(*), ZKAM(*), WORK(LWORK) 55#include "ccorb.h" 56#include "infind.h" 57#include "inftap.h" 58#include "blocks.h" 59#include "ccsdinp.h" 60#include "ccsdsym.h" 61#include "ccsdio.h" 62#include "distcl.h" 63#include "cbieri.h" 64#include "eribuf.h" 65#include "ccfop.h" 66#include "ccfro.h" 67C 68 CHARACTER MODEL*10,MODDUM*10 69 CHARACTER LLIST*(*) 70 CHARACTER*(*) FNDPTIA, FNDPTAB, FNDPTIJ, FNDPTIA2 71 LOGICAL NO 72C 73C---------------------------------------- 74C Initialization of timing parameter. 75C---------------------------------------- 76C 77 CALL QENTER('CC_D1AO') 78 TIMTOT = ZERO 79 TIMTOT = SECOND() 80C 81 IF ((CC2) .AND. (RELORB)) THEN 82C 83 CALL CC2_D1AO(AODEN,ZKAM,WORK,LWORK) 84 CALL QEXIT('CC_D1AO') 85 RETURN 86C 87 ELSE IF (MP2) THEN 88C 89 KONEAI = 1 90 KONEAB = KONEAI + NT1AMX 91 KONEIJ = KONEAB + NMATAB(1) 92 KONEIA = KONEIJ + NMATIJ(1) 93 KT1AM = KONEIA + NT1AMX 94 KLAMDH = KT1AM + NT1AMX 95 KLAMDP = KLAMDH + NLAMDT 96 KRMAT = KLAMDP + NLAMDT 97 KEND1 = KRMAT + NMATIJ(1) 98 LWRK1 = LWORK - KEND1 99C 100 IF (LWRK1 .LT. 0) THEN 101 WRITE(LUPRI,*) 'MP2 - Available:', LWORK, 'Needed:', KEND1 102 CALL QUIT('Insufficient memory for work '// 103 & 'allocation in CC_D1AO for MP2') 104 ENDIF 105C 106C----------------------------------------------------------------- 107C Initialize arrays (note that the t1-amplitudes are zero). 108C----------------------------------------------------------------- 109C 110 CALL DZERO(WORK(KONEAI),NT1AMX) 111 CALL DZERO(WORK(KONEAB),NMATAB(1)) 112 CALL DZERO(WORK(KONEIJ),NMATIJ(1)) 113 CALL DZERO(WORK(KONEIA),NT1AMX) 114 CALL DZERO(WORK(KT1AM),NT1AMX) 115 CALL DZERO(WORK(KRMAT),NMATIJ(1)) 116C 117C-------------------------------------------------- 118C Set up the relaxation part of the density. 119C-------------------------------------------------- 120C 121 CALL DCOPY(NMATIJ(1),ZKAM(1),1,WORK(KONEIJ),1) 122 CALL DCOPY(NMATAB(1),ZKAM(NMATIJ(1)+1),1,WORK(KONEAB),1) 123 CALL DCOPY(NT1AMX,ZKAM(NMATIJ(1)+NMATAB(1)+1),1,WORK(KONEAI),1) 124 CALL DCOPY(NT1AMX,WORK(KONEAI),1,WORK(KONEIA),1) 125C 126C------------------------------------- 127C Add the Hartree-Fock density. 128C Add R12/A' contribution (Elena) 129C------------------------------------- 130C 131celena 132 IF (R12PRP .AND. (IPDD .EQ. 3 .OR. IPDD .EQ. 5)) THEN 133 lunit = -1 134 IF (IPDD .EQ. 3) THEN 135 call gpopen(lunit,'CCR12YIJ','unknown',' ','unformatted', 136 & idummy,.false.) 137 ELSEIF (IPDD .EQ. 5) THEN 138 call gpopen(lunit,'CCR12YIJB','unknown',' ', 139 & 'unformatted',idummy,.false.) 140 ENDIF 141 142 do isymaj = 1,nsym 143 isymij = isymaj 144 ncvij = 0 145 do isyma = 1,nsym 146 isymj = muld2h(isymaj,isyma) 147 isymi = muld2h(isymij,isymj) 148 ncvij = ncvij + nrhf(isyma)*nrhf(isymi) 149 enddo 150 enddo 151 152 read(lunit)(work(krmat+i-1),i=1,ncvij) 153 call gpclose(lunit,'KEEP') 154 155 DO ISYMI = 1,NSYM 156 ISYMJ = ISYMI 157 DO J = 1,NRHF(ISYMJ) 158 DO I = J+1,NRHF(ISYMI) 159 NIJ = IMATIJ(ISYMI,ISYMJ) + NRHF(ISYMI) 160 * *(J - 1) + I 161 NJI = IMATIJ(ISYMJ,ISYMI) + NRHF(ISYMJ) 162 * *(I - 1) + J 163 WORK(KONEIJ + NIJ - 1) = WORK(KONEIJ + NIJ - 1) 164 * + work(krmat + nij - 1) 165 WORK(KONEIJ + NJI - 1) = WORK(KONEIJ + NJI - 1) 166 * + work(krmat + nji - 1) 167C 168 ENDDO 169 ENDDO 170 ENDDO 171 172 ENDIF 173 174celena 175 176 DO 80 ISYM = 1,NSYM 177 DO 90 I = 1,NRHF(ISYM) 178 NII = IMATIJ(ISYM,ISYM) + NRHF(ISYM)*(I - 1) + I 179 WORK(KONEIJ + NII - 1) = WORK(KONEIJ + NII - 1) + TWO 180 if (r12prp .and. (ipdd .eq. 3 .or. ipdd .eq. 5)) then 181 WORK(KONEIJ + NII - 1) = WORK(KONEIJ + NII - 1) + 182 & 1.0d0*work(krmat + nii - 1) 183 endif 184 90 CONTINUE 185 80 CONTINUE 186C 187C---------------------------------------- 188C Calculate MO coefficient matrix. 189C---------------------------------------- 190C 191 CALL LAMMAT(WORK(KLAMDP),WORK(KLAMDH),WORK(KT1AM), 192 * WORK(KEND1),LWRK1) 193C 194C-------------------------------------- 195C Transform density to AO basis. 196C-------------------------------------- 197C 198 CALL DZERO(AODEN,N2BST(1)) 199C 200 ISDEN = 1 201 CALL CC_DENAO(AODEN,1,WORK(KONEAI),WORK(KONEAB), 202 * WORK(KONEIJ),WORK(KONEIA),ISDEN,WORK(KLAMDP),1, 203 * WORK(KLAMDH),1,WORK(KEND1),LWRK1) 204C 205C---------------------------------------------------- 206C Add frozen core contributions to AO density. 207C---------------------------------------------------- 208C 209 IF (FROIMP) THEN 210C 211 KOFFAI = NMATIJ(1) + NMATAB(1) + 2*NT1AMX + 1 212 KOFFIA = KOFFAI + NT1FRO(1) 213 KOFFIJ = KOFFIA + NT1FRO(1) 214 KOFFJI = KOFFIJ + NCOFRO(1) 215C 216 ISDEN = 1 217 ICON = 1 218 CALL CC_D1FCB(AODEN,ZKAM(KOFFIJ),ZKAM(KOFFJI),ZKAM(KOFFAI), 219 * ZKAM(KOFFIA),WORK(KEND1),LWRK1,ISDEN,ICON) 220C 221 ENDIF 222C 223 ELSE IF ((CCSD .OR. CCPT .OR. CC3 .OR. CCD) 224 * .OR. ((CC2) .AND. (.NOT. RELORB))) THEN 225C 226C------------------------------------------------------- 227C Both zeta- and t-vectors are totally symmetric. 228C------------------------------------------------------- 229C 230 ISYMTR = 1 231 ISYMOP = 1 232C 233C-------------------------------------- 234C Initial work space allocation. 235C-------------------------------------- 236C 237 NHELP = NT2AMX 238C 239 KZ2AM = 1 240 KT2AM = KZ2AM + NT2SQ(1) 241 KT2AMT = KT2AM + NT2AMX 242 KLAMDP = KT2AMT + NHELP 243 KLAMDH = KLAMDP + NLAMDT 244 KT1AM = KLAMDH + NLAMDT 245 KZ1AM = KT1AM + NT1AMX 246 KEND1 = KZ1AM + NT1AMX 247 LWRK1 = LWORK - KEND1 248C 249 IF (LWRK1 .LT. 0) THEN 250 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 251 CALL QUIT('Insufficient core for initial '// 252 & 'allocation in CC_D1AO') 253 ENDIF 254C 255C------------------------------------------- 256C Read zero'th order zeta amplitudes. 257C------------------------------------------- 258C 259 ISYML = 1 260 IOPT = 3 261 CALL CC_RDRSP(LLIST,ILLNR,ISYML,IOPT,MODDUM, 262 & WORK(KZ1AM),WORK(KT2AM)) 263 IF ( IPRINT .GT. 10 ) THEN 264 RHO1N = DDOT(NT1AM(ISYML),WORK(KZ1AM),1,WORK(KZ1AM),1) 265 RHO2N = DDOT(NT2AM(ISYML),WORK(KT2AM),1,WORK(KT2AM),1) 266 WRITE(LUPRI,*) 'CC_D1AO: Norm of L. vector :',RHO1N 267 WRITE(LUPRI,*) 'CC_D1AO: Norm of L. vector :',RHO2N 268 ENDIF 269 270C 271C----------------------------------- 272C Square up zeta2 amplitudes. 273C----------------------------------- 274C 275 CALL CC_T2SQ(WORK(KT2AM),WORK(KZ2AM),1) 276C 277C---------------------------------------------- 278C Read zero'th order cluster amplitudes. 279C---------------------------------------------- 280C 281 IOPT = 3 282 CALL CC_RDRSP('R0',0,1,IOPT,MODEL,WORK(KT1AM),WORK(KT2AM)) 283C 284C--------------------------------------------------- 285C Zero out single vectors in CCD-calculation. 286C--------------------------------------------------- 287C 288 IF (CCD) THEN 289 CALL DZERO(WORK(KT1AM),NT1AMX) 290 CALL DZERO(WORK(KZ1AM),NT1AMX) 291 ENDIF 292C 293C------------------------------------- 294C Calculate the lamda matrices. 295C------------------------------------- 296C 297 CALL LAMMAT(WORK(KLAMDP),WORK(KLAMDH),WORK(KT1AM),WORK(KEND1), 298 * LWRK1) 299C 300C------------------------------------------ 301C Set up 2C-E of cluster amplitudes. 302C------------------------------------------ 303C 304 CALL DCOPY(NT2AMX,WORK(KT2AM),1,WORK(KT2AMT),1) 305 ISYOPE = 1 306 IOPTTCME = 1 307 CALL CCSD_TCMEPK(WORK(KT2AMT),1.0D0,ISYOPE,IOPTTCME) 308C 309C---------------------------------- 310C Work space allocation one. 311C Note that D(ai)=ZETA(ai) & 312C D(ia) is stored transposed 313C---------------------------------- 314C 315 KONEAI = KZ1AM 316 KONEAB = KONEAI + NT1AMX 317 KONEIJ = KONEAB + NMATAB(1) 318 KONEIA = KONEIJ + NMATIJ(1) 319 KXMAT = KONEIA + NT1AMX 320 KYMAT = KXMAT + NMATIJ(1) 321 KCMO = KYMAT + NMATAB(1) 322 KEND1 = KCMO + NLAMDS 323 LWRK1 = LWORK - KEND1 324C 325 IF (LWRK1 .LT. 0) THEN 326 WRITE(LUPRI,*) 'CC_D1AO; Available:',LWORK, 'Needed:',KEND1 327 CALL QUIT( 328 & 'Insufficient memory for first allocation in CC_D1AO') 329 ENDIF 330C 331C--------------------------------------------------------- 332C Initialize remaining one electron density arrays. 333C--------------------------------------------------------- 334C 335 CALL DZERO(WORK(KONEAB),NMATAB(1)) 336 CALL DZERO(WORK(KONEIJ),NMATIJ(1)) 337 CALL DZERO(WORK(KONEIA),NT1AMX) 338C 339C----------------------------------------------------------- 340C Calculate X-intermediate of zeta- and t-amplitudes. 341C----------------------------------------------------------- 342C 343 CALL CC_XI(WORK(KXMAT),WORK(KZ2AM),ISYMTR,WORK(KT2AM),ISYMOP, 344 * WORK(KEND1),LWRK1) 345C 346C----------------------------------------------------------- 347C Calculate Y-intermediate of zeta- and t-amplitudes. 348C----------------------------------------------------------- 349C 350 CALL CC_YI(WORK(KYMAT),WORK(KZ2AM),ISYMTR,WORK(KT2AM),ISYMOP, 351 * WORK(KEND1),LWRK1) 352C 353C----------------------------------------------------------------- 354C Construct three remaining blocks of one electron density. 355C----------------------------------------------------------------- 356C 357 CALL DCOPY(NMATAB(1),WORK(KYMAT),1,WORK(KONEAB),1) 358 CALL CC_EITR(WORK(KONEAB),WORK(KONEIJ),WORK(KEND1),LWRK1,1) 359 CALL DIJGEN(WORK(KONEIJ),WORK(KXMAT)) 360 CALL DIAGEN(WORK(KONEIA),WORK(KT2AMT),WORK(KONEAI)) 361C 362C------------------------------------------------------------ 363C Get the contribution from CC3. 364C------------------------------------------------------------ 365C 366 IF (CC3) THEN 367C 368 IF (LWRK1 .LT. MAX(NT1AMX,NMATIJ(1),NMATAB(1))) THEN 369 CALL QUIT('Out of memory in CC_D1AO (CC3)') 370 ENDIF 371C 372 LUPTIA = -1 373 CALL WOPEN2(LUPTIA,FNDPTIA,64,0) 374 IOFF = 1 + NT1AMX*(ILLNR-1) 375 CALL GETWA2(LUPTIA,FNDPTIA,WORK(KEND1),IOFF,NT1AMX) 376 CALL WCLOSE2(LUPTIA,FNDPTIA,'KEEP') 377 CALL DAXPY(NT1AMX,1.0D0,WORK(KEND1),1,WORK(KONEIA),1) 378C 379 LUPTAB = -1 380 CALL WOPEN2(LUPTAB,FNDPTAB,64,0) 381 IOFF = 1 + NMATAB(1)*(ILLNR-1) 382 CALL GETWA2(LUPTAB,FNDPTAB,WORK(KEND1),IOFF,NMATAB(1)) 383 CALL WCLOSE2(LUPTAB,FNDPTAB,'KEEP') 384 CALL DAXPY(NMATAB(1),1.0D0,WORK(KEND1),1,WORK(KONEAB),1) 385C 386 LUPTIJ = -1 387 CALL WOPEN2(LUPTIJ,FNDPTIJ,64,0) 388 IOFF = 1 + NMATIJ(1)*(ILLNR-1) 389 CALL GETWA2(LUPTIJ,FNDPTIJ,WORK(KEND1),IOFF,NMATIJ(1)) 390 CALL WCLOSE2(LUPTIJ,FNDPTIJ,'KEEP') 391 CALL DAXPY(NMATIJ(1),1.0D0,WORK(KEND1),1,WORK(KONEIJ),1) 392C 393 LUPTIA2 = -1 394 CALL WOPEN2(LUPTIA2,FNDPTIA2,64,0) 395 IOFF = 1 + NT1AMX*(ILLNR-1) 396 CALL GETWA2(LUPTIA2,FNDPTIA2,WORK(KEND1),IOFF,NT1AMX) 397 CALL WCLOSE2(LUPTIA2,FNDPTIA2,'KEEP') 398 CALL DAXPY(NT1AMX,1.0D0,WORK(KEND1),1,WORK(KONEIA),1) 399 ENDIF 400C 401C----------------------------------------------------------- 402C Backtransform the one electron density to AO-basis. 403C----------------------------------------------------------- 404C 405 CALL DZERO(AODEN,N2BST(1)) 406C 407 ISDEN = 1 408 CALL CC_DENAO(AODEN,1,WORK(KONEAI),WORK(KONEAB), 409 * WORK(KONEIJ),WORK(KONEIA),ISDEN,WORK(KLAMDP),1, 410 * WORK(KLAMDH),1,WORK(KEND1),LWRK1) 411C 412 IF (RELORB) THEN 413C 414C---------------------------------------------------------------- 415C Read MO-coefficients from interface file and reorder. 416C---------------------------------------------------------------- 417C 418 LUSIFC = -1 419 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY, 420 & .FALSE.) 421 REWIND LUSIFC 422 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 423 READ (LUSIFC) 424 READ (LUSIFC) 425 READ (LUSIFC) (WORK(KCMO+I-1), I=1,NLAMDS) 426 CALL GPCLOSE(LUSIFC,'KEEP') 427C 428 CALL CMO_REORDER(WORK(KCMO),WORK(KEND1),LWRK1) 429C 430C------------------------------------- 431C Add orbital relaxation terms. 432C------------------------------------- 433C 434c CALL CC_D1ORRE(AODEN,ZKAM(NMATIJ(1)+NMATAB(1)+1), 435c * WORK(KEND1),LWRK1) 436C 437 KOFDIJ = 1 438 KOFDAB = KOFDIJ + NMATIJ(1) 439 KOFDAI = KOFDAB + NMATAB(1) 440 KOFDIA = KOFDAI + NT1AMX 441C 442 ISDEN = 1 443 CALL CC_DENAO(AODEN,1,ZKAM(KOFDAI),ZKAM(KOFDAB), 444 * ZKAM(KOFDIJ),ZKAM(KOFDIA),ISDEN,WORK(KCMO),1, 445 * WORK(KCMO),1,WORK(KEND1),LWRK1) 446C 447C 448C------------------------------------------------------- 449C Add frozen core contributions to AO density. 450C------------------------------------------------------- 451C 452 IF (FROIMP) THEN 453C 454 KOFFAI = NMATIJ(1) + NMATAB(1) + 2*NT1AMX + 1 455 KOFFIA = KOFFAI + NT1FRO(1) 456 KOFFIJ = KOFFIA + NT1FRO(1) 457 KOFFJI = KOFFIJ + NCOFRO(1) 458C 459 ISDEN = 1 460 ICON = 2 461 CALL CC_D1FCB(AODEN,ZKAM(KOFFIJ),ZKAM(KOFFJI), 462 * ZKAM(KOFFAI),ZKAM(KOFFIA),WORK(KEND1), 463 * LWRK1,ISDEN,ICON) 464C 465 ENDIF 466 ENDIF 467C 468C 469 ELSE IF ((RCCD).or.(DRCCD)) THEN 470C 471C------------------------------------------------------- 472C Both zeta- and t-vectors are totally symmetric. 473C------------------------------------------------------- 474C 475 ISYMTR = 1 476 ISYMOP = 1 477C 478C-------------------------------------- 479C Initial work space allocation. 480C-------------------------------------- 481C 482 NHELP = NT2AMX 483C 484 KZ2AM = 1 485 KT2AM = KZ2AM + NT2SQ(1) 486 KT2AMT = KT2AM + NT2AMX 487 KLAMDP = KT2AMT + NHELP 488 KLAMDH = KLAMDP + NLAMDT 489 KT1AM = KLAMDH + NLAMDT 490 KZ1AM = KT1AM + NT1AMX 491 KEND1 = KZ1AM + NT1AMX 492 LWRK1 = LWORK - KEND1 493C 494 IF (LWRK1 .LT. 0) THEN 495 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 496 CALL QUIT('Insufficient core for initial '// 497 & 'allocation in CC_D1AO') 498 ENDIF 499C 500C------------------------------------------- 501C Read zero'th order zeta amplitudes. 502C------------------------------------------- 503C 504 ISYML = 1 505 IOPT = 2 506 CALL CC_RDRSP(LLIST,ILLNR,ISYML,IOPT,MODDUM, 507 & WORK(KZ1AM),WORK(KT2AM)) 508 IF ( IPRINT .GT. 10 ) THEN 509 RHO1N = DDOT(NT1AM(ISYML),WORK(KZ1AM),1,WORK(KZ1AM),1) 510 RHO2N = DDOT(NT2AM(ISYML),WORK(KT2AM),1,WORK(KT2AM),1) 511 WRITE(LUPRI,*) 'CC_D1AO: Norm of L. vector :',RHO1N 512 WRITE(LUPRI,*) 'CC_D1AO: Norm of L. vector :',RHO2N 513 ENDIF 514 515C 516C----------------------------------- 517C Square up zeta2 amplitudes. 518C----------------------------------- 519C 520 CALL CC_T2SQ(WORK(KT2AM),WORK(KZ2AM),1) 521C 522C---------------------------------------------- 523C Read zero'th order cluster amplitudes. 524C---------------------------------------------- 525C 526 IOPT = 2 527 CALL CC_RDRSP('R0',0,1,IOPT,MODEL,WORK(KT1AM),WORK(KT2AM)) 528C 529C--------------------------------------------------- 530C Zero out single vectors (do no exist) 531C--------------------------------------------------- 532C 533 CALL DZERO(WORK(KT1AM),NT1AMX) 534 CALL DZERO(WORK(KZ1AM),NT1AMX) 535C 536C------------------------------------------------- 537C Calculate the lambda matrices. 538C FIXME Useless (since T1=0) but I am lazy. Sonia 539C------------------------------------------------- 540C 541 CALL LAMMAT(WORK(KLAMDP),WORK(KLAMDH),WORK(KT1AM),WORK(KEND1), 542 * LWRK1) 543C 544C------------------------------------------ 545C Set up 2C-E of cluster amplitudes. 546C------------------------------------------ 547C 548 CALL DCOPY(NT2AMX,WORK(KT2AM),1,WORK(KT2AMT),1) 549 ISYOPE = 1 550 IOPTTCME = 1 551 CALL CCSD_TCMEPK(WORK(KT2AMT),1.0D0,ISYOPE,IOPTTCME) 552C 553C---------------------------------- 554C Work space allocation one. 555C Note that D(ai)=ZETA(ai) & 556C D(ia) is stored transposed 557C---------------------------------- 558C 559 KONEAI = KZ1AM 560 KONEAB = KONEAI + NT1AMX 561 KONEIJ = KONEAB + NMATAB(1) 562 KONEIA = KONEIJ + NMATIJ(1) 563 KXMAT = KONEIA + NT1AMX 564 KYMAT = KXMAT + NMATIJ(1) 565 KCMO = KYMAT + NMATAB(1) 566 KEND1 = KCMO + NLAMDS 567 LWRK1 = LWORK - KEND1 568C 569 IF (LWRK1 .LT. 0) THEN 570 WRITE(LUPRI,*) 'CC_D1AO; Available:',LWORK, 'Needed:',KEND1 571 CALL QUIT( 572 & 'Insufficient memory for first allocation in CC_D1AO') 573 ENDIF 574C 575C--------------------------------------------------------- 576C Initialize remaining one electron density arrays. 577C--------------------------------------------------------- 578C 579 CALL DZERO(WORK(KONEAB),NMATAB(1)) 580 CALL DZERO(WORK(KONEIJ),NMATIJ(1)) 581 CALL DZERO(WORK(KONEIA),NT1AMX) 582 CALL DZERO(WORK(KONEAI),NT1AMX) 583C 584C----------------------------------------------------------- 585C Calculate X-intermediate of zeta- and t-amplitudes. 586C----------------------------------------------------------- 587C 588 CALL CC_XI(WORK(KXMAT),WORK(KZ2AM),ISYMTR,WORK(KT2AM),ISYMOP, 589 * WORK(KEND1),LWRK1) 590 CALL DSCAL(NMATIJ(1),TWO,WORK(KXMAT),1) 591C 592C----------------------------------------------------------- 593C Calculate Y-intermediate of zeta- and t-amplitudes. 594C----------------------------------------------------------- 595C 596 CALL CC_YI(WORK(KYMAT),WORK(KZ2AM),ISYMTR,WORK(KT2AM),ISYMOP, 597 * WORK(KEND1),LWRK1) 598 CALL DSCAL(NMATAB(1),TWO,WORK(KYMAT),1) 599C 600C----------------------------------------------------------------- 601C Construct nonzero blocks of one electron density. 602C----------------------------------------------------------------- 603C 604 CALL DCOPY(NMATAB(1),WORK(KYMAT),1,WORK(KONEAB),1) 605 CALL CC_EITR(WORK(KONEAB),WORK(KONEIJ),WORK(KEND1),LWRK1,1) 606 !pure Hartree-Fock contribution to the density is added inside DIJGEN 607 CALL DIJGEN(WORK(KONEIJ),WORK(KXMAT)) 608C 609C----------------------------------------------------------- 610C Backtransform the one electron density to AO-basis. 611C----------------------------------------------------------- 612C 613 CALL DZERO(AODEN,N2BST(1)) 614C 615 ISDEN = 1 616 CALL CC_DENAO(AODEN,1,WORK(KONEAI),WORK(KONEAB), 617 * WORK(KONEIJ),WORK(KONEIA),ISDEN,WORK(KLAMDP),1, 618 * WORK(KLAMDH),1,WORK(KEND1),LWRK1) 619C 620 IF (RELORB) THEN 621C 622C---------------------------------------------------------------- 623C Read MO-coefficients from interface file and reorder. 624C---------------------------------------------------------------- 625C 626 LUSIFC = -1 627 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY, 628 & .FALSE.) 629 REWIND LUSIFC 630 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 631 READ (LUSIFC) 632 READ (LUSIFC) 633 READ (LUSIFC) (WORK(KCMO+I-1), I=1,NLAMDS) 634 CALL GPCLOSE(LUSIFC,'KEEP') 635C 636 CALL CMO_REORDER(WORK(KCMO),WORK(KEND1),LWRK1) 637C 638C------------------------------------- 639C Add orbital relaxation terms. 640C (backtransform kappabar_pq passed from outside) 641C------------------------------------- 642C 643c CALL CC_D1ORRE(AODEN,ZKAM(NMATIJ(1)+NMATAB(1)+1), 644c * WORK(KEND1),LWRK1) 645 KOFDIJ = 1 646 KOFDAB = KOFDIJ + NMATIJ(1) 647 KOFDAI = KOFDAB + NMATAB(1) 648 KOFDIA = KOFDAI + NT1AMX 649C 650! write(lupri,*)'The IJ kappabar of ', MODEL(1:5) 651! CALL OUTPUT(ZKAM(KOFDIJ),1,NRHF(1),1,NRHF(1), 652! * NRHF(1),NRHF(1),1,LUPRI) 653! write(lupri,*)'The AI kappabar of ', MODEL(1:5) 654! CALL OUTPUT(ZKAM(KOFDAI),1,NVIR(1),1,NRHF(1), 655! * NVIR(1),NRHF(1),1,LUPRI) 656! 657!! CALL DSCAL(NT1AMX,-ONE,ZKAM(KOFDIA),1) 658! write(lupri,*)'The IA kappabar of ', MODEL(1:5) 659! CALL OUTPUT(ZKAM(KOFDIA),1,NVIR(1),1,NRHF(1), 660! * NVIR(1),NRHF(1),1,LUPRI) 661! 662! write(lupri,*)'The AB kappabar of ', MODEL(1:5) 663! CALL OUTPUT(ZKAM(KOFDAB),1,NVIR(1),1,NVIR(1), 664! * NVIR(1),NVIR(1),1,LUPRI) 665 666 ISDEN = 1 667 CALL CC_DENAO(AODEN,1,ZKAM(KOFDAI),ZKAM(KOFDAB), 668 * ZKAM(KOFDIJ),ZKAM(KOFDIA),ISDEN,WORK(KCMO),1, 669 * WORK(KCMO),1,WORK(KEND1),LWRK1) 670C 671C------------------------------------------------------- 672C Add frozen core contributions to AO density. 673C------------------------------------------------------- 674C 675! IF (FROIMP) THEN 676C 677! KOFFAI = NMATIJ(1) + NMATAB(1) + 2*NT1AMX + 1 678! KOFFIA = KOFFAI + NT1FRO(1) 679! KOFFIJ = KOFFIA + NT1FRO(1) 680! KOFFJI = KOFFIJ + NCOFRO(1) 681C 682! ISDEN = 1 683! ICON = 2 684! CALL CC_D1FCB(AODEN,ZKAM(KOFFIJ),ZKAM(KOFFJI), 685! * ZKAM(KOFFAI),ZKAM(KOFFIA),WORK(KEND1), 686! * LWRK1,ISDEN,ICON) 687C 688! ENDIF 689 ENDIF 690C 691 ENDIF 692C 693 IF (FROIMP) THEN 694 CALL QEXIT('CC_D1AO') 695 RETURN 696 END IF 697C 698C----------------------------------------------------------------- 699C Calculate additional terms to the one electron density 700C needed for the evaluation of the natural occupation numbers. 701C Only if NO, else skip. 702C----------------------------------------------------------------- 703C 704 IF (.NOT. NO) THEN 705 CALL QEXIT('CC_D1AO') 706 RETURN 707 END IF 708C 709 IF (.NOT. MP2) THEN 710C 711 CALL CCD1T1CO(WORK(KONEAI),WORK(KONEAB),WORK(KONEIJ), 712 * WORK(KONEIA),WORK(KT1AM),WORK(KYMAT), 713 * WORK(KEND1),LWRK1) 714C 715 ENDIF 716C 717C-------------------------------------------------- 718C Reorder and diagonalize one electron density. 719C-------------------------------------------------- 720C 721 KDENFO = KEND1 722 KEND2 = KDENFO + N2BST(1) 723 LWRK2 = LWORK - KEND2 724C 725 IF (LWRK2 .LT. 0) THEN 726 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND2 727 CALL QUIT('Insufficient memory for second '// 728 & 'allocation in CC_D1AO') 729 ENDIF 730C 731 CALL DZERO(WORK(KDENFO),N2BST(1)) 732C 733 CALL CCD1REOR(WORK(KDENFO),WORK(KONEAI),WORK(KONEAB), 734 * WORK(KONEIJ),WORK(KONEIA)) 735C 736 IF (IPRINT .GT. 10) THEN 737C 738 CALL AROUND('Nondiagonalised CC one electron density') 739C 740 DO 100 ISYM = 1,NSYM 741C 742 WRITE(LUPRI,333) 'Symmetry block number:', ISYM 743 333 FORMAT(3X,A22,2X,I1) 744C 745 KOFFAI = KONEAI + IT1AM(ISYM,ISYM) 746 KOFFAB = KONEAB + IMATAB(ISYM,ISYM) 747 KOFFIJ = KONEIJ + IMATIJ(ISYM,ISYM) 748 KOFFIA = KONEIA + IT1AM(ISYM,ISYM) 749C 750 CALL AROUND('DAI') 751 CALL OUTPUT(WORK(KOFFAI),1,NVIR(ISYM),1,NRHF(ISYM), 752 * NVIR(ISYM),NRHF(ISYM),1,LUPRI) 753C 754 CALL AROUND('DAB') 755 CALL OUTPUT(WORK(KOFFAB),1,NVIR(ISYM),1,NVIR(ISYM), 756 * NVIR(ISYM),NVIR(ISYM),1,LUPRI) 757C 758 CALL AROUND('DIJ') 759 CALL OUTPUT(WORK(KOFFIJ),1,NRHF(ISYM),1,NRHF(ISYM), 760 * NRHF(ISYM),NRHF(ISYM),1,LUPRI) 761C 762 CALL AROUND('DIA') 763 CALL OUTPUT(WORK(KOFFIA),1,NVIR(ISYM),1,NRHF(ISYM), 764 * NVIR(ISYM),NRHF(ISYM),1,LUPRI) 765C 766 100 CONTINUE 767C 768 ENDIF 769C 770 CALL CC_DEDIAN(WORK(KDENFO),MODEL,WORK(KEND2),LWRK2) 771C 772C----------------------- 773C write out timings. 774C----------------------- 775C 776 TIMTOT = SECOND() -TIMTOT 777C 778 IF (IPRINT .GT. 9) THEN 779 WRITE(LUPRI,*) ' ' 780 WRITE(LUPRI,*) ' ' 781 WRITE(LUPRI,*) 'CC one electron density '// 782 & 'calculated & diagonalized' 783 WRITE(LUPRI,*) 'Total time used in CC_D1AO:', TIMTOT 784 ENDIF 785C 786 CALL QEXIT('CC_D1AO') 787 RETURN 788 END 789C /* Deck cc_denao */ 790 SUBROUTINE CC_DENAO(AODEN,ISYDAO,DENAI,DENAB,DENIJ,DENIA,ISYDMO, 791 * XLAMDP,ISYMLP,XLAMDH,ISYMLH,WORK,LWORK) 792C 793C Written by Asger Halkier 26/1 - 1996 794C 795C Version: 1.0 796C 797C Symmetries: ISYDAO -- AODEN 798C ISYDMO -- DENAI,DENAB,DENIJ,DENIA 799C ISYMLP -- XLAMDP 800C ISYMLH -- XLAMDH 801C 802C Purpose: To transform all four blocks of the Coupled 803C Cluster one electron density to AO-basis and add 804C the results! 805C 806#include "implicit.h" 807 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 808 DIMENSION AODEN(*), DENAI(*), DENAB(*), DENIJ(*), DENIA(*) 809 DIMENSION XLAMDP(*), XLAMDH(*), WORK(LWORK) 810#include "priunit.h" 811#include "ccorb.h" 812#include "ccsdsym.h" 813C 814 CALL QENTER('CC_DENAO') 815C 816 DO ISYM1 = 1,NSYM 817C 818 ISYM2 = MULD2H(ISYDMO,ISYM1) 819 ISYMA = MULD2H(ISYMLP,ISYM1) 820 ISYMB = MULD2H(ISYMLH,ISYM2) 821C 822 IF (ISYDAO.NE.MULD2H(ISYMA,ISYMB)) THEN 823 CALL QUIT('Symmetry mismatch in CC_DENAO') 824 END IF 825C 826 LVIR = MAX(NBAS(ISYMA)*NVIR(ISYM2),NBAS(ISYMB)*NVIR(ISYM1)) 827 LRHF = MAX(NBAS(ISYMA)*NRHF(ISYM2),NBAS(ISYMB)*NRHF(ISYM2)) 828 LNEED = MAX(LVIR,LRHF) 829C 830 IF (LWORK .LT. LNEED) THEN 831 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', LNEED 832 CALL QUIT('Insufficient work space in CC_DENAO') 833 ENDIF 834C 835 CALL DZERO(WORK,LNEED) 836C 837 NBASA = MAX(NBAS(ISYMA),1) 838 NBASB = MAX(NBAS(ISYMB),1) 839 NTOVI1 = MAX(NVIR(ISYM1),1) 840 NTOVI2 = MAX(NVIR(ISYM2),1) 841 NTORH1 = MAX(NRHF(ISYM1),1) 842C 843C----------------------------------------- 844C Transform virtual-occupied block. 845C----------------------------------------- 846C 847 KOFF1 = IGLMVI(ISYMA,ISYM1) + 1 848 KOFF2 = IT1AM(ISYM1,ISYM2) + 1 849C 850 CALL DGEMM('N','N',NBAS(ISYMA),NRHF(ISYM2),NVIR(ISYM1), 851 * ONE,XLAMDP(KOFF1),NBASA,DENAI(KOFF2),NTOVI1, 852 * ZERO,WORK,NBASA) 853C 854 KOFF3 = IGLMRH(ISYMB,ISYM2) + 1 855 KOFF4 = IAODIS(ISYMA,ISYMB) + 1 856C 857 CALL DGEMM('N','T',NBAS(ISYMA),NBAS(ISYMB),NRHF(ISYM2), 858 * ONE,WORK,NBASA,XLAMDH(KOFF3),NBASB,ONE, 859 * AODEN(KOFF4),NBASA) 860C 861C---------------------------------------- 862C Transform virtual-virtual block. 863C---------------------------------------- 864C 865 KOFF1 = IGLMVI(ISYMA,ISYM1) + 1 866 KOFF2 = IMATAB(ISYM1,ISYM2) + 1 867C 868 CALL DGEMM('N','N',NBAS(ISYMA),NVIR(ISYM2),NVIR(ISYM1), 869 * ONE,XLAMDP(KOFF1),NBASA,DENAB(KOFF2),NTOVI1, 870 * ZERO,WORK,NBASA) 871C 872 KOFF3 = IGLMVI(ISYMB,ISYM2) + 1 873 KOFF4 = IAODIS(ISYMA,ISYMB) + 1 874C 875 CALL DGEMM('N','T',NBAS(ISYMA),NBAS(ISYMB),NVIR(ISYM2), 876 * ONE,WORK,NBASA,XLAMDH(KOFF3),NBASB,ONE, 877 * AODEN(KOFF4),NBASA) 878C 879C------------------------------------------ 880C Transform occupied-occupied block. 881C------------------------------------------ 882C 883 KOFF1 = IGLMRH(ISYMA,ISYM1) + 1 884 KOFF2 = IMATIJ(ISYM1,ISYM2) + 1 885C 886C 887 CALL DGEMM('N','N',NBAS(ISYMA),NRHF(ISYM2),NRHF(ISYM1), 888 * ONE,XLAMDP(KOFF1),NBASA,DENIJ(KOFF2),NTORH1, 889 * ZERO,WORK,NBASA) 890C 891 KOFF3 = IGLMRH(ISYMB,ISYM2) + 1 892 KOFF4 = IAODIS(ISYMA,ISYMB) + 1 893C 894 CALL DGEMM('N','T',NBAS(ISYMA),NBAS(ISYMB),NRHF(ISYM2), 895 * ONE,WORK,NBASA,XLAMDH(KOFF3),NBASB,ONE, 896 * AODEN(KOFF4),NBASA) 897C 898C------------------------------------------------------------- 899C Transform occupied-virtual block (stored transposed). 900C------------------------------------------------------------- 901C 902 KOFF1 = IGLMVI(ISYMB,ISYM2) + 1 903 KOFF2 = IT1AM(ISYM2,ISYM1) + 1 904C 905 CALL DGEMM('N','N',NBAS(ISYMB),NRHF(ISYM1),NVIR(ISYM2), 906 * ONE,XLAMDH(KOFF1),NBASB,DENIA(KOFF2),NTOVI2, 907 * ZERO,WORK,NBASB) 908C 909 KOFF3 = IGLMRH(ISYMA,ISYM1) + 1 910 KOFF4 = IAODIS(ISYMA,ISYMB) + 1 911C 912 CALL DGEMM('N','T',NBAS(ISYMA),NBAS(ISYMB),NRHF(ISYM1), 913 * ONE,XLAMDP(KOFF3),NBASA,WORK,NBASB,ONE, 914 * AODEN(KOFF4),NBASA) 915C 916 END DO 917C 918 CALL QEXIT('CC_DENAO') 919C 920 RETURN 921 END 922C /* Deck dijgen */ 923 SUBROUTINE DIJGEN(DIJ,XMAT) 924C 925C Written by Asger Halkier 26/1 - 1996 926C 927C Version: 1.0 928C 929C Purpose: To set up the (occ,occ) part of the Coupled 930C Cluster one electron density. 931C 932#include "implicit.h" 933 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 934 DIMENSION DIJ(*), XMAT(*) 935#include "priunit.h" 936#include "ccorb.h" 937#include "ccsdsym.h" 938C 939C-------------------------- 940C Set up density block. 941C-------------------------- 942C 943 CALL QENTER('DIJGEN') 944C 945 CALL DAXPY(NMATIJ(1),-ONE,XMAT,1,DIJ,1) 946C 947 DO 100 ISYMI = 1,NSYM 948C 949 ISYMJ = ISYMI 950C 951 DO 110 I = 1,NRHF(ISYMI) 952C 953 NIJ = IMATIJ(ISYMI,ISYMJ) + NRHF(ISYMI)*(I - 1) + I 954C 955 DIJ(NIJ) = DIJ(NIJ) + TWO 956C 957 110 CONTINUE 958C 959 100 CONTINUE 960C 961 CALL QEXIT('DIJGEN') 962C 963 RETURN 964 END 965C /* Deck diagen */ 966 SUBROUTINE DIAGEN(DIA,T2AMT,Z1AM) 967C 968C Written by Asger Halkier 26/1 - 1996 969C 970C Version: 1.0 971C 972C Purpose: To set up the (occ,vir) part of the Coupled 973C Cluster one electron density. 974C 975C NOTE: This block is stored transposed, i.e. like a t1-amplitude! 976C 977#include "implicit.h" 978 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 979 DIMENSION DIA(*), T2AMT(*), Z1AM(*) 980#include "priunit.h" 981#include "ccorb.h" 982#include "ccsdsym.h" 983C 984 INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 985C 986 CALL QENTER('DIAGEN') 987C 988 ISYMDK = 1 989 ISYMAI = ISYMDK 990C 991C-------------------------- 992C Set up density block. 993C-------------------------- 994C 995 DO 100 NAI = 1,NT1AM(ISYMAI) 996C 997 DO 110 NDK = 1,NT1AM(ISYMDK) 998C 999 NDKAI = IT2AM(ISYMDK,ISYMAI) + INDEX(NDK,NAI) 1000C 1001 DIA(NAI) = DIA(NAI) + T2AMT(NDKAI)*Z1AM(NDK) 1002C 1003 110 CONTINUE 1004 100 CONTINUE 1005C 1006 CALL QEXIT('DIAGEN') 1007C 1008 RETURN 1009 END 1010C /* Deck dijk01 */ 1011 SUBROUTINE DIJK01(D2IJK,ISYIJK,XLAMDH,ISYMLH,IDEL,ISYMD) 1012C 1013C Written by Asger Halkier 29/1 - 1996 1014C 1015C Version: 1.0 1016C 1017C Purpose: To calculate the first contribution to D2IJK! 1018C 1019#include "implicit.h" 1020 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, FOUR = 4.0D0) 1021 DIMENSION D2IJK(*), XLAMDH(*) 1022#include "priunit.h" 1023#include "ccorb.h" 1024#include "ccsdsym.h" 1025C 1026 CALL QENTER('DIJK01') 1027C 1028 ISYMK = MULD2H(ISYMD,ISYMLH) 1029 ISYMIJ = MULD2H(ISYIJK,ISYMK) 1030C 1031 IF (ISYMIJ.NE.1) CALL QUIT('Symmetry mismatch in DIJK01') 1032C 1033C----------------------------------------------- 1034C Calculate special adresses and add result. 1035C----------------------------------------------- 1036C 1037 DO 100 K = 1,NRHF(ISYMK) 1038C 1039 DO 110 ISYMJ = 1,NSYM 1040C 1041 DO 120 J = 1,NRHF(ISYMJ) 1042C 1043 NIJK = IMAIJK(ISYMIJ,ISYMK) + NMATIJ(ISYMIJ)*(K - 1) 1044 * + IMATIJ(ISYMJ,ISYMJ) + NRHF(ISYMJ)*(J - 1) + J 1045 NDEL = IGLMRH(ISYMD,ISYMK) + NBAS(ISYMD)*(K - 1) 1046 * + IDEL - IBAS(ISYMD) 1047C 1048 D2IJK(NIJK) = D2IJK(NIJK) + FOUR*XLAMDH(NDEL) 1049C 1050 120 CONTINUE 1051 110 CONTINUE 1052 100 CONTINUE 1053C 1054 CALL QEXIT('DIJK01') 1055C 1056 RETURN 1057 END 1058C /* Deck dijk02 */ 1059 SUBROUTINE DIJK02(D2IJK,ISYIJK,XLAMDH,ISYMLH,IDEL,ISYMD) 1060C 1061C Written by Asger Halkier 29/1 - 1996 1062C 1063C Version: 1.0 1064C 1065C Purpose: To calculate the second contribution to D2IJK! 1066C 1067#include "implicit.h" 1068 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 1069 DIMENSION D2IJK(*), XLAMDH(*) 1070#include "priunit.h" 1071#include "ccorb.h" 1072#include "ccsdsym.h" 1073C 1074 CALL QENTER('DIJK02') 1075C 1076 ISYMI = MULD2H(ISYMD,ISYMLH) 1077C 1078 IF (ISYMI.NE.ISYIJK) CALL QUIT('Symmetry mismatch in DIJK02') 1079C 1080C----------------------------------------------- 1081C Calculate special adresses and add result. 1082C----------------------------------------------- 1083C 1084 DO 100 ISYMK = 1,NSYM 1085C 1086 ISYMIJ = MULD2H(ISYMI,ISYMK) 1087C 1088 DO 110 K = 1,NRHF(ISYMK) 1089C 1090 DO 120 I = 1,NRHF(ISYMI) 1091C 1092 NIJK = IMAIJK(ISYMIJ,ISYMK) + NMATIJ(ISYMIJ)*(K - 1) 1093 * + IMATIJ(ISYMI,ISYMK) + NRHF(ISYMI)*(K - 1) + I 1094 NDEL = IGLMRH(ISYMD,ISYMI) + NBAS(ISYMD)*(I - 1) 1095 * + IDEL - IBAS(ISYMD) 1096C 1097 D2IJK(NIJK) = D2IJK(NIJK) - TWO*XLAMDH(NDEL) 1098C 1099 120 CONTINUE 1100 110 CONTINUE 1101 100 CONTINUE 1102C 1103 CALL QEXIT('DIJK02') 1104C 1105 RETURN 1106 END 1107C /* Deck cc_d2pao */ 1108 SUBROUTINE CC_D2PAO(D2AOIJ,D2AOAI,D2AOAB,D2AOIA,XLAMDP,D2IJK, 1109 * D2AIJ,D2IAJ,D2ABI,G,ISYMG,ISYDEN) 1110C 1111C Written by Asger Halkier 30/1 - 1996 1112C 1113C Version: 1.0 1114C 1115C Purpose: To backtransform the four 2 electron densities 1116C d(pq,k;del) to AO basis for a given gamma d(pq;gam;del)! 1117C 1118#include "implicit.h" 1119 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 1120 DIMENSION D2AOIJ(*), D2AOAI(*), D2AOAB(*), D2AOIA(*), XLAMDP(*) 1121 DIMENSION D2IJK(*), D2AIJ(*), D2IAJ(*), D2ABI(*) 1122#include "priunit.h" 1123#include "ccorb.h" 1124#include "ccsdsym.h" 1125C 1126 CALL QENTER('CC_D2PAO') 1127C 1128 ISYMK = ISYMG 1129 ISYMPQ = MULD2H(ISYMG,ISYDEN) 1130C 1131 KOFFGK = IGLMRH(ISYMG,ISYMK) + G 1132C 1133C------------------------------- 1134C Transform (occ,occ) block. 1135C------------------------------- 1136C 1137 KOFF1 = IMAIJK(ISYMPQ,ISYMK) + 1 1138C 1139 NTOTIJ = MAX(NMATIJ(ISYMPQ),1) 1140C 1141 CALL DGEMV('N',NMATIJ(ISYMPQ),NRHF(ISYMK),ONE,D2IJK(KOFF1), 1142 * NTOTIJ,XLAMDP(KOFFGK),NBAS(ISYMG),ZERO,D2AOIJ,1) 1143C 1144C---------------------------------------------- 1145C Transform (vir,occ) and (occ,vir) blocks. 1146C---------------------------------------------- 1147C 1148 KOFF2 = IT2BCD(ISYMPQ,ISYMK) + 1 1149C 1150 NTOTAI = MAX(NT1AM(ISYMPQ),1) 1151C 1152 CALL DGEMV('N',NT1AM(ISYMPQ),NRHF(ISYMK),ONE,D2AIJ(KOFF2), 1153 * NTOTAI,XLAMDP(KOFFGK),NBAS(ISYMG),ZERO,D2AOAI,1) 1154C 1155 CALL DGEMV('N',NT1AM(ISYMPQ),NRHF(ISYMK),ONE,D2IAJ(KOFF2), 1156 * NTOTAI,XLAMDP(KOFFGK),NBAS(ISYMG),ZERO,D2AOIA,1) 1157C 1158C------------------------------- 1159C Transform (vir,vir) block. 1160C------------------------------- 1161C 1162 KOFF3 = ICKASR(ISYMPQ,ISYMK) + 1 1163C 1164 NTOTAB = MAX(NMATAB(ISYMPQ),1) 1165C 1166 CALL DGEMV('N',NMATAB(ISYMPQ),NRHF(ISYMK),ONE,D2ABI(KOFF3), 1167 * NTOTAB,XLAMDP(KOFFGK),NBAS(ISYMG),ZERO,D2AOAB,1) 1168C 1169 CALL QEXIT('CC_D2PAO') 1170C 1171 RETURN 1172 END 1173C /* Deck ccsd_d2en */ 1174 SUBROUTINE CCSD_D2EN(ECCSD,XINT,XLAMDH,XLAMDP,D2AOIJ,D2AOAI, 1175 * D2AOAB,D2AOIA,G,ISYMG,ISYMD,WORK,LWORK) 1176C 1177C Written by Asger Halkier 30/1 - 1996 1178C 1179C Version: 1.0 1180C 1181C Purpose: To calculate the contribution from the 2 electron 1182C Coupled Cluster density to the ccsd energy! 1183C 1184#include "implicit.h" 1185 PARAMETER(ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0) 1186 DIMENSION ECCSD(*), XINT(*), XLAMDH(*), XLAMDP(*) 1187 DIMENSION D2AOIJ(*), D2AOAI(*), D2AOAB(*), D2AOIA(*), WORK(LWORK) 1188#include "priunit.h" 1189#include "ccorb.h" 1190#include "ccsdsym.h" 1191C 1192 CALL QENTER('CCSD_D2EN') 1193C 1194 ISYINT = ISYMD 1195 ISYDEN = ISYMD 1196 ISYMPQ = MULD2H(ISYMG,ISYDEN) 1197 ISALBE = ISYMPQ 1198C 1199C------------------------------- 1200C Work space allocation one. 1201C------------------------------- 1202C 1203 KAOINT = 1 1204 KEND1 = KAOINT + N2BST(ISALBE) 1205 LWRK1 = LWORK - KEND1 1206C 1207 IF (LWRK1. LT. 0) THEN 1208 WRITE(LUPRI,*) 'Needed:', KEND1, 'Available:', LWORK 1209 CALL QUIT('Insufficient memory for allocation in CCSD_D2EN') 1210 ENDIF 1211C 1212C------------------------------------- 1213C Square up integral distribution. 1214C------------------------------------- 1215C 1216 KOFF1 = IDSAOG(ISYMG,ISYINT) + NNBST(ISALBE)*(G - 1) + 1 1217C 1218 CALL CCSD_SYMSQ(XINT(KOFF1),ISALBE,WORK(KAOINT)) 1219C 1220C------------------------------- 1221C Work space allocation two. 1222C------------------------------- 1223C 1224 LINTMO = MAX(NT1AM(ISYMPQ),NMATIJ(ISYMPQ),NMATAB(ISYMPQ)) 1225C 1226 KINTMO = KEND1 1227 KEND2 = KINTMO + LINTMO 1228 LWRK2 = LWORK - KEND2 1229C 1230 IF (LWRK2. LT. 0) THEN 1231 WRITE(LUPRI,*) 'Needed:', KEND2, 'Available:', LWORK 1232 CALL QUIT('Insufficient memory for allocation in CCSD_D2EN') 1233 ENDIF 1234C 1235C-------------------------------------- 1236C Calculate (occ,occ) contribution. 1237C-------------------------------------- 1238C 1239 CALL DZERO(WORK(KINTMO),NMATIJ(ISYMPQ)) 1240C 1241 DO 100 ISYMAL = 1,NSYM 1242C 1243 ISYMBE = MULD2H(ISYMAL,ISALBE) 1244 ISYMJ = ISYMBE 1245 ISYMI = ISYMAL 1246C 1247 LNEED = LWORK - KEND2 - NBAS(ISYMAL)*NRHF(ISYMJ) 1248C 1249 IF (LNEED .LT. 0) THEN 1250 WRITE(LUPRI,*) 'Work exceeded by:', -LNEED 1251 CALL QUIT('Insufficient work space in CCSD_D2EN') 1252 ENDIF 1253C 1254 CALL DZERO(WORK(KEND2),NBAS(ISYMAL)*NRHF(ISYMJ)) 1255C 1256 KOFF2 = KAOINT + IAODIS(ISYMAL,ISYMBE) 1257 KOFF3 = IGLMRH(ISYMBE,ISYMJ) + 1 1258C 1259 NTOTAL = MAX(NBAS(ISYMAL),1) 1260 NTOTBE = MAX(NBAS(ISYMBE),1) 1261C 1262 CALL DGEMM('N','N',NBAS(ISYMAL),NRHF(ISYMJ),NBAS(ISYMBE),ONE, 1263 * WORK(KOFF2),NTOTAL,XLAMDH(KOFF3),NTOTBE,ZERO, 1264 * WORK(KEND2),NTOTAL) 1265C 1266 KOFF4 = IGLMRH(ISYMAL,ISYMI) + 1 1267 KOFF5 = KINTMO + IMATIJ(ISYMI,ISYMJ) 1268C 1269 NTOTAL = MAX(NBAS(ISYMAL),1) 1270 NTOTI = MAX(NRHF(ISYMI),1) 1271C 1272 CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),NBAS(ISYMAL),HALF, 1273 * XLAMDP(KOFF4),NTOTAL,WORK(KEND2),NTOTAL,ZERO, 1274 * WORK(KOFF5),NTOTI) 1275C 1276 100 CONTINUE 1277C 1278 ECCSD(1) = ECCSD(1) + DDOT(NMATIJ(ISYMPQ),WORK(KINTMO),1, 1279 * D2AOIJ,1) 1280C 1281C-------------------------------------- 1282C Calculate (vir,vir) contribution. 1283C-------------------------------------- 1284C 1285 CALL DZERO(WORK(KINTMO),NMATAB(ISYMPQ)) 1286C 1287 DO 110 ISYMAL = 1,NSYM 1288C 1289 ISYMBE = MULD2H(ISYMAL,ISALBE) 1290 ISYMA = ISYMAL 1291 ISYMB = ISYMBE 1292C 1293 LNEED = LWORK - KEND2 - NBAS(ISYMAL)*NVIR(ISYMB) 1294C 1295 IF (LNEED .LT. 0) THEN 1296 WRITE(LUPRI,*) 'Work exceeded by:', -LNEED 1297 CALL QUIT('Insufficient work space in CCSD_D2EN') 1298 ENDIF 1299C 1300 CALL DZERO(WORK(KEND2),NBAS(ISYMAL)*NVIR(ISYMB)) 1301C 1302 KOFF6 = KAOINT + IAODIS(ISYMAL,ISYMBE) 1303 KOFF7 = IGLMVI(ISYMBE,ISYMB) + 1 1304C 1305 NTOTAL = MAX(NBAS(ISYMAL),1) 1306 NTOTBE = MAX(NBAS(ISYMBE),1) 1307C 1308 CALL DGEMM('N','N',NBAS(ISYMAL),NVIR(ISYMB),NBAS(ISYMBE),ONE, 1309 * WORK(KOFF6),NTOTAL,XLAMDH(KOFF7),NTOTBE,ZERO, 1310 * WORK(KEND2),NTOTAL) 1311C 1312 KOFF8 = IGLMVI(ISYMAL,ISYMA) + 1 1313 KOFF9 = KINTMO + IMATAB(ISYMA,ISYMB) 1314C 1315 NTOTAL = MAX(NBAS(ISYMAL),1) 1316 NTOTA = MAX(NVIR(ISYMA),1) 1317C 1318 CALL DGEMM('T','N',NVIR(ISYMA),NVIR(ISYMB),NBAS(ISYMAL),HALF, 1319 * XLAMDP(KOFF8),NTOTAL,WORK(KEND2),NTOTAL,ZERO, 1320 * WORK(KOFF9),NTOTA) 1321C 1322 110 CONTINUE 1323C 1324 ECCSD(1) = ECCSD(1) + DDOT(NMATAB(ISYMPQ),WORK(KINTMO),1, 1325 * D2AOAB,1) 1326C 1327C------------------------------------------ 1328C Calculate the (vir,occ) contribution. 1329C------------------------------------------ 1330C 1331 CALL DZERO(WORK(KINTMO),NT1AM(ISYMPQ)) 1332C 1333 DO 120 ISYMAL = 1,NSYM 1334C 1335 ISYMBE = MULD2H(ISYMAL,ISALBE) 1336 ISYMI = ISYMBE 1337 ISYMA = ISYMAL 1338C 1339 LNEED = LWORK - KEND2 - NBAS(ISYMAL)*NRHF(ISYMI) 1340C 1341 IF (LNEED .LT. 0) THEN 1342 WRITE(LUPRI,*) 'Work exceeded by:', -LNEED 1343 CALL QUIT('Insufficient work space in CCSD_D2EN') 1344 ENDIF 1345C 1346 CALL DZERO(WORK(KEND2),NBAS(ISYMAL)*NRHF(ISYMI)) 1347C 1348 KOFF10 = KAOINT + IAODIS(ISYMAL,ISYMBE) 1349 KOFF11 = IGLMRH(ISYMBE,ISYMI) + 1 1350C 1351 NTOTAL = MAX(NBAS(ISYMAL),1) 1352 NTOTBE = MAX(NBAS(ISYMBE),1) 1353C 1354 CALL DGEMM('N','N',NBAS(ISYMAL),NRHF(ISYMI),NBAS(ISYMBE),ONE, 1355 * WORK(KOFF10),NTOTAL,XLAMDH(KOFF11),NTOTBE,ZERO, 1356 * WORK(KEND2),NTOTAL) 1357C 1358 KOFF12 = IGLMVI(ISYMAL,ISYMA) + 1 1359 KOFF13 = KINTMO + IT1AM(ISYMA,ISYMI) 1360C 1361 NTOTAL = MAX(NBAS(ISYMAL),1) 1362 NTOTA = MAX(NVIR(ISYMA),1) 1363C 1364 CALL DGEMM('T','N',NVIR(ISYMA),NRHF(ISYMI),NBAS(ISYMAL),HALF, 1365 * XLAMDP(KOFF12),NTOTAL,WORK(KEND2),NTOTAL,ZERO, 1366 * WORK(KOFF13),NTOTA) 1367C 1368 120 CONTINUE 1369C 1370 ECCSD(1) = ECCSD(1) + DDOT(NT1AM(ISYMPQ),WORK(KINTMO),1, 1371 * D2AOAI,1) 1372C 1373C------------------------------------------ 1374C Calculate the (occ,vir) contribution. 1375C------------------------------------------ 1376C 1377 CALL DZERO(WORK(KINTMO),NT1AM(ISYMPQ)) 1378C 1379 DO 130 ISYMAL = 1,NSYM 1380C 1381 ISYMBE = MULD2H(ISYMAL,ISALBE) 1382 ISYMI = ISYMAL 1383 ISYMA = ISYMBE 1384C 1385 LNEED = LWORK - KEND2 - NBAS(ISYMBE)*NRHF(ISYMI) 1386C 1387 IF (LNEED .LT. 0) THEN 1388 WRITE(LUPRI,*) 'Work exceeded by:', -LNEED 1389 CALL QUIT('Insufficient work space in CCSD_D2EN') 1390 ENDIF 1391C 1392 CALL DZERO(WORK(KEND2),NBAS(ISYMBE)*NRHF(ISYMI)) 1393C 1394 KOFF14 = KAOINT + IAODIS(ISYMAL,ISYMBE) 1395 KOFF15 = IGLMRH(ISYMAL,ISYMI) + 1 1396C 1397 NTOTAL = MAX(NBAS(ISYMAL),1) 1398 NTOTBE = MAX(NBAS(ISYMBE),1) 1399C 1400 CALL DGEMM('T','N',NBAS(ISYMBE),NRHF(ISYMI),NBAS(ISYMAL),ONE, 1401 * WORK(KOFF14),NTOTAL,XLAMDP(KOFF15),NTOTAL,ZERO, 1402 * WORK(KEND2),NTOTBE) 1403C 1404 KOFF16 = IGLMVI(ISYMBE,ISYMA) + 1 1405 KOFF17 = KINTMO + IT1AM(ISYMA,ISYMI) 1406C 1407 NTOTBE = MAX(NBAS(ISYMBE),1) 1408 NTOTA = MAX(NVIR(ISYMA),1) 1409C 1410 CALL DGEMM('T','N',NVIR(ISYMA),NRHF(ISYMI),NBAS(ISYMBE),HALF, 1411 * XLAMDH(KOFF16),NTOTBE,WORK(KEND2),NTOTBE,ZERO, 1412 * WORK(KOFF17),NTOTA) 1413C 1414 130 CONTINUE 1415C 1416 ECCSD(1) = ECCSD(1) + DDOT(NT1AM(ISYMPQ),WORK(KINTMO),1, 1417 * D2AOIA,1) 1418C 1419 CALL QEXIT('CCSD_D2EN') 1420C 1421 RETURN 1422 END 1423C /* Deck dijk11 */ 1424 SUBROUTINE DIJK11(D2IJK,ISYIJK,XLAMDH,ISYMLH,DIA,ISYD1, 1425 * WORK,LWORK,IDEL,ISYMD) 1426C 1427C Written by Asger Halkier 2/2 - 1996 1428C 1429C Version: 1.0 1430C 1431C Purpose: To calculate the first two contribution to D2IJK 1432C from projection against singles! 1433C 1434#include "implicit.h" 1435 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 1436 DIMENSION D2IJK(*), XLAMDH(*), DIA(*), WORK(LWORK) 1437#include "priunit.h" 1438#include "ccorb.h" 1439#include "ccsdsym.h" 1440C 1441 CALL QENTER('DIJK11') 1442C 1443 ISYMA = MULD2H(ISYMD,ISYMLH) 1444 ISYMK = MULD2H(ISYMA,ISYD1) 1445 ISYMIJ = MULD2H(ISYMK,ISYIJK) 1446C 1447 IF (ISYMK.NE.ISYIJK) CALL QUIT('Symmetry mismatch in DIJK11 (1)') 1448C 1449 IF (LWORK .LT. NRHF(ISYMK)) THEN 1450 WRITE(LUPRI,*) 'Available:', LWORK, 'NEEDED:', NRHF(ISYMK) 1451 CALL QUIT('Insufficient space for allocation in DIJK11') 1452 ENDIF 1453C 1454 CALL DZERO(WORK,NRHF(ISYMK)) 1455C 1456C------------------------------------------- 1457C Contract D(ia) with lamda hole matrix. 1458C------------------------------------------- 1459C 1460 KOFF1 = IT1AM(ISYMA,ISYMK) + 1 1461 KOFF2 = IGLMVI(ISYMD,ISYMA) + IDEL - IBAS(ISYMD) 1462C 1463 NTOTA = MAX(NVIR(ISYMA),1) 1464C 1465 CALL DGEMV('T',NVIR(ISYMA),NRHF(ISYMK),ONE,DIA(KOFF1),NTOTA, 1466 * XLAMDH(KOFF2),NBAS(ISYMD),ZERO,WORK,1) 1467C 1468C--------------------------------------- 1469C Calculate the DIJK11 contribution. 1470C--------------------------------------- 1471C 1472 DO 100 K = 1,NRHF(ISYMK) 1473C 1474 DO 110 ISYMI = 1,NSYM 1475C 1476 ISYMJ = ISYMI 1477C 1478 DO 120 I = 1,NRHF(ISYMI) 1479C 1480 J = I 1481 NIJK = IMAIJK(ISYMIJ,ISYMK) + NMATIJ(ISYMIJ)*(K - 1) 1482 * + IMATIJ(ISYMI,ISYMJ) + NRHF(ISYMI)*(J - 1) + I 1483C 1484 D2IJK(NIJK) = D2IJK(NIJK) + TWO*WORK(K) 1485C 1486 120 CONTINUE 1487 110 CONTINUE 1488 100 CONTINUE 1489C 1490C--------------------------------------- 1491C Calculate the DIJK12 contribution. 1492C--------------------------------------- 1493C 1494 ISYMI = MULD2H(ISYMA,ISYD1) 1495 ISYMJK = MULD2H(ISYMI,ISYIJK) 1496 1497 IF (ISYMI.NE.ISYIJK) CALL QUIT('Symmetry mismatch in DIJK11 (2)') 1498C 1499 DO 130 ISYMK = 1,NSYM 1500C 1501 ISYMJ = ISYMK 1502 ISYMIJ = MULD2H(ISYMI,ISYMJ) 1503C 1504 DO 140 K = 1,NRHF(ISYMK) 1505C 1506 J = K 1507C 1508 DO 150 I = 1,NRHF(ISYMI) 1509C 1510 NIJK = IMAIJK(ISYMIJ,ISYMK) + NMATIJ(ISYMIJ)*(K - 1) 1511 * + IMATIJ(ISYMI,ISYMJ) + NRHF(ISYMI)*(J - 1) + I 1512C 1513 D2IJK(NIJK) = D2IJK(NIJK) - WORK(I) 1514C 1515 150 CONTINUE 1516 140 CONTINUE 1517 130 CONTINUE 1518C 1519 CALL QEXIT('DIJK11') 1520C 1521 RETURN 1522 END 1523C /* Deck dijk13 */ 1524 SUBROUTINE DIJK13(D2IJK,ISYIJK,DENAI,ISYD1,T2TIN,ISYTIN) 1525C 1526C Written by Asger Halkier 2/2 - 1996 1527C 1528C Version: 1.0 1529C 1530C Purpose: To calculate the third contribution to D2IJK 1531C from projection against singles! 1532C 1533#include "implicit.h" 1534 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 1535 DIMENSION D2IJK(*), DENAI(*), T2TIN(*) 1536#include "ccorb.h" 1537#include "ccsdsym.h" 1538C 1539 CALL QENTER('DIJK13') 1540C 1541 IF (MULD2H(ISYTIN,ISYD1).NE.ISYIJK) 1542 & CALL QUIT('Symmetry mismatch in DIJK13') 1543C 1544C---------------------------- 1545C Calculate contribution. 1546C---------------------------- 1547C 1548 DO 100 ISYMK = 1,NSYM 1549C 1550 ISYMEI = MULD2H(ISYMK,ISYTIN) 1551 ISYMIJ = MULD2H(ISYMK,ISYIJK) 1552C 1553 DO 110 K = 1,NRHF(ISYMK) 1554C 1555 DO 120 ISYMI = 1,NSYM 1556C 1557 ISYME = MULD2H(ISYMI,ISYMEI) 1558 ISYMJ = MULD2H(ISYME,ISYD1) 1559C 1560 KOFF1 = IT2BCD(ISYMEI,ISYMK) + NT1AM(ISYMEI)*(K - 1) 1561 * + IT1AM(ISYME,ISYMI) + 1 1562 KOFF2 = IT1AM(ISYME,ISYMJ) + 1 1563 KOFF3 = IMAIJK(ISYMIJ,ISYMK) + NMATIJ(ISYMIJ)*(K - 1) 1564 * + IMATIJ(ISYMI,ISYMJ) + 1 1565C 1566 NTOTE = MAX(NVIR(ISYME),1) 1567 NTOTI = MAX(NRHF(ISYMI),1) 1568C 1569 CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),NVIR(ISYME), 1570 * -ONE,T2TIN(KOFF1),NTOTE,DENAI(KOFF2),NTOTE, 1571 * ONE,D2IJK(KOFF3),NTOTI) 1572C 1573 120 CONTINUE 1574 110 CONTINUE 1575 100 CONTINUE 1576C 1577 CALL QEXIT('DIJK13') 1578C 1579 RETURN 1580 END 1581C /* Deck daij11 */ 1582 SUBROUTINE DAIJ11(D2AIJ,ISYAIJ,DAI,ISYD1,XLAMDH,ISYMLH,IDEL,ISYMD) 1583C 1584C Written by Asger Halkier 2/2 - 1996 1585C 1586C Version: 1.0 1587C 1588C Purpose: To calculate the first contribution to D2AIJ 1589C from projection against singles! 1590C 1591#include "implicit.h" 1592 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 1593 DIMENSION D2AIJ(*), XLAMDH(*), DAI(*) 1594#include "ccorb.h" 1595#include "ccsdsym.h" 1596C 1597 CALL QENTER('DAIJ11') 1598C 1599 ISYMJ = MULD2H(ISYMD,ISYMLH) 1600 ISYMAI = ISYD1 1601 IF (MULD2H(ISYMJ,ISYMAI).NE.ISYAIJ) THEN 1602 CALL QUIT('Symmetry mismatch in DAIJ11') 1603 END IF 1604C 1605C---------------------------- 1606C Calculate contribution. 1607C---------------------------- 1608C 1609 DO 100 J = 1,NRHF(ISYMJ) 1610C 1611 NDJ = IGLMRH(ISYMD,ISYMJ) + NBAS(ISYMD)*(J - 1) + 1612 * IDEL - IBAS(ISYMD) 1613C 1614 FACT = TWO*XLAMDH(NDJ) 1615C 1616 KOFF1 = IT2BCD(ISYMAI,ISYMJ) + NT1AM(ISYMAI)*(J - 1) + 1 1617C 1618 CALL DAXPY(NT1AM(ISYMAI),FACT,DAI,1,D2AIJ(KOFF1),1) 1619C 1620 100 CONTINUE 1621C 1622 CALL QEXIT('DAIJ11') 1623C 1624 RETURN 1625 END 1626C /* Deck daij12 */ 1627 SUBROUTINE DAIJ12(D2AIJ,ISYAIJ,DAI,ISYD1,XLAMDH,ISYMLH, 1628 * WORK,LWORK,IDEL,ISYMD) 1629C 1630C Written by Asger Halkier 2/2 - 1996 1631C 1632C Version: 1.0 1633C 1634C Purpose: To calculate the second contribution to D2AIJ 1635C from projection against singles! 1636C 1637#include "implicit.h" 1638 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 1639 DIMENSION D2AIJ(*), XLAMDH(*), DAI(*), WORK(LWORK) 1640#include "priunit.h" 1641#include "ccorb.h" 1642#include "ccsdsym.h" 1643C 1644 CALL QENTER('DAIJ12') 1645C 1646 ISYMK = MULD2H(ISYMD,ISYMLH) 1647 ISYMA = MULD2H(ISYMK,ISYD1) 1648C 1649 IF (ISYMA.NE.ISYAIJ) THEN 1650 CALL QUIT('Symmetry mismatch in DAIJ12') 1651 END IF 1652C 1653 IF (LWORK .LT. NVIR(ISYMA)) THEN 1654 WRITE(LUPRI,*) 'Available:', LWORK, 'NEEDED:', NVIR(ISYMA) 1655 CALL QUIT('Insufficient space for allocation in DAIJ12') 1656 ENDIF 1657C 1658 CALL DZERO(WORK,NVIR(ISYMA)) 1659C 1660C---------------------------------------------------------- 1661C Calculate contraction of D(ai) and lamda hole matrix. 1662C---------------------------------------------------------- 1663C 1664 KOFF1 = IT1AM(ISYMA,ISYMK) + 1 1665 KOFF2 = IGLMRH(ISYMD,ISYMK) + IDEl - IBAS(ISYMD) 1666C 1667 NTOTA = MAX(NVIR(ISYMA),1) 1668C 1669 CALL DGEMV('N',NVIR(ISYMA),NRHF(ISYMK),ONE,DAI(KOFF1),NTOTA, 1670 * XLAMDH(KOFF2),NBAS(ISYMD),ZERO,WORK,1) 1671C 1672C---------------------------- 1673C Calculate contribution. 1674C---------------------------- 1675C 1676 DO 100 ISYMI = 1,NSYM 1677C 1678 ISYMJ = ISYMI 1679 ISYMAI = MULD2H(ISYMI,ISYMA) 1680C 1681 DO 110 I = 1,NRHF(ISYMI) 1682C 1683 J = I 1684 NAIJ = IT2BCD(ISYMAI,ISYMJ) + NT1AM(ISYMAI)*(J - 1) 1685 * + IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I - 1) + 1 1686C 1687 CALL DAXPY(NVIR(ISYMA),-ONE,WORK,1,D2AIJ(NAIJ),1) 1688C 1689 110 CONTINUE 1690 100 CONTINUE 1691C 1692 CALL QEXIT('DAIJ12') 1693C 1694 RETURN 1695 END 1696C /* Deck diaj11 */ 1697 SUBROUTINE DIAJ11(D2IAJ,ISYIAJ,DIA,ISYD1,XLAMDH,ISYMLH,IDEL,ISYMD) 1698C 1699C Written by Asger Halkier 4/2 - 1996 1700C 1701C Version: 1.0 1702C 1703C Purpose: To calculate the first and second contribution to D2IAJ 1704C from projection against singles! 1705C 1706#include "implicit.h" 1707 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 1708 DIMENSION D2IAJ(*), XLAMDH(*), DIA(*) 1709#include "ccorb.h" 1710#include "ccsdsym.h" 1711C 1712 CALL QENTER('DIAJ11') 1713C 1714 ISYMJ = MULD2H(ISYMD,ISYMLH) 1715 ISYMAI = ISYD1 1716 1717 IF (MULD2H(ISYMJ,ISYMAI).NE.ISYIAJ) 1718 * CALL QUIT('Symmetry mismatch in DIAJ11. (1)') 1719C 1720C-------------------------------------- 1721C Calculate the first contribution. 1722C-------------------------------------- 1723C 1724 DO 100 J = 1,NRHF(ISYMJ) 1725C 1726 NDJ = IGLMRH(ISYMD,ISYMJ) + NBAS(ISYMD)*(J - 1) + 1727 * IDEL - IBAS(ISYMD) 1728C 1729 FACT = TWO*XLAMDH(NDJ) 1730C 1731 KOFF = IT2BCD(ISYMAI,ISYMJ) + NT1AM(ISYMAI)*(J - 1) + 1 1732C 1733 CALL DAXPY(NT1AM(ISYMAI),FACT,DIA,1,D2IAJ(KOFF),1) 1734C 1735 100 CONTINUE 1736C 1737 ISYMI = MULD2H(ISYMD,ISYMLH) 1738 ISYMAJ = ISYD1 1739 1740 IF (MULD2H(ISYMI,ISYMAJ).NE.ISYIAJ) 1741 * CALL QUIT('Symmetry mismatch in DIAJ11. (2)') 1742C 1743C--------------------------------------- 1744C Calculate the second contribution. 1745C--------------------------------------- 1746C 1747 DO 110 ISYMJ = 1,NSYM 1748C 1749 ISYMA = MULD2H(ISYMJ,ISYMAJ) 1750 ISYMAI = MULD2H(ISYMA,ISYMI) 1751C 1752 DO 120 J = 1,NRHF(ISYMJ) 1753C 1754 DO 130 I = 1,NRHF(ISYMI) 1755C 1756 NDI = IGLMRH(ISYMD,ISYMI) + NBAS(ISYMD)*(I- 1) + IDEL 1757 * - IBAS(ISYMD) 1758 NAJ = IT1AM(ISYMA,ISYMJ) + NVIR(ISYMA)*(J - 1) + 1 1759 NIAJ = IT2BCD(ISYMAI,ISYMJ) + NT1AM(ISYMAI)*(J - 1) 1760 * + IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I - 1) + 1 1761C 1762 CALL DAXPY(NVIR(ISYMA),-XLAMDH(NDI),DIA(NAJ),1, 1763 * D2IAJ(NIAJ),1) 1764C 1765 130 CONTINUE 1766 120 CONTINUE 1767 110 CONTINUE 1768C 1769 CALL QEXIT('DIAJ11') 1770C 1771 RETURN 1772 END 1773C /* Deck dabi11 */ 1774 SUBROUTINE DABI11(D2ABI,ISYABI,DAI,ISYD1,T2TDEL,ISYMTI) 1775C 1776C Written by Asger Halkier 4/2 - 1996 1777C 1778C Version: 1.0 1779C 1780C Purpose: To calculate the one and only contribution to D2ABI 1781C from projection against singles! 1782C 1783#include "implicit.h" 1784 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 1785 DIMENSION D2ABI(*), DAI(*), T2TDEL(*) 1786#include "ccorb.h" 1787#include "ccsdsym.h" 1788C 1789 CALL QENTER('DABI11') 1790C 1791 IF (MULD2H(ISYMTI,ISYD1).NE.ISYABI) 1792 & CALL QUIT('Symmetry in DABI11') 1793C 1794C---------------------------- 1795C Calculate contribution. 1796C---------------------------- 1797C 1798 DO 100 ISYMI = 1,NSYM 1799C 1800 ISYMBM = MULD2H(ISYMI,ISYMTI) 1801 ISYMAB = MULD2H(ISYMI,ISYABI) 1802C 1803 DO 110 I = 1,NRHF(ISYMI) 1804C 1805 DO 120 ISYMA = 1,NSYM 1806C 1807 ISYMM = MULD2H(ISYMA,ISYD1) 1808 ISYMB = MULD2H(ISYMM,ISYMBM) 1809C 1810 KOFF1 = IT1AM(ISYMA,ISYMM) + 1 1811 KOFF2 = IT2BCD(ISYMBM,ISYMI) + NT1AM(ISYMBM)*(I - 1) 1812 * + IT1AM(ISYMB,ISYMM) + 1 1813 KOFF3 = ICKASR(ISYMAB,ISYMI) + NMATAB(ISYMAB)*(I - 1) 1814 * + IMATAB(ISYMA,ISYMB) + 1 1815C 1816 NTOTA = MAX(NVIR(ISYMA),1) 1817 NTOTB = MAX(NVIR(ISYMB),1) 1818C 1819 CALL DGEMM('N','T',NVIR(ISYMA),NVIR(ISYMB),NRHF(ISYMM), 1820 * ONE,DAI(KOFF1),NTOTA,T2TDEL(KOFF2),NTOTB,ONE, 1821 * D2ABI(KOFF3),NTOTA) 1822C 1823 120 CONTINUE 1824 110 CONTINUE 1825 100 CONTINUE 1826C 1827 CALL QEXIT('DABI11') 1828C 1829 RETURN 1830 END 1831C /* Deck dijg11 */ 1832 SUBROUTINE DIJG11(D2AOIJ,XLAMDH,Z1AO,IDEL,ISYMD,G,ISYMG) 1833C 1834C Written by Asger Halkier 4/2 - 1996 1835C 1836C Version: 1.0 1837C 1838C Purpose: To calculate the first and second contribution to D2AOIJ 1839C from projection against singles! 1840C 1841#include "implicit.h" 1842 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 1843 DIMENSION D2AOIJ(*), XLAMDH(*), Z1AO(*) 1844#include "ccorb.h" 1845#include "ccsdsym.h" 1846C 1847 CALL QENTER('DIJG11') 1848C 1849C-------------------------------------- 1850C Calculate the first contribution. 1851C-------------------------------------- 1852C 1853 IF (ISYMG .EQ. ISYMD) THEN 1854C 1855 ISYMK = ISYMG 1856C 1857 NDK = IGLMRH(ISYMD,ISYMK) + IDEL - IBAS(ISYMD) 1858 NGK = IGLMRH(ISYMG,ISYMK) + G 1859C 1860 FACT = DDOT(NRHF(ISYMK),Z1AO(NGK),NBAS(ISYMG),XLAMDH(NDK), 1861 * NBAS(ISYMD)) 1862C 1863 DO 100 ISYMI = 1,NSYM 1864C 1865 ISYMJ = ISYMI 1866C 1867 DO 110 I = 1,NRHF(ISYMI) 1868C 1869 J = I 1870 NIJ = IMATIJ(ISYMI,ISYMJ) + NRHF(ISYMI)*(J - 1) + I 1871C 1872 D2AOIJ(NIJ) = D2AOIJ(NIJ) + TWO*FACT 1873C 1874 110 CONTINUE 1875 100 CONTINUE 1876C 1877 ENDIF 1878C 1879C----------------------------------- 1880C Calculate second contribution. 1881C----------------------------------- 1882C 1883 ISYMI = ISYMD 1884 ISYMJ = ISYMG 1885C 1886 DO 120 J = 1,NRHF(ISYMJ) 1887C 1888 NGJ = IGLMRH(ISYMG,ISYMJ) + NBAS(ISYMG)*(J - 1) + G 1889C 1890 DO 130 I = 1,NRHF(ISYMI) 1891C 1892 NDI = IGLMRH(ISYMD,ISYMI) + NBAS(ISYMD)*(I - 1) + IDEL 1893 * - IBAS(ISYMD) 1894 NIJ = IMATIJ(ISYMI,ISYMJ) + NRHF(ISYMI)*(J - 1) + I 1895C 1896 D2AOIJ(NIJ) = D2AOIJ(NIJ) - XLAMDH(NDI)*Z1AO(NGJ) 1897C 1898 130 CONTINUE 1899 120 CONTINUE 1900C 1901 CALL QEXIT('DIJG11') 1902C 1903 RETURN 1904 END 1905C /* Deck diaj13 */ 1906 SUBROUTINE DIAJ13(D2IAJ,ISYIAJ,DAI,ISYD1,T2AMT,XLAMDH,ISYMLH, 1907 * WORK,LWORK,IDEL,ISYMD) 1908C 1909C Written by Asger Halkier 5/2 - 1996 1910C 1911C Version: 1.0 1912C 1913C Purpose: To calculate the third contribution to D2IAJ from 1914C projection against singles! 1915C 1916#include "implicit.h" 1917 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 1918 DIMENSION D2IAJ(*), DAI(*), T2AMT(*), XLAMDH(*), WORK(LWORK) 1919#include "priunit.h" 1920#include "ccorb.h" 1921#include "ccsdsym.h" 1922C 1923 INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 1924C 1925 CALL QENTER('DIAJ13') 1926C 1927 IF (MULD2H(ISYMLH,ISYMD).NE.MULD2H(ISYIAJ,ISYD1)) THEN 1928 CALL QUIT('Symmetry mismatch in DIAJ13.') 1929 END IF 1930C 1931 ISYMK = MULD2H(ISYMD,ISYMLH) 1932 ISYME = MULD2H(ISYMK,ISYD1) 1933C 1934C----------------------------------------- 1935C Spaghetti goto's if no contribution. 1936C----------------------------------------- 1937C 1938 IF (NRHF(ISYMK) .EQ. 0) GOTO 101 1939 IF (NVIR(ISYME) .EQ. 0) GOTO 101 1940 IF (NBAS(ISYMD) .EQ. 0) GOTO 101 1941C 1942C------------------------------- 1943C Work space allocation one. 1944C------------------------------- 1945C 1946 KVECE = 1 1947 KEND1 = KVECE + NVIR(ISYME) 1948 LWRK1 = LWORK - KEND1 1949C 1950 IF (LWRK1 .LT. 0) THEN 1951 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 1952 CALL QUIT('Insufficient work space for allocation in DIAJ13') 1953 ENDIF 1954C 1955 CALL DZERO(WORK(KVECE),NVIR(ISYME)) 1956C 1957C----------------------------------------------- 1958C Calculate contraction of zeta1 and xlamdh. 1959C----------------------------------------------- 1960C 1961 KOFF1 = IT1AM(ISYME,ISYMK) + 1 1962 KOFF2 = IGLMRH(ISYMD,ISYMK) + IDEL - IBAS(ISYMD) 1963C 1964 NTOTE = MAX(NVIR(ISYME),1) 1965C 1966 CALL DGEMV('N',NVIR(ISYME),NRHF(ISYMK),ONE,DAI(KOFF1),NTOTE, 1967 * XLAMDH(KOFF2),NBAS(ISYMD),ZERO,WORK(KVECE),1) 1968C 1969 DO 100 ISYMJ = 1,NSYM 1970C 1971 ISYMAI = MULD2H(ISYMJ,ISYIAJ) 1972 ISYMEJ = ISYMAI 1973C 1974C---------------------------------- 1975C Work space allocation two. 1976C---------------------------------- 1977C 1978 KT2SUB = KEND1 1979 KEND2 = KT2SUB + NT1AM(ISYMAI) 1980 LWRK2 = LWORK - KEND2 1981C 1982 IF (LWRK2 .LT. 0) THEN 1983 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND2 1984 CALL QUIT('Insufficient work space for '// 1985 & 'allocation in DIAJ13') 1986 ENDIF 1987C 1988 DO 110 J = 1,NRHF(ISYMJ) 1989C 1990 DO 120 E = 1,NVIR(ISYME) 1991C 1992 NEJ = IT1AM(ISYME,ISYMJ) + NVIR(ISYME)*(J - 1) + E 1993C 1994C-------------------------------------------------- 1995C Copy t1-submatrix (ai) out of t2amt. 1996C-------------------------------------------------- 1997C 1998 DO 130 NAI = 1,NT1AM(ISYMAI) 1999C 2000 NAIEJ = IT2AM(ISYMAI,ISYMEJ) + INDEX(NAI,NEJ) 2001C 2002 WORK(KT2SUB + NAI - 1) = T2AMT(NAIEJ) 2003C 2004 130 CONTINUE 2005C 2006C------------------------------------------------------ 2007C Final contraction and storage in result. 2008C------------------------------------------------------ 2009C 2010 FACT = - WORK(KVECE + E - 1) 2011C 2012 KOFF3 = IT2BCD(ISYMAI,ISYMJ) + NT1AM(ISYMAI)*(J - 1) + 1 2013C 2014 CALL DAXPY(NT1AM(ISYMAI),FACT,WORK(KT2SUB),1, 2015 * D2IAJ(KOFF3),1) 2016C 2017 120 CONTINUE 2018 110 CONTINUE 2019 100 CONTINUE 2020C 2021 101 CONTINUE 2022C 2023 CALL QEXIT('DIAJ13') 2024C 2025 RETURN 2026 END 2027C /* Deck diag11 */ 2028 SUBROUTINE DIAG11(D2AOIA,Z1AO,T2TIN,ISYMD,G,ISYMG) 2029C 2030C Written by Asger Halkier 7/2 - 1996 2031C 2032C Version: 1.0 2033C 2034C Purpose: To calculate the one and only contribution to D2AOIA 2035C from projection against singles! 2036C 2037#include "implicit.h" 2038 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 2039 DIMENSION D2AOIA(*), Z1AO(*), T2TIN(*) 2040#include "ccorb.h" 2041#include "ccsdsym.h" 2042C 2043 CALL QENTER('DIAG11') 2044C 2045 ISYAIM = ISYMD 2046 ISYMM = ISYMG 2047 ISYMAI = MULD2H(ISYMM,ISYAIM) 2048C 2049C---------------------------- 2050C Calculate contribution. 2051C---------------------------- 2052C 2053 KOFF1 = IT2BCD(ISYMAI,ISYMM) + 1 2054 KOFF2 = IGLMRH(ISYMG,ISYMM) + G 2055C 2056 NTOTAI = MAX(NT1AM(ISYMAI),1) 2057C 2058 CALL DGEMV('N',NT1AM(ISYMAI),NRHF(ISYMM),ONE,T2TIN(KOFF1), 2059 * NTOTAI,Z1AO(KOFF2),NBAS(ISYMG),ONE,D2AOIA,1) 2060C 2061 CALL QEXIT('DIAG11') 2062C 2063 RETURN 2064 END 2065C /* Deck dijk25 */ 2066 SUBROUTINE DIJK25(D2IJK,ISYIJK,XMINT,XLAMDH,ISYMLH,IDEL,ISYMD) 2067C 2068C Written by Asger Halkier 10/2 - 1996 2069C 2070C Version: 1.0 2071C 2072C Purpose: To calculate the fifth contribution to D2IJK 2073C from projection against doubles! 2074C 2075#include "implicit.h" 2076 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 2077 DIMENSION D2IJK(*), XMINT(*), XLAMDH(*) 2078#include "ccorb.h" 2079#include "ccsdsym.h" 2080C 2081 CALL QENTER('DIJK25') 2082C 2083 ISYML = MULD2H(ISYMD,ISYMLH) 2084 IF (ISYIJK.NE.ISYML) CALL QUIT('Symmetry mismatch in DIJK25') 2085C 2086C-------------------------------- 2087C Calculate the contribution. 2088C-------------------------------- 2089C 2090 KOFF1 = I3ORHF(ISYIJK,ISYML) + 1 2091 KOFF2 = IGLMRH(ISYMD,ISYML) + IDEL - IBAS(ISYMD) 2092C 2093 NTOIJK = MAX(NMAIJK(ISYIJK),1) 2094C 2095 CALL DGEMV('N',NMAIJK(ISYIJK),NRHF(ISYML),ONE,XMINT(KOFF1), 2096 * NTOIJK,XLAMDH(KOFF2),NBAS(ISYMD),ONE,D2IJK,1) 2097C 2098 CALL QEXIT('DIJK25') 2099C 2100 RETURN 2101 END 2102C /* Deck dijk24 */ 2103 SUBROUTINE DIJK24(D2IJK,ISYIJK,XMAT,XLAMDH,ISYMLH,IDEL,ISYMD) 2104C 2105C Written by Asger Halkier 10/2 - 1996 2106C 2107C Version: 1.0 2108C 2109C Purpose: To calculate the fourth contribution to D2IJK 2110C from projection against doubles! 2111C 2112#include "implicit.h" 2113 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 2114 DIMENSION D2IJK(*), XMAT(*), XLAMDH(*) 2115#include "ccorb.h" 2116#include "ccsdsym.h" 2117C 2118 CALL QENTER('DIJK24') 2119C 2120 ISYMI = MULD2H(ISYMD,ISYMLH) 2121 ISYMKJ = 1 2122C 2123 IF (ISYMI.NE.ISYIJK) CALL QUIT('Symmetry mismatch in DIJK24') 2124C 2125C-------------------------------- 2126C Calculate the contribution. 2127C-------------------------------- 2128C 2129 DO 100 ISYMK = 1,NSYM 2130C 2131 ISYMJ = MULD2H(ISYMK,ISYMKJ) 2132 ISYMIJ = MULD2H(ISYMI,ISYMJ) 2133C 2134 DO 110 K = 1,NRHF(ISYMK) 2135C 2136 DO 120 J = 1,NRHF(ISYMJ) 2137C 2138 NKJ = IMATIJ(ISYMK,ISYMJ) + NRHF(ISYMK)*(J - 1) + K 2139 NDEI = IGLMRH(ISYMD,ISYMI) + IDEL - IBAS(ISYMD) 2140 NIJK = IMAIJK(ISYMIJ,ISYMK) + NMATIJ(ISYMIJ)*(K - 1) 2141 * + IMATIJ(ISYMI,ISYMJ) + NRHF(ISYMI)*(J - 1) + 1 2142C 2143 CALL DAXPY(NRHF(ISYMI),XMAT(NKJ),XLAMDH(NDEI), 2144 * NBAS(ISYMD),D2IJK(NIJK),1) 2145C 2146 120 CONTINUE 2147 110 CONTINUE 2148 100 CONTINUE 2149C 2150 CALL QEXIT('DIJK24') 2151C 2152 RETURN 2153 END 2154C /* Deck dijk23 */ 2155 SUBROUTINE DIJK23(D2IJK,ISYIJK,XMAT,XLAMDH,ISYMLH,IDEL,ISYMD) 2156C 2157C Written by Asger Halkier 10/2 - 1996 2158C 2159C Version: 1.0 2160C 2161C Purpose: To calculate the third contribution to D2IJK 2162C from projection against doubles! 2163C 2164#include "implicit.h" 2165 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 2166 DIMENSION D2IJK(*), XMAT(*), XLAMDH(*) 2167#include "ccorb.h" 2168#include "ccsdsym.h" 2169C 2170 CALL QENTER('DIJK23') 2171C 2172 ISYMK = MULD2H(ISYMD,ISYMLH) 2173 ISYMIJ = 1 2174C 2175 IF (ISYMK.NE.ISYIJK) CALL QUIT('Symmetry mismatch in DIJK23') 2176C 2177C-------------------------------- 2178C Calculate the contribution. 2179C-------------------------------- 2180C 2181 DO 100 K = 1,NRHF(ISYMK) 2182C 2183 NDEK = IGLMRH(ISYMD,ISYMK) + NBAS(ISYMD)*(K - 1) 2184 * + IDEL - IBAS(ISYMD) 2185 NIJK = IMAIJK(ISYMIJ,ISYMK) + NMATIJ(ISYMIJ)*(K - 1) + 1 2186C 2187 FACT = -TWO*XLAMDH(NDEK) 2188C 2189 CALL DAXPY(NMATIJ(ISYMIJ),FACT,XMAT,1,D2IJK(NIJK),1) 2190C 2191 100 CONTINUE 2192C 2193 CALL QEXIT('DIJK23') 2194C 2195 RETURN 2196 END 2197C /* Deck dijk21 */ 2198 SUBROUTINE DIJK21(D2IJK,ISYIJK,XMAT,XLAMDH,ISYMLH, 2199 * WORK,LWORK,IDEL,ISYMD) 2200C 2201C Written by Asger Halkier 10/2 - 1996 2202C 2203C Version: 1.0 2204C 2205C Purpose: To calculate the first and second contribution to D2IJK 2206C from projection against doubles! 2207C 2208#include "implicit.h" 2209 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 2210 DIMENSION D2IJK(*), XMAT(*), XLAMDH(*), WORK(LWORK) 2211#include "priunit.h" 2212#include "ccorb.h" 2213#include "ccsdsym.h" 2214C 2215 CALL QENTER('DIJK21') 2216C 2217 ISYML = MULD2H(ISYMD,ISYMLH) 2218 ISYMX1 = ISYML 2219C 2220 IF (ISYML.NE.ISYIJK) CALL QUIT('Symmetry mismatch in DIJK21') 2221C 2222C----------------------------------------------------------------- 2223C Calculate contraction of X-intermediate & lamda hole matrix. 2224C----------------------------------------------------------------- 2225C 2226 IF (LWORK .LT. NRHF(ISYMX1)) THEN 2227 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', NRHF(ISYMX1) 2228 CALL QUIT('Insufficient work space in DIJK21') 2229 ENDIF 2230C 2231 CALL DZERO(WORK,NRHF(ISYMX1)) 2232C 2233 KOFF1 = IMATIJ(ISYMX1,ISYML) + 1 2234 KOFF2 = IGLMRH(ISYMD,ISYML) + IDEL - IBAS(ISYMD) 2235C 2236 NTOTI = MAX(NRHF(ISYMX1),1) 2237C 2238 CALL DGEMV('N',NRHF(ISYMX1),NRHF(ISYML),ONE,XMAT(KOFF1),NTOTI, 2239 * XLAMDH(KOFF2),NBAS(ISYMD),ZERO,WORK,1) 2240C 2241C---------------------------------- 2242C Calculate first contribution. 2243C---------------------------------- 2244C 2245 ISYMI = ISYMX1 2246C 2247 DO 100 ISYMK = 1,NSYM 2248C 2249 ISYMJ = ISYMK 2250 ISYMIJ = MULD2H(ISYMJ,ISYMI) 2251C 2252 DO 110 K = 1,NRHF(ISYMK) 2253C 2254 J = K 2255 NIJK = IMAIJK(ISYMIJ,ISYMK) + NMATIJ(ISYMIJ)*(K - 1) 2256 * + IMATIJ(ISYMI,ISYMJ) + NRHF(ISYMI)*(J - 1) + 1 2257C 2258 CALL DAXPY(NRHF(ISYMI),ONE,WORK,1,D2IJK(NIJK),1) 2259C 2260 110 CONTINUE 2261 100 CONTINUE 2262C 2263C----------------------------------- 2264C Calculate second contribution. 2265C----------------------------------- 2266C 2267 ISYMK = ISYMX1 2268 ISYMIJ = 1 2269C 2270 DO 120 K = 1,NRHF(ISYMK) 2271C 2272 DO 130 ISYMJ = 1,NSYM 2273C 2274 ISYMI = ISYMJ 2275C 2276 DO 140 J = 1,NRHF(ISYMJ) 2277C 2278 I = J 2279 NIJK = IMAIJK(ISYMIJ,ISYMK) + NMATIJ(ISYMIJ)*(K - 1) 2280 * + IMATIJ(ISYMI,ISYMJ) + NRHF(ISYMI)*(J - 1) + I 2281C 2282 D2IJK(NIJK) = D2IJK(NIJK) - TWO*WORK(K) 2283C 2284 140 CONTINUE 2285 130 CONTINUE 2286 120 CONTINUE 2287C 2288 CALL QEXIT('DIJK21') 2289C 2290 RETURN 2291 END 2292C /* Deck daij22 */ 2293 SUBROUTINE DAIJ22(D2AIJ,ISYAIJ,DAB,ISYD1,XLAMDH,ISYMLH, 2294 * WORK,LWORK,IDEL,ISYMD) 2295C 2296C Written by Asger Halkier 10/2 - 1996 2297C 2298C Version: 1.0 2299C 2300C Purpose: To calculate the second contribution to D2AIJ 2301C from projection against doubles! 2302C 2303#include "implicit.h" 2304 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 2305 DIMENSION D2AIJ(*), DAB(*), XLAMDH(*), WORK(LWORK) 2306#include "priunit.h" 2307#include "ccorb.h" 2308#include "ccsdsym.h" 2309C 2310 CALL QENTER('DAIJ22') 2311C 2312 ISYMB = MULD2H(ISYMD,ISYMLH) 2313 ISYMA = MULD2H(ISYMB,ISYD1) 2314 ISYMIJ = 1 2315C 2316 IF (ISYMA.NE.ISYAIJ) CALL QUIT('Symmetry mismatch in DAIJ22') 2317C 2318C----------------------------------------------------------------- 2319C Calculate contraction of Y-intermediate & lamda hole matrix. 2320C----------------------------------------------------------------- 2321C 2322 IF (LWORK .LT. NVIR(ISYMA)) THEN 2323 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', NVIR(ISYMA) 2324 CALL QUIT('Insufficient work space in DAIJ22') 2325 ENDIF 2326C 2327 CALL DZERO(WORK,NVIR(ISYMA)) 2328C 2329 KOFF1 = IMATAB(ISYMA,ISYMB) + 1 2330 KOFF2 = IGLMVI(ISYMD,ISYMB) + IDEL - IBAS(ISYMD) 2331C 2332 NTOTA = MAX(NVIR(ISYMA),1) 2333C 2334 CALL DGEMV('N',NVIR(ISYMA),NVIR(ISYMB),ONE,DAB(KOFF1),NTOTA, 2335 * XLAMDH(KOFF2),NBAS(ISYMD),ZERO,WORK,1) 2336C 2337C------------------------------------- 2338C Calculate contribution to D2AIJ. 2339C------------------------------------- 2340C 2341 DO 100 ISYMJ = 1,NSYM 2342C 2343 ISYMI = ISYMJ 2344 ISYMAI = MULD2H(ISYMI,ISYMA) 2345C 2346 DO 110 J = 1,NRHF(ISYMJ) 2347C 2348 I = J 2349 NAIJ = IT2BCD(ISYMAI,ISYMJ) + NT1AM(ISYMAI)*(J - 1) 2350 * + IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I - 1) + 1 2351C 2352 CALL DAXPY(NVIR(ISYMA),-ONE,WORK,1,D2AIJ(NAIJ),1) 2353C 2354 110 CONTINUE 2355 100 CONTINUE 2356C 2357 CALL QEXIT('DAIJ22') 2358C 2359 RETURN 2360 END 2361C /* Deck dabi21 */ 2362 SUBROUTINE DABI21(D2ABI,ISYABI,DAB,ISYD1,XLAMDH,ISYMLH,IDEL,ISYMD) 2363C 2364C Written by Asger Halkier 11/2 - 1996 2365C 2366C Version: 1.0 2367C 2368C Purpose: To calculate the first contribution to D2ABI 2369C from projection against doubles! 2370C 2371#include "implicit.h" 2372 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 2373 DIMENSION D2ABI(*), DAB(*), XLAMDH(*) 2374#include "ccorb.h" 2375#include "ccsdsym.h" 2376C 2377 CALL QENTER('DABI21') 2378C 2379 ISYMI = MULD2H(ISYMD,ISYMLH) 2380 ISYMAB = ISYD1 2381C 2382 IF (MULD2H(ISYMI,ISYMAB).NE.ISYABI) 2383 * CALL QUIT('Symmetry mismatch in DABI21') 2384C 2385C-------------------------------- 2386C Calculate the contribution. 2387C-------------------------------- 2388C 2389 DO 100 I = 1,NRHF(ISYMI) 2390C 2391 NDEI = IGLMRH(ISYMD,ISYMI) + NBAS(ISYMD)*(I - 1) 2392 * + IDEL - IBAS(ISYMD) 2393 NABI = ICKASR(ISYMAB,ISYMI) + NMATAB(ISYMAB)*(I - 1) + 1 2394C 2395 FACT = TWO*XLAMDH(NDEI) 2396C 2397 CALL DAXPY(NMATAB(ISYMAB),FACT,DAB,1,D2ABI(NABI),1) 2398C 2399 100 CONTINUE 2400C 2401 CALL QEXIT('DABI21') 2402C 2403 RETURN 2404 END 2405C /* Deck diaj21 */ 2406 SUBROUTINE DIAJ21(D2IAJ,ISYAIJ,YMAT,T2TINT,ISYMTI) 2407C 2408C Written by Asger Halkier 11/2 - 1996 2409C 2410C Version: 1.0 2411C 2412C Purpose: To calculate the first contribution to D2IAJ 2413C from projection against doubles! 2414C 2415#include "implicit.h" 2416 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 2417 DIMENSION D2IAJ(*), YMAT(*), T2TINT(*) 2418#include "ccorb.h" 2419#include "ccsdsym.h" 2420C 2421 CALL QENTER('DIAJ21') 2422C 2423 ISYMAE = 1 2424C 2425 IF (ISYMTI.NE.ISYAIJ) CALL QUIT('Symmetry mismatch in DIAJ21.') 2426C 2427C-------------------------------- 2428C Calculate the contribution. 2429C-------------------------------- 2430C 2431 DO 100 ISYMJ = 1,NSYM 2432C 2433 ISYMEI = MULD2H(ISYMJ,ISYMTI) 2434 ISYMAI = MULD2H(ISYMJ,ISYAIJ) 2435C 2436 DO 110 J = 1,NRHF(ISYMJ) 2437C 2438 DO 120 ISYMA = 1,NSYM 2439C 2440 ISYME = MULD2H(ISYMA,ISYMAE) 2441 ISYMI = MULD2H(ISYMA,ISYMAI) 2442C 2443 KOFF1 = IMATAB(ISYMA,ISYME) + 1 2444 KOFF2 = IT2BCD(ISYMEI,ISYMJ) + NT1AM(ISYMEI)*(J - 1) 2445 * + IT1AM(ISYME,ISYMI) + 1 2446 KOFF3 = IT2BCD(ISYMAI,ISYMJ) + NT1AM(ISYMAI)*(J - 1) 2447 * + IT1AM(ISYMA,ISYMI) + 1 2448C 2449 NTOTA = MAX(NVIR(ISYMA),1) 2450 NTOTE = MAX(NVIR(ISYME),1) 2451C 2452 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI),NVIR(ISYME), 2453 * -ONE,YMAT(KOFF1),NTOTA,T2TINT(KOFF2),NTOTE, 2454 * ONE,D2IAJ(KOFF3),NTOTA) 2455C 2456 120 CONTINUE 2457 110 CONTINUE 2458 100 CONTINUE 2459C 2460 CALL QEXIT('DIAJ21') 2461C 2462 RETURN 2463 END 2464C /* Deck diaj22 */ 2465 SUBROUTINE DIAJ22(D2IAJ,ISYAIJ,XMAT,T2TINT,ISYMTI) 2466C 2467C Written by Asger Halkier 11/2 - 1996 2468C 2469C Version: 1.0 2470C 2471C Purpose: To calculate the second contribution to D2IAJ 2472C from projection against doubles! 2473C 2474#include "implicit.h" 2475 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 2476 DIMENSION D2IAJ(*), XMAT(*), T2TINT(*) 2477#include "ccorb.h" 2478#include "ccsdsym.h" 2479C 2480 CALL QENTER('DIAJ22') 2481C 2482 IF (ISYAIJ.NE.ISYMTI) CALL QUIT('Symmetry mismatch in DIAJ22.') 2483C 2484C-------------------------------- 2485C Calculate the contribution. 2486C-------------------------------- 2487C 2488 DO 100 ISYMJ = 1,NSYM 2489C 2490 ISYMAM = MULD2H(ISYMJ,ISYMTI) 2491 ISYMAI = MULD2H(ISYMJ,ISYAIJ) 2492C 2493 DO 110 J = 1,NRHF(ISYMJ) 2494C 2495 DO 120 ISYMA = 1,NSYM 2496C 2497 ISYMM = MULD2H(ISYMA,ISYMAM) 2498 ISYMI = MULD2H(ISYMA,ISYMAI) 2499C 2500 KOFF1 = IT2BCD(ISYMAM,ISYMJ) + NT1AM(ISYMAM)*(J - 1) 2501 * + IT1AM(ISYMA,ISYMM) + 1 2502 KOFF2 = IMATIJ(ISYMI,ISYMM) + 1 2503 KOFF3 = IT2BCD(ISYMAI,ISYMJ) + NT1AM(ISYMAI)*(J - 1) 2504 * + IT1AM(ISYMA,ISYMI) + 1 2505C 2506 NTOTA = MAX(NVIR(ISYMA),1) 2507 NTOTI = MAX(NRHF(ISYMI),1) 2508C 2509 CALL DGEMM('N','T',NVIR(ISYMA),NRHF(ISYMI),NRHF(ISYMM), 2510 * -ONE,T2TINT(KOFF1),NTOTA,XMAT(KOFF2),NTOTI, 2511 * ONE,D2IAJ(KOFF3),NTOTA) 2512C 2513 120 CONTINUE 2514 110 CONTINUE 2515 100 CONTINUE 2516C 2517 CALL QEXIT('DIAJ22') 2518C 2519 RETURN 2520 END 2521C /* Deck diaj23 */ 2522 SUBROUTINE DIAJ23(D2IAJ,ISYAIJ,XMAT,T2TINT,ISYMTI) 2523C 2524C Written by Asger Halkier 11/2 - 1996 2525C 2526C Version: 1.0 2527C 2528C Purpose: To calculate the third contribution to D2IAJ 2529C from projection against doubles! 2530C 2531#include "implicit.h" 2532 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 2533 DIMENSION D2IAJ(*), XMAT(*), T2TINT(*) 2534#include "ccorb.h" 2535#include "ccsdsym.h" 2536C 2537 CALL QENTER('DIAJ23') 2538C 2539 ISYMJM = 1 2540C 2541 IF (ISYMTI.NE.ISYAIJ) CALL QUIT('Symmetry mismatch in DIAJ23') 2542C 2543C-------------------------------- 2544C Calculate the contribution. 2545C-------------------------------- 2546C 2547 DO 100 ISYMAI = 1,NSYM 2548C 2549 ISYMM = MULD2H(ISYMAI,ISYMTI) 2550 ISYMJ = MULD2H(ISYMAI,ISYAIJ) 2551C 2552 KOFF1 = IT2BCD(ISYMAI,ISYMM) + 1 2553 KOFF2 = IMATIJ(ISYMJ,ISYMM) + 1 2554 KOFF3 = IT2BCD(ISYMAI,ISYMJ) + 1 2555C 2556 NTOAI = MAX(NT1AM(ISYMAI),1) 2557 NTOTJ = MAX(NRHF(ISYMJ),1) 2558C 2559 CALL DGEMM('N','T',NT1AM(ISYMAI),NRHF(ISYMJ),NRHF(ISYMM),-ONE, 2560 * T2TINT(KOFF1),NTOAI,XMAT(KOFF2),NTOTJ,ONE, 2561 * D2IAJ(KOFF3),NTOAI) 2562C 2563 100 CONTINUE 2564C 2565 CALL QEXIT('DIAJ23') 2566C 2567 RETURN 2568 END 2569C /* Deck diaj24 */ 2570 SUBROUTINE DIAJ24(D2IAJ,ISYIAJ,T2AMT,DAB,ISYD1,XLAMDH,ISYMLH, 2571 * WORK,LWORK,IDEL,ISYMD) 2572C 2573C Written by Asger Halkier 13/2 - 1996 2574C 2575C Version: 1.0 2576C 2577C Purpose: To calculate the fourth contribution to D2IAJ 2578C from projection against doubles! 2579C 2580#include "implicit.h" 2581 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 2582 DIMENSION D2IAJ(*), T2AMT(*), DAB(*), XLAMDH(*), WORK(LWORK) 2583#include "priunit.h" 2584#include "ccorb.h" 2585#include "ccsdsym.h" 2586C 2587 INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 2588C 2589 CALL QENTER('DIAJ24') 2590C 2591 ISYMB = MULD2H(ISYMD,ISYMLH) 2592 ISYMC = MULD2H(ISYMB,ISYD1) 2593C 2594 IF (MULD2H(ISYMLH,ISYMD).NE.MULD2H(ISYD1,ISYIAJ)) THEN 2595 CALL QUIT('Symmetry mismatch in DIAJ24.') 2596 END IF 2597C 2598 IF (NVIR(ISYMB) .EQ. 0) THEN 2599 CALL QEXIT('DIAJ24') 2600 RETURN 2601 END IF 2602C 2603C------------------------------- 2604C Work space allocation one. 2605C------------------------------- 2606C 2607 KVECC = 1 2608 KEND1 = KVECC + NVIR(ISYMC) 2609 LWRK1 = LWORK - KEND1 2610C 2611 IF (LWRK1 .LT. 0) THEN 2612 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 2613 CALL QUIT('Insufficient work space for allocation in DIAJ24') 2614 ENDIF 2615C 2616 CALL DZERO(WORK(KVECC),NVIR(ISYMC)) 2617C 2618C------------------------------------------------------------------- 2619C Calculate contraction of transposed Y-matrix (DAB) and XLAMDH. 2620C------------------------------------------------------------------- 2621C 2622 KOFF1 = IMATAB(ISYMC,ISYMB) + 1 2623 KOFF2 = IGLMVI(ISYMD,ISYMB) + IDEL - IBAS(ISYMD) 2624C 2625 NTOTC = MAX(NVIR(ISYMC),1) 2626C 2627 CALL DGEMV('N',NVIR(ISYMC),NVIR(ISYMB),ONE,DAB(KOFF1),NTOTC, 2628 * XLAMDH(KOFF2),NBAS(ISYMD),ZERO,WORK(KVECC),1) 2629C 2630C-------------------------------- 2631C Calculate the contribution. 2632C-------------------------------- 2633C 2634 DO 100 ISYMJ = 1,NSYM 2635C 2636 ISYMAI = MULD2H(ISYMJ,ISYIAJ) 2637 ISYMCJ = ISYMAI 2638C 2639 DO 110 J = 1,NRHF(ISYMJ) 2640C 2641 DO 120 C = 1,NVIR(ISYMC) 2642C 2643 NCJ = IT1AM(ISYMC,ISYMJ) + NVIR(ISYMC)*(J - 1) + C 2644C 2645 DO 130 NAI = 1,NT1AM(ISYMAI) 2646C 2647 NAICJ = IT2AM(ISYMAI,ISYMCJ) + INDEX(NAI,NCJ) 2648 NAIJ = IT2BCD(ISYMAI,ISYMJ) + NT1AM(ISYMAI)*(J - 1) 2649 * + NAI 2650C 2651 D2IAJ(NAIJ) = D2IAJ(NAIJ) 2652 * - T2AMT(NAICJ)*WORK(KVECC + C - 1) 2653C 2654 130 CONTINUE 2655 120 CONTINUE 2656 110 CONTINUE 2657 100 CONTINUE 2658C 2659 CALL QEXIT('DIAJ24') 2660C 2661 RETURN 2662 END 2663C /* Deck ccd1reor */ 2664 SUBROUTINE CCD1REOR(DFOCK,DAI,DAB,DIJ,DIA) 2665C 2666C Written by Asger Halkier 28/2 - 1996 2667C 2668C Version: 1.0 2669C 2670C Purpose: To reorder the four separate blocks DAI through DIA 2671C of the CC one electron density to one matrix stored 2672C as a directly diagonalisable Fock-matrix (DFOCK)! 2673C 2674#include "implicit.h" 2675 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 2676 DIMENSION DFOCK(*), DAI(*), DAB(*), DIJ(*), DIA(*) 2677#include "ccorb.h" 2678#include "ccsdsym.h" 2679C 2680 CALL QENTER('CCD1REOR') 2681C 2682 DO 100 ISYM = 1,NSYM 2683C 2684C--------------------------------------- 2685C Do the occupied occupied block. 2686C--------------------------------------- 2687C 2688 DO 110 J = 1,NRHF(ISYM) 2689C 2690 NIJOR = IMATIJ(ISYM,ISYM) + NRHF(ISYM)*(J - 1) + 1 2691 NIJNE = IFCKDO(ISYM) + NORB(ISYM)*(J - 1) + 1 2692C 2693 CALL DCOPY(NRHF(ISYM),DIJ(NIJOR),1,DFOCK(NIJNE),1) 2694C 2695 110 CONTINUE 2696C 2697C-------------------------------------- 2698C Do the virtual occupied block. 2699C-------------------------------------- 2700C 2701 DO 120 I = 1,NRHF(ISYM) 2702C 2703 NAIOR = IT1AM(ISYM,ISYM) + NVIR(ISYM)*(I - 1) + 1 2704 NAINE = IFCKDO(ISYM) + NORB(ISYM)*(I - 1) 2705 * + NRHF(ISYM) + 1 2706C 2707 CALL DCOPY(NVIR(ISYM),DAI(NAIOR),1,DFOCK(NAINE),1) 2708C 2709 120 CONTINUE 2710C 2711C-------------------------------------- 2712C Do the occupied virtual block. 2713C-------------------------------------- 2714C 2715 DO 130 A = 1,NVIR(ISYM) 2716C 2717 NIAOR = IT1AM(ISYM,ISYM) + A 2718 NIANE = IFCKDV(ISYM) + NORB(ISYM)*(A - 1) + 1 2719C 2720 CALL DCOPY(NRHF(ISYM),DIA(NIAOR),NVIR(ISYM), 2721 * DFOCK(NIANE),1) 2722C 2723 130 CONTINUE 2724C 2725C------------------------------------- 2726C Do the virtual virtual block. 2727C------------------------------------- 2728C 2729 DO 140 B = 1,NVIR(ISYM) 2730C 2731 NABOR = IMATAB(ISYM,ISYM) + NVIR(ISYM)*(B - 1) + 1 2732 NABNE = IFCKDV(ISYM) + NORB(ISYM)*(B - 1) 2733 * + NRHF(ISYM) + 1 2734C 2735 CALL DCOPY(NVIR(ISYM),DAB(NABOR),1,DFOCK(NABNE),1) 2736C 2737 140 CONTINUE 2738 100 CONTINUE 2739C 2740 CALL QEXIT('CCD1REOR') 2741C 2742 RETURN 2743 END 2744C /* Deck ccd1t1co */ 2745 SUBROUTINE CCD1T1CO(DAI,DAB,DIJ,DIA,T1AM,YMAT,WORK,LWORK) 2746C 2747C Written by Asger Halkier 29/2 - 1996 2748C 2749C Version: 1.0 2750C 2751C Purpose: To add contributions from the t1-amplitudes to 2752C the CC one electron density in MO-basis. 2753C -- 2754C 2755#include "implicit.h" 2756 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 2757 DIMENSION DAI(*), DAB(*), DIJ(*), DIA(*) 2758 DIMENSION T1AM(*), YMAT(*), WORK(LWORK) 2759#include "priunit.h" 2760#include "ccorb.h" 2761#include "ccsdsym.h" 2762C 2763 CALL QENTER('CCD1T1CO') 2764C 2765C---------------------------------------------- 2766C Add contributions to the (occ,vir) block. 2767C---------------------------------------------- 2768C 2769 DO 100 ISYMA = 1,NSYM 2770C 2771 ISYMI = ISYMA 2772 ISYMK = ISYMA 2773 ISYMC = ISYMA 2774C 2775 KOFF1 = IT1AM(ISYMA,ISYMK) + 1 2776 KOFF2 = IMATIJ(ISYMI,ISYMK) + 1 2777 KOFF3 = IT1AM(ISYMA,ISYMI) + 1 2778C 2779 NTOTA = MAX(NVIR(ISYMA),1) 2780 NTOTI = MAX(NRHF(ISYMI),1) 2781C 2782 CALL DGEMM('N','T',NVIR(ISYMA),NRHF(ISYMI),NRHF(ISYMK), 2783 * ONE,T1AM(KOFF1),NTOTA,DIJ(KOFF2),NTOTI,ONE, 2784 * DIA(KOFF3),NTOTA) 2785C 2786 KOFF4 = IMATAB(ISYMA,ISYMC) + 1 2787 KOFF5 = IT1AM(ISYMC,ISYMI) + 1 2788 KOFF6 = IT1AM(ISYMA,ISYMI) + 1 2789C 2790 NTOTA = MAX(NVIR(ISYMA),1) 2791 NTOTC = MAX(NVIR(ISYMC),1) 2792C 2793 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI),NVIR(ISYMC), 2794 * -ONE,YMAT(KOFF4),NTOTA,T1AM(KOFF5),NTOTC,ONE, 2795 * DIA(KOFF6),NTOTA) 2796C 2797C------------------------------ 2798C Work space allocation. 2799C------------------------------ 2800C 2801 KSCR = 1 2802 KEND1 = KSCR + NRHF(ISYMK)*NRHF(ISYMI) 2803 LWRK1 = LWORK - KEND1 2804C 2805 IF (LWRK1 .LT. 0) THEN 2806 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 2807 CALL QUIT('Insufficient work space for '// 2808 & 'allocation in CCD1T1CO') 2809 ENDIF 2810C 2811 CALL DZERO(WORK(KSCR),NRHF(ISYMK)*NRHF(ISYMI)) 2812C 2813 KOFF7 = IT1AM(ISYMC,ISYMK) + 1 2814 KOFF8 = IT1AM(ISYMC,ISYMI) + 1 2815C 2816 NTOTC = MAX(NVIR(ISYMC),1) 2817 NTOTK = MAX(NRHF(ISYMK),1) 2818C 2819 CALL DGEMM('T','N',NRHF(ISYMK),NRHF(ISYMI),NVIR(ISYMC), 2820 * ONE,DAI(KOFF7),NTOTC,T1AM(KOFF8),NTOTC,ZERO, 2821 * WORK(KSCR),NTOTK) 2822C 2823 KOFF9 = IT1AM(ISYMA,ISYMK) + 1 2824 KOFF10 = IT1AM(ISYMA,ISYMI) + 1 2825C 2826 NTOTA = MAX(NVIR(ISYMA),1) 2827 NTOTK = MAX(NRHF(ISYMK),1) 2828C 2829 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI),NRHF(ISYMK), 2830 * -ONE,T1AM(KOFF9),NTOTA,WORK(KSCR),NTOTK,ONE, 2831 * DIA(KOFF10),NTOTA) 2832C 2833 100 CONTINUE 2834C 2835C---------------------------------------------- 2836C Add contributions to the (vir,vir) block. 2837C---------------------------------------------- 2838C 2839 DO 110 ISYMA = 1,NSYM 2840C 2841 ISYMB = ISYMA 2842 ISYMK = ISYMA 2843C 2844 KOFF11 = IT1AM(ISYMA,ISYMK) + 1 2845 KOFF12 = IT1AM(ISYMB,ISYMK) + 1 2846 KOFF13 = IMATAB(ISYMA,ISYMB) + 1 2847C 2848 NTOTA = MAX(NVIR(ISYMA),1) 2849 NTOTB = MAX(NVIR(ISYMB),1) 2850C 2851 CALL DGEMM('N','T',NVIR(ISYMA),NVIR(ISYMB),NRHF(ISYMK), 2852 * ONE,DAI(KOFF11),NTOTA,T1AM(KOFF12),NTOTB,ONE, 2853 * DAB(KOFF13),NTOTA) 2854C 2855 110 CONTINUE 2856C 2857C---------------------------------------------- 2858C Add contributions to the (occ,occ) block. 2859C---------------------------------------------- 2860C 2861 DO 120 ISYMI = 1,NSYM 2862C 2863 ISYMJ = ISYMI 2864 ISYMC = ISYMI 2865C 2866 KOFF14 = IT1AM(ISYMC,ISYMI) + 1 2867 KOFF15 = IT1AM(ISYMC,ISYMJ) + 1 2868 KOFF16 = IMATIJ(ISYMI,ISYMJ) + 1 2869C 2870 NTOTC = MAX(NVIR(ISYMC),1) 2871 NTOTI = MAX(NRHF(ISYMI),1) 2872C 2873 CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),NVIR(ISYMC), 2874 * -ONE,T1AM(KOFF14),NTOTC,DAI(KOFF15),NTOTC,ONE, 2875 * DIJ(KOFF16),NTOTI) 2876C 2877 120 CONTINUE 2878C 2879 CALL QEXIT('CCD1T1CO') 2880C 2881 RETURN 2882 END 2883C /* Deck sortash */ 2884 SUBROUTINE SORTASH(XARR,YARR,N) 2885C 2886C Written by Asger Halkier 4/2 - 1996 2887C 2888C Version: 1.0 2889C 2890C Purpose: To reorder the array xarr after size, and 2891C do the same reordering of the concomitant 2892C array yarr! 2893C (This routine is based on Bubble-sort from NUM.REC.) 2894C 2895#include "implicit.h" 2896 DIMENSION XARR(*), YARR(*) 2897#include "ccorb.h" 2898#include "ccsdsym.h" 2899C 2900 CALL QENTER('SORTASH') 2901C 2902C---------------------- 2903C Do the resorting. 2904C---------------------- 2905C 2906 DO 100 J = 2,N 2907C 2908 XA = XARR(J) 2909 XB = YARR(J) 2910C 2911 DO 110 I = J-1,1,-1 2912C 2913 IF (XARR(I) .LE. XA) GOTO 120 2914C 2915 XARR(I + 1) = XARR(I) 2916 YARR(I + 1) = YARR(I) 2917C 2918 110 CONTINUE 2919C 2920 I=0 2921C 2922 120 XARR(I + 1) = XA 2923 YARR(I + 1) = XB 2924C 2925 100 CONTINUE 2926C 2927 CALL QEXIT('SORTASH') 2928C 2929 RETURN 2930 END 2931C /* Deck ccnaocan */ 2932 SUBROUTINE CCNAOCAN(XARR,YARR) 2933C 2934C Written by Asger Halkier 4/2 - 1996 2935C 2936C Version: 1.0 2937C 2938C Purpose: To analyze and check the natural occupation 2939C numbers in xarr and yarr (real and imag)! 2940C 2941#include "implicit.h" 2942 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 2943 DIMENSION XARR(*), YARR(*) 2944#include "priunit.h" 2945#include "ccorb.h" 2946#include "ccsdsym.h" 2947#include "ccsdinp.h" 2948C 2949 CALL QENTER('CCNAOCAN') 2950C 2951 SUMRHF = ZERO 2952 SUMALL = ZERO 2953 SUMIMA = ZERO 2954C 2955 DO 100 I = 1,NORBT 2956C 2957 SUMALL = SUMALL + XARR(I) 2958 SUMIMA = SUMIMA + YARR(I) 2959C 2960 100 CONTINUE 2961C 2962 CHECK = 2*NRHFT - SUMALL 2963 CHEPR = NRHFT*2.0D0 2964 THR = 1.0D-06 2965C 2966 IF (CHECK .LT. THR) THEN 2967C 2968 WRITE(LUPRI,*) ' ' 2969 WRITE(LUPRI,*) ' ' 2970 WRITE(LUPRI,111) 'Total Sum of natural occupation numbers:', 2971 & CHEPR 2972C 2973 ELSE 2974C 2975 WRITE(LUPRI,*) ' ' 2976 WRITE(LUPRI,*) ' ' 2977 WRITE(LUPRI,111) 'Total Sum of natural occupation numbers:', 2978 & SUMALL 2979C 2980 ENDIF 2981C 2982 IF (IPRINT .GT. 20) THEN 2983C 2984 WRITE(LUPRI,*) ' ' 2985 WRITE(LUPRI,111) 'Total Sum of imaginary parts of nat.occ:', 2986 & SUMIMA 2987C 2988 ENDIF 2989C 2990 111 FORMAT(3X,A40,2X,F9.6) 2991C 2992 DO 110 I = 1,NRHFT 2993C 2994 SUMRHF = SUMRHF + XARR(NORBT - I + 1) 2995C 2996 110 CONTINUE 2997C 2998 TOMOVE = 2*NRHFT - SUMRHF 2999C 3000 WRITE(LUPRI,*) ' ' 3001 WRITE (LUPRI,222) 'Dynamical correlation move:', 3002 * TOMOVE,'electrons' 3003C 3004 222 FORMAT(3X,A27,F9.6,1X,A9) 3005C 3006 CALL QEXIT('CCNAOCAN') 3007C 3008 RETURN 3009 END 3010C /* Deck cc_fcd1ao */ 3011 SUBROUTINE CC_FCD1AO(AODEN,WORK,LWORK,MODEL) 3012C 3013C Written by Ove Christiansen 2-may 1996. 3014C 3015C Version: 1.0 3016C 3017C Purpose: To calculate the Frozen core contribution to the one electron Coupled Cluster 3018C density matrix and add it to AODEN. 3019C 3020#include "implicit.h" 3021#include "priunit.h" 3022#include "dummy.h" 3023#include "maxash.h" 3024#include "maxorb.h" 3025#include "mxcent.h" 3026#include "aovec.h" 3027#include "iratdef.h" 3028 PARAMETER (ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, TWO = 2.0D0) 3029 DIMENSION AODEN(*), WORK(LWORK) 3030#include "ccorb.h" 3031#include "infind.h" 3032#include "inftap.h" 3033#include "blocks.h" 3034#include "ccsdinp.h" 3035#include "ccsdsym.h" 3036#include "ccsdio.h" 3037#include "distcl.h" 3038#include "cbieri.h" 3039#include "eribuf.h" 3040C 3041 CHARACTER MODEL*10 3042C 3043 CALL QENTER('CC_FCD1AO') 3044C 3045 IF (IPRINT .GT. 10) THEN 3046 CALL AROUND( 'Calculate core contribution to density') 3047 ENDIF 3048C 3049 KCMO = 1 3050 KEND = KCMO + NLAMDS 3051 LWRK1 = LWORK - KEND 3052C 3053 IF (LWRK1 .LT. 0) THEN 3054 CALL QUIT('Insufficient spaces in CC_FCD1AO') 3055 ENDIF 3056C 3057C---------------------------------------------- 3058C Read MO-coefficients from interface file. 3059C---------------------------------------------- 3060C 3061 LUSIFC = -1 3062 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY, 3063 & .FALSE.) 3064 REWIND LUSIFC 3065C 3066 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 3067 READ (LUSIFC) 3068C 3069 READ (LUSIFC) 3070 READ (LUSIFC) (WORK(KCMO+I-1), I=1,NLAMDS) 3071C 3072 CALL GPCLOSE(LUSIFC,'KEEP') 3073C 3074 ICMOI = KCMO - 1 3075 DO ISYMI = 1, NSYM 3076 ISYMA = ISYMI 3077 ISYMB = ISYMI 3078 DO IFR = 1, NRHFFR(ISYMI) 3079 I = KFRRHF(IFR,ISYMI) 3080 DO IAL = 1, NBAS(ISYMA) 3081 ICMOIA = ICMOI + NBAS(ISYMA)*(I-1) + IAL 3082 DO IBE = 1, NBAS(ISYMA) 3083 IDAB = IAODIS(ISYMA,ISYMB)+NBAS(ISYMA)*(IBE-1)+IAL 3084 ICMOIB = ICMOI + NBAS(ISYMB)*(I-1) + IBE 3085 AODEN(IDAB) = AODEN(IDAB) + TWO*WORK(ICMOIA)*WORK(ICMOIB) 3086 ENDDO 3087 ENDDO 3088 ENDDO 3089 ICMOI = ICMOI + NBAS(ISYMI)*NRHFS(ISYMI) 3090 ICMOI = ICMOI + NBAS(ISYMI)*NVIRS(ISYMI) 3091 ENDDO 3092C 3093 CALL QEXIT('CC_FCD1AO') 3094C 3095 END 3096C /* Deck cc2_d1ao */ 3097 SUBROUTINE CC2_D1AO(AODEN,ZKAM,WORK,LWORK) 3098C 3099C Written by A. Halkier & S. Coriani 13/01-2000. 3100C 3101C Version: 1.0 3102C 3103C Purpose: To calculate the relaxed one electron CC2 3104C density matrix and return it backtransformed 3105C to AO-basis in AODEN! ZKAM is kappa-bar-0 3106C (input), which contributes to the relaxed density. 3107C 3108C Important note: NO FROZEN CORE SO FAR!!! 3109C 3110#include "implicit.h" 3111#include "priunit.h" 3112#include "dummy.h" 3113#include "ccinftap.h" 3114#include "maxash.h" 3115#include "maxorb.h" 3116#include "mxcent.h" 3117#include "aovec.h" 3118#include "iratdef.h" 3119 PARAMETER (ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, TWO = 2.0D0) 3120 DIMENSION INDEXA(MXCORB_CC) 3121 DIMENSION AODEN(*), ZKAM(*), WORK(LWORK) 3122 CHARACTER*10 MODEL 3123#include "ccorb.h" 3124#include "infind.h" 3125#include "blocks.h" 3126#include "ccsdinp.h" 3127#include "ccsdsym.h" 3128#include "ccsdio.h" 3129#include "distcl.h" 3130#include "cbieri.h" 3131#include "eribuf.h" 3132#include "ccfop.h" 3133#include "ccfro.h" 3134C 3135 LOGICAL LOCDBG 3136 PARAMETER ( LOCDBG=.FALSE. ) 3137C 3138 CALL QENTER('CC2_D1AO') 3139C 3140C------------------------------------------------------------------ 3141C Both t-vectors and tbar-vectors (zeta) are totally symmetric. 3142C------------------------------------------------------------------ 3143C 3144 ISYMTR = 1 3145 ISYMOP = 1 3146C 3147C------------------------------- 3148C Work space allocation one. 3149C------------------------------- 3150C 3151 KT2AM = 1 3152 KZ2AM = KT2AM + NT2AMX 3153 KLAMDP = KZ2AM + NT2SQ(1) 3154 KLAMDH = KLAMDP + NLAMDT 3155 KT1AM = KLAMDH + NLAMDT 3156 KZ1AM = KT1AM + NT1AMX 3157 KEND1 = KZ1AM + NT1AMX 3158 LWRK1 = LWORK - KEND1 3159C 3160 IF (LWRK1 .LT. 0) THEN 3161 WRITE(LUPRI,*) 'Available:', 3162 * LWORK, 'Needed:', KEND1 3163 CALL QUIT( 3164 * 'Insufficient memory for initial allocation in CC2_D1AO') 3165 ENDIF 3166C 3167C---------------------------------------- 3168C Read zero'th order zeta amplitudes. 3169C---------------------------------------- 3170C 3171 IOPT = 3 3172 CALL CC_RDRSP('L0',0,1,IOPT,MODEL,WORK(KZ1AM),WORK(KZ2AM)) 3173C 3174C-------------------------------- 3175C Square up zeta2 amplitudes. 3176C-------------------------------- 3177C 3178 CALL DCOPY(NT2AMX,WORK(KZ2AM),1,WORK(KT2AM),1) 3179 CALL CC_T2SQ(WORK(KT2AM),WORK(KZ2AM),1) 3180C 3181C 3182C------------------------------------------- 3183C Read zero'th order cluster amplitudes. 3184C------------------------------------------- 3185C 3186 IOPT = 3 3187 CALL CC_RDRSP('R0',0,1,IOPT,MODEL,WORK(KT1AM),WORK(KT2AM)) 3188C 3189C---------------------------------- 3190C Calculate the lamda matrices. 3191C---------------------------------- 3192C 3193 CALL LAMMAT(WORK(KLAMDP),WORK(KLAMDH),WORK(KT1AM),WORK(KEND1), 3194 * LWRK1) 3195C 3196C----------------------------------------------- 3197C Set up 2C-E of cluster amplitudes and save 3198C in KT2AM, as we only need T(2c-e) below. 3199C----------------------------------------------- 3200C 3201 ISYOPE = 1 3202 IOPTTCME = 1 3203 CALL CCSD_TCMEPK(WORK(KT2AM),1.0D0,ISYOPE,IOPTTCME) 3204 KT2AMT = KT2AM !for safety 3205C 3206c write(6,*) 'after initial stuff' 3207c CALL FLSHFO(LUPRI) 3208C 3209C------------------------------- 3210C Work space allocation one. 3211C Note that D(ai) = ZETA(ai) 3212C and both D(ia) and h(ia) 3213C are stored transposed! 3214C------------------------------- 3215C 3216 KONEAI = KZ1AM 3217 KONEAB = KONEAI + NT1AMX 3218 KONEIJ = KONEAB + NMATAB(1) 3219 KONEIA = KONEIJ + NMATIJ(1) 3220 KCMO = KONEIA + NT1AMX 3221 KDUM = KCMO + NLAMDS 3222 KEND1 = KDUM + NMATIJ(1) 3223 LWRK1 = LWORK - KEND1 3224C 3225 IF (LWRK1 .LT. 0) THEN 3226 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 3227 CALL QUIT('Insufficient memory for allocation 1 CC2_D1AO') 3228 ENDIF 3229C 3230C 3231C------------------------------------------------------ 3232C Initialize remaining one electron density arrays. 3233C------------------------------------------------------ 3234C 3235 CALL DZERO(WORK(KONEAB),NMATAB(1)) 3236 CALL DZERO(WORK(KONEIJ),NMATIJ(1)) 3237 CALL DZERO(WORK(KONEIA),NT1AMX) 3238C 3239C-------------------------------------------------------- 3240C Construct remaining blocks of one electron density. 3241C-------------------------------------------------------- 3242C 3243 CALL DZERO(WORK(KDUM),NMATIJ(1)) 3244 CALL DIJGEN(WORK(KONEIJ),WORK(KDUM)) 3245 CALL DIAGEN(WORK(KONEIA),WORK(KT2AMT),WORK(KONEAI)) 3246C 3247c write(6,*) 'after D-lambda' 3248c CALL FLSHFO(LUPRI) 3249C 3250 IF (LOCDBG) THEN 3251 DIJNO = DDOT(NMATIJ(1),WORK(KONEIJ),1,WORK(KONEIJ),1) 3252 DAINO = DDOT(NT1AMX,WORK(KONEAI),1,WORK(KONEAI),1) 3253 DIANO = DDOT(NT1AMX,WORK(KONEIA),1,WORK(KONEIA),1) 3254 DABNO = DDOT(NMATAB(1),WORK(KONEAB),1,WORK(KONEAB),1) 3255 WRITE(LUPRI,*) ' ' 3256 WRITE(LUPRI,*) 'Norm of one electron densities in MO-basis:' 3257 WRITE(LUPRI,*) DIJNO, DAINO, DIANO, DABNO 3258 ENDIF 3259C 3260C-------------------------------------------------------- 3261C Backtransform the one electron density to AO-basis. 3262C-------------------------------------------------------- 3263C 3264 CALL DZERO(AODEN,N2BST(1)) 3265C 3266 ISDEN = 1 3267 CALL CC_DENAO(AODEN,1,WORK(KONEAI),WORK(KONEAB), 3268 * WORK(KONEIJ),WORK(KONEIA),ISDEN,WORK(KLAMDP),1, 3269 * WORK(KLAMDH),1,WORK(KEND1),LWRK1) 3270C 3271C---------------------------------------------------------- 3272C Read MO-coefficients from interface file and reorder. 3273C---------------------------------------------------------- 3274C 3275 LUSIFC = -1 3276 CALL GPOPEN(LUSIFC,'SIRIFC','UNKNOWN',' ','UNFORMATTED', 3277 * IDUMMY,.FALSE.) 3278 REWIND LUSIFC 3279 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 3280 READ (LUSIFC) 3281 READ (LUSIFC) 3282 READ (LUSIFC) (WORK(KCMO+I-1), I=1,NLAMDS) 3283 CALL GPCLOSE (LUSIFC,'KEEP') 3284C 3285 CALL CMO_REORDER(WORK(KCMO),WORK(KEND1),LWRK1) 3286C 3287C---------------------------------- 3288C Add orbital relaxation terms. 3289C---------------------------------- 3290C 3291 KOFDIJ = 1 3292 KOFDAB = KOFDIJ + NMATIJ(1) 3293 KOFDAI = KOFDAB + NMATAB(1) 3294 KOFDIA = KOFDAI + NT1AMX 3295C 3296 ISDEN = 1 3297 CALL CC_DENAO(AODEN,1,ZKAM(KOFDAI),ZKAM(KOFDAB), 3298 * ZKAM(KOFDIJ),ZKAM(KOFDIA),ISDEN,WORK(KCMO),1, 3299 * WORK(KCMO),1,WORK(KEND1),LWRK1) 3300C 3301 CALL QEXIT('CC2_D1AO') 3302 RETURN 3303 END 3304