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 ccsdpt_dens2_sc */ 20 SUBROUTINE CCSDPT_DENS2_SC(T1AM,ISYMT1,T2TP,ISYMT2,MODEL, 21 * L1AM,ISYML1,L2TP,ISYML2,WORK,LWORK, 22 * LUDELD,FNDELD,LUCKJD,FNCKJD,LUDKBC,FNDKBC, 23 * LUTOC,FNTOC,LU3VI,FN3VI,LUDKBC3,FNDKBC3, 24 * LU3FOP,FN3FOP,LU3FOP2,FN3FOP2, 25 * LU3FOPX,FN3FOPX,LU3FOP2X,FN3FOP2X) 26C 27C Written by K. Hald, Fall 2001. 28C Modified (temporary) version for CCSD(T) gradients with symmetry. 29C S. Coriani, Fall 2002 30C 31C Calculate the triples contribution to the electronic densities 32C in the MO basis and store them on file. 33C Calculate also the diagonal kappabar multipliers if (RELORB). 34C 35C ISYMT2 is symmetry of T2TP 36C ISYMT1 is symmetry of T1AM 37C Isyres = isymt1*isymt2*isymop 38C 39C For CCSD(T) LUDKBC3, FNDKBC3 is actually LU3VI2, FN3VI2 40C For CCSD(T) we do not use LUDELD,FNDELD,LUCKJD,FNCKJD,LUDKBC,FNDKBC 41C 42 IMPLICIT NONE 43C 44#include "priunit.h" 45#include "dummy.h" 46#include "iratdef.h" 47#include "ccsdsym.h" 48#include "inftap.h" 49#include "ccinftap.h" 50#include "ccorb.h" 51#include "ccsdinp.h" 52#include "ccfop.h" 53#include "second.h" 54C 55 INTEGER ISYMTR, ISYMT1, ISYMT2, ISYML1, ISYML2, LWORK 56 INTEGER ISYRES, ISINT1, ISINT2, ISYMIM, KFOCKD 57 INTEGER KOMG22, KCMO, KEND0, LWRK0, KTROC2 58 INTEGER KTROC0, KXIAJB, KEND1, LWRK1, KINTOC, KEND2, LWRK2 59 INTEGER LENGTH, ISYOPE, IOPTTCME, IOFF, ISYMD, ISAIJ1, ISYCKB 60 INTEGER ISCKB1, ISCKB2, KTRVI1, KTRVI2, KRMAT1, KTRVI0 61 INTEGER KTRVI3, KEND3, LWRK3, KINTVI, KEND4, LWRK4, ISYMB 62 INTEGER ISYALJ, ISAIJ2, ISYMBD, ISCKIJ, KSMAT2, KSMAT, KQMAT 63 INTEGER KDIAG, ISYMC, ISYMK, KOFF1, KOFF2, KOFF3 64 INTEGER KINDSQ, KINDEX, KTMAT, KRMAT2, LENSQ 65 INTEGER LUFCK, KFCKBA, KT2TCME, IOPTT2, KTRVI4, KTRVI5 66 INTEGER KTRVI6, KQMAT2, KVIR1, KVIR2, KVIR3, KVIR4, LUPTIA 67 INTEGER LUPTIAJB, LUABI1, LUABI2, LUABI3, LUABI4, ISYAIB 68 INTEGER ISYMAI, ISYAID, KOCC1, KOCC2, KOMG1, KUMAT, KUMAT2 69 INTEGER LUAIJK, LUIAJK, LUPTAB, LUPTIJ, LUPTIA2 70 INTEGER KTROC02, KTROC22, KTRVI7, KTRVI8, KTRVI9, KTRVI10 71 INTEGER KTRVI11, KTRVI12, KTRVI13, KDENSAB, KDENSIJ 72 INTEGER ISYCKD, ISCKD2, KSMAT3, KUMAT3, KEND5, LWRK5 73 INTEGER KKAPAA, KKAPII, LUKAPAB, LUKAPIJ, ISYALJ2, KINDEX2 74 INTEGER KSMAT4, KUMAT4, KOMG12, KLAMDP, KLAMDH 75 INTEGER KTRVI14, KTRVI15, KTRVI16, KTRVI17, KTRVI18, KTRVI19 76 INTEGER KTRVI20, ISYTMP, KTROC01, KTROC21, KTROC03, KTROC23 77 INTEGER LUDELD, LUCKJD, LUDKBC, LUTOC, LU3VI, LUDKBC3, LU3FOP 78 INTEGER LU3FOP2, LU3FOPX, LU3FOP2X 79 integer KENDSV,NTOT 80C 81#if defined (SYS_CRAY) 82 REAL T1AM(*), T2TP(*) 83 REAL L1AM(*), L2TP(*) 84 REAL WORK(LWORK), ONE 85 REAL TITRAN, TISORT, TISMAT, TIQMAT, TIOME1 86 REAL RHO1N, RHO2N 87 REAL XT2TP, DDOT, XIAJB, XINT, XTROC, XTROC1, XTROC0 88 REAL XTRVI0, XTRVI2, XTRVI3, XTRVI, XTRVI1, XDIA 89 REAL XSMAT, XTMAT, XQMAT, XRMAT, ZERO, TWO, HALF 90 REAL DTIME 91#else 92 DOUBLE PRECISION T1AM(*), T2TP(*) 93 DOUBLE PRECISION L1AM(*), L2TP(*) 94 DOUBLE PRECISION WORK(LWORK), ONE 95 DOUBLE PRECISION TITRAN, TISORT, TISMAT, TIQMAT, TIOME1 96 DOUBLE PRECISION RHO1N, RHO2N 97 DOUBLE PRECISION XT2TP, DDOT, XIAJB, XINT, XTROC, XTROC1, XTROC0 98 DOUBLE PRECISION XTRVI0, XTRVI2, XTRVI3, XTRVI, XTRVI1, XDIA 99 DOUBLE PRECISION XSMAT, XTMAT, XQMAT, XRMAT, ZERO, TWO, HALF 100 DOUBLE PRECISION DTIME 101#endif 102C 103 LOGICAL C3LRSV, CC1ASV, CC1BSV, LDEBUG 104 LOGICAL LOCDBG 105 PARAMETER (LOCDBG = .FALSE.) 106 CHARACTER*(*) FNDELD, FNCKJD, FNDKBC, FNTOC, FN3VI, FNDKBC3 107 CHARACTER*(*) FN3FOP, FN3FOP2 108 CHARACTER*(*) FN3FOPX, FN3FOP2X 109 CHARACTER*5 FNDPTIA, FNDPTAB, FNDPTIJ, FNKAPAB, FNKAPIJ 110 CHARACTER*6 FNDPTIA2 111 CHARACTER*7 FNDIAJB, FNDAIJK, FNDIAJK 112 CHARACTER*8 FNDABI1, FNDABI2, FNDABI3, FNDABI4 113 CHARACTER*10 MODEL 114C 115 PARAMETER (ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, TWO = 2.0D0) 116C 117 CALL QENTER('CCSDPT_DENS2_SC') 118C 119 LDEBUG = .FALSE. 120C 121C------------------------------------------------------------- 122C Set symmetry flags. 123C 124C omega = int1*T2*int2 125C isymres is symmetry of result(omega) 126C isint1 is symmetry of integrals in contraction.(int1) 127C isint2 is symmetry of integrals in the triples equation.(int2) 128C isymim is symmetry of S and Q intermediates.(t2*int2) 129C (sym is for all index of S and Q (cbd,klj) 130C thus cklj=b*d*isymim) 131C------------------------------------------------------------- 132C 133 IPRCC = IPRINT 134 ISYMTR = MULD2H(ISYMT1,ISYMT2) 135 ISYRES = MULD2H(ISYMTR,ISYMOP) 136 ISINT1 = ISYMOP 137 ISINT2 = MULD2H(ISYMT1,ISYMOP) 138 ISYMIM = MULD2H(ISYMTR,ISYMOP) 139C 140C-------------------- 141C Time variables. 142C-------------------- 143C 144 TITRAN = 0.0D0 145 TISORT = 0.0D0 146 TISMAT = 0.0D0 147 TIQMAT = 0.0D0 148 TIOME1 = 0.0D0 149C 150C-------------------------------------- 151C Reorder the t2-amplitudes i T2TP. 152C-------------------------------------- 153C 154 IF (LWORK .LT. NT2SQ(ISYMT2)) THEN 155 CALL QUIT('Not enough memory to construct T2TP (CCSDPT_DENS2)') 156 ENDIF 157C 158 CALL DCOPY(NT2SQ(ISYMT2),T2TP,1,WORK,1) 159 CALL CC3_T2TP(T2TP,WORK,ISYMT2) 160C 161 IF (IPRINT .GT. 55) THEN 162 XT2TP = DDOT(NT2SQ(ISYMT2),T2TP,1,T2TP,1) 163 WRITE(LUPRI,*) 'Norm of T2TP ',XT2TP 164 ENDIF 165C 166C----------------------------------------------- 167C Reorder the l2-amplitudes i L2TP if CC3. 168C----------------------------------------------- 169C 170 IF (CC3) THEN 171C 172 IF (LWORK .LT. NT2SQ(ISYML2)) THEN 173 CALL QUIT('Not enough memory to construct L2TP') 174 ENDIF 175C 176 CALL DCOPY(NT2SQ(ISYML2),L2TP,1,WORK,1) 177 CALL CC3_T2TP(L2TP,WORK,ISYML2) 178C 179 IF (IPRINT .GT. 55) THEN 180 XT2TP = DDOT(NT2SQ(ISYML2),L2TP,1,L2TP,1) 181 WRITE(LUPRI,*) 'Norm of L2TP ',XT2TP 182 ENDIF 183C 184 ENDIF 185C 186C--------------------------------------------------------- 187C Read canonical orbital energies and MO coefficients. 188C--------------------------------------------------------- 189C 190 KFOCKD = 1 191 KOMG1 = KFOCKD + NORBTS 192 KOMG22 = KOMG1 + NT1AM(ISYMOP) 193 KFCKBA = KOMG22 + NT2AM(ISYMOP) 194 KEND0 = KFCKBA + N2BST(ISYMOP) 195C 196 IF (CC3) THEN 197 KLAMDP = KEND0 198 KLAMDH = KLAMDP + NLAMDT 199 KEND0 = KLAMDH + NLAMDT 200 ELSE 201 KCMO = KEND0 202 KEND0 = KCMO + NLAMDS 203 ENDIF 204C 205 LWRK0 = LWORK - KEND0 206C 207 IF (LWRK0 .LT. 0) THEN 208 WRITE(LUPRI,*) 'Memory available : ',LWORK 209 WRITE(LUPRI,*) 'Memory needed : ',KEND0 210 CALL QUIT('Insufficient space in CCSDPT_DENS2') 211 END IF 212C 213 CALL DZERO(WORK(KOMG1),NT1AM(ISYMOP)) 214 CALL DZERO(WORK(KOMG22),NT2AM(ISYMOP)) 215C 216 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY, 217 & .FALSE.) 218 REWIND LUSIFC 219C 220 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 221 READ (LUSIFC) 222 READ (LUSIFC) (WORK(KFOCKD+I-1), I=1,NORBTS) 223C 224 IF (.NOT. CC3) THEN 225 READ (LUSIFC) (WORK(KCMO+I-1), I=1,NLAMDS) 226 CALL CMO_REORDER(WORK(KCMO),WORK(KEND0),LWRK0) 227 ENDIF 228C 229 CALL GPCLOSE(LUSIFC,'KEEP') 230C 231C--------------------------------------------- 232C Delete frozen orbitals in Fock diagonal. 233C--------------------------------------------- 234C 235 IF (FROIMP .OR. FROEXP) 236 * CALL CCSD_DELFRO(WORK(KFOCKD),WORK(KEND0),LWRK0) 237C 238C---------------------------------------------- 239C Calculate the lambda matrices for CC3 240C---------------------------------------------- 241C 242 IF (CC3) THEN 243 CALL LAMMAT(WORK(KLAMDP),WORK(KLAMDH),T1AM, 244 * WORK(KEND0),LWRK0) 245C 246 IF (IPRINT .GT.100) THEN 247 CALL AROUND('Usual Lambda matrices ') 248 CALL CC_PRLAM(WORK(KLAMDP),WORK(KLAMDH),1) 249 ENDIF 250 ENDIF 251C 252C----------------------------------------------------- 253C Construct the transformed Fock matrix 254C----------------------------------------------------- 255C 256 LUFCK = -1 257C 258 IF (CC3) THEN 259C This AO Fock matrix is constructed from the T1 transformed density 260 CALL GPOPEN(LUFCK,'CC_FCKH','UNKNOWN',' ','UNFORMATTED', 261 * IDUMMY,.FALSE.) 262 ELSE 263C This AO Fock matrix is constructed from the CMO transformed density 264 CALL GPOPEN(LUFCK,'CC_FCKREF','UNKNOWN',' ','UNFORMATTED', 265 * IDUMMY,.FALSE.) 266 ENDIF 267C 268 REWIND(LUFCK) 269 READ(LUFCK)(WORK(KFCKBA + I-1),I = 1,N2BST(ISYMOP)) 270 CALL GPCLOSE(LUFCK,'KEEP' ) 271C 272 IF (IPRINT .GT. 140) THEN 273 CALL AROUND( 'Usual Fock AO matrix' ) 274 CALL CC_PRFCKAO(WORK(KFCKBA),ISYMOP) 275 ENDIF 276C 277 ! SCF Fock matrix in transformed using CMO vector 278 IF (CC3) THEN 279 CALL CC_FCKMO(WORK(KFCKBA),WORK(KLAMDP),WORK(KLAMDH), 280 * WORK(KEND0),LWRK0,1,1,1) 281 ELSE 282 CALL CC_FCKMO(WORK(KFCKBA),WORK(KCMO),WORK(KCMO), 283 * WORK(KEND0),LWRK0,1,1,1) 284 ENDIF 285C 286 IF (IPRINT .GT. 50) THEN 287 CALL AROUND( 'In CC_ETA: Triples Fock MO matrix' ) 288 CALL CC_PRFCKMO(WORK(KFCKBA),ISYMOP) 289 ENDIF 290C 291C Sort the fock matrix 292C 293C 294 CALL DCOPY(N2BST(ISINT1),WORK(KFCKBA),1,WORK(KEND0),1) 295C 296 DO ISYMC = 1,NSYM 297C 298 ISYMK = MULD2H(ISYMC,ISINT1) 299C 300 DO K = 1,NRHF(ISYMK) 301C 302 DO C = 1,NVIR(ISYMC) 303C 304 KOFF1 = KEND0 + IFCVIR(ISYMK,ISYMC) + 305 * NORB(ISYMK)*(C - 1) + K - 1 306 KOFF2 = KFCKBA + IT1AM(ISYMC,ISYMK) 307 * + NVIR(ISYMC)*(K - 1) + C - 1 308C 309 WORK(KOFF2) = WORK(KOFF1) 310C 311 ENDDO 312 ENDDO 313 ENDDO 314C 315 IF (IPRINT .GT. 50) THEN 316 CALL AROUND('In CCSDPT_DENS2: Triples Fock MO matrix (sort)') 317 CALL CC_PRFCKMO(WORK(KFCKBA),ISYMOP) 318 ENDIF 319C 320C------------------------------------------------------------------ 321C Read in another T2 amplitude, and transform it to 2*C-E 322C Square up to full matrix and reorder the index 323C------------------------------------------------------------------ 324C 325 IF (.NOT. CC3) THEN 326 KT2TCME = KEND0 327 KEND0 = KT2TCME + NT2SQ(1) 328 LWRK0 = LWORK - KEND0 329C 330 IF (LWRK0 .LT. NT2SQ(1)) 331 * CALL QUIT('Too litlle workspace CCSDPT_DENS2 T2TCME') 332C 333 IOPTT2 = 2 334 CALL CC_RDRSP('R0',0,1,IOPTT2,MODEL,DUMMY,WORK(KEND0)) 335C 336 ISYOPE = ISYMOP 337 IOPTT2 = 1 338 CALL CCSD_TCMEPK(WORK(KEND0),1.0D0,ISYOPE,IOPTT2) 339C 340 CALL CC_T2SQ(WORK(KEND0),WORK(KT2TCME),1) 341C 342 CALL DCOPY(NT2SQ(1),WORK(KT2TCME),1,WORK(KEND0),1) 343 CALL CC3_T2TP(WORK(KT2TCME),WORK(KEND0),1) 344C 345 IF (IPRINT .GT. 55) THEN 346 XT2TP = DDOT(NT2SQ(1),WORK(KT2TCME),1,WORK(KT2TCME),1) 347 WRITE(LUPRI,*) 'Norm of 2*C-E T2 amplitudes after resort ', 348 * XT2TP 349 ENDIF 350 ENDIF 351C 352C----------------------------- 353C Read occupied integrals. 354C----------------------------- 355C 356C Memory allocation. 357C 358 KTROC0 = KEND0 359 KTROC02= KTROC0 + NTRAOC(ISINT2) 360 KTROC2 = KTROC02+ NTRAOC(ISINT2) 361 KTROC22= KTROC2 + NTRAOC(ISINT2) 362 KXIAJB = KTROC22+ NTRAOC(ISINT2) 363 KOCC1 = KXIAJB + NT2AM(ISYMOP) 364 KOCC2 = KOCC1 + NCKIJ(ISYRES) 365 KKAPAA = KOCC2 + NCKIJ(ISYRES) 366 KKAPII = KKAPAA + NVIRT 367 KEND1 = KKAPII + NRHFT 368 LWRK1 = LWORK - KEND1 369C 370 IF (CC3) THEN 371 KTROC01 = KEND1 372 KTROC21 = KTROC01 + NTRAOC(ISINT2) 373 KTROC03 = KTROC21 + NTRAOC(ISINT2) 374 KTROC23 = KTROC03 + NTRAOC(ISINT2) 375 KEND1 = KTROC23 + NTRAOC(ISINT2) 376 LWRK1 = LWORK - KEND1 377 ENDIF 378C 379 IF (.NOT. RELORB) THEN 380 KOMG12 = KEND1 381 KDENSAB = KOMG12 + NT1AM(ISYRES) 382 KDENSIJ = KDENSAB + NMATAB(ISYRES) 383 KEND1 = KDENSIJ + NMATIJ(ISYRES) 384 LWRK1 = LWORK - KEND1 385C 386 CALL DZERO(WORK(KOMG12),NT1AM(ISYRES)) 387 CALL DZERO(WORK(KDENSAB),NMATAB(ISYRES)) 388 CALL DZERO(WORK(KDENSIJ),NMATIJ(ISYRES)) 389 ENDIF 390C 391 KINTOC = KEND1 392 KEND2 = KINTOC + MAX(NTOTOC(ISYMOP),NTOTOC(ISINT2)) 393 LWRK2 = LWORK - KEND2 394C 395 IF (LWRK2 .LT. 0) THEN 396 WRITE(LUPRI,*) 'Memory available : ',LWORK 397 WRITE(LUPRI,*) 'Memory needed : ',KEND2 398 CALL QUIT('Insufficient space in CCSDPT_DENS2') 399 END IF 400C 401C 402C---------------------------------- 403C Initialize result vectors 404C---------------------------------- 405C 406 CALL DZERO(WORK(KOCC1),NCKIJ(ISYRES)) 407 CALL DZERO(WORK(KOCC2),NCKIJ(ISYRES)) 408 CALL DZERO(WORK(KKAPAA),NVIRT) 409 CALL DZERO(WORK(KKAPII),NRHFT) 410C 411C------------------------ 412C Construct L(ia,jb). 413C------------------------ 414C 415 LENGTH = IRAT*NT2AM(ISYMOP) 416C 417 REWIND(LUIAJB) 418 CALL READI(LUIAJB,LENGTH,WORK(KXIAJB)) 419C 420 ISYOPE = ISYMOP 421 IOPTTCME = 1 422 CALL CCSD_TCMEPK(WORK(KXIAJB),1.0D0,ISYOPE,IOPTTCME) 423C 424 IF ( IPRINT .GT. 55) THEN 425 XIAJB = DDOT(NT2AM(ISYMOP),WORK(KXIAJB),1, 426 * WORK(KXIAJB),1) 427 WRITE(LUPRI,*) 'Norm of IAJB ',XIAJB 428 ENDIF 429C 430C------------------------ 431C Occupied integrals. 432C------------------------ 433C 434 IF (CC3) THEN 435 IOFF = 1 436 IF (NTOTOC(ISYMOP) .GT. 0) THEN 437 CALL GETWA2(LUCKJD,FNCKJD,WORK(KINTOC),IOFF,NTOTOC(ISYMOP)) 438 ENDIF 439 ELSE 440 IOFF = 1 441 IF (NTOTOC(ISYMOP) .GT. 0) THEN 442 CALL GETWA2(LUTOC,FNTOC,WORK(KINTOC),IOFF,NTOTOC(ISYMOP)) 443 ENDIF 444 ENDIF 445C 446C---------------------------------- 447C Write out norms of Integrals. 448C---------------------------------- 449C 450 IF (IPRINT .GT. 55) THEN 451 XINT = DDOT(NTOTOC(ISYMOP),WORK(KINTOC),1, 452 * WORK(KINTOC),1) 453 WRITE(LUPRI,*) 'Norm of OCC-INT ',XINT 454 ENDIF 455C 456C---------------------------------------------------------------------- 457C Transform (ia|j delta) integrals to (ia|j k) and sort as (ij,k,a) 458C---------------------------------------------------------------------- 459C 460 DTIME = SECOND() 461 IF (CC3) THEN 462 CALL CC3_TROCC(WORK(KINTOC),WORK(KTROC0),WORK(KLAMDP), 463 * WORK(KEND2),LWRK2,ISINT2) 464 ELSE 465 CALL CCSDT_TROCC(WORK(KINTOC),WORK(KTROC0),WORK(KCMO), 466 * WORK(KEND2),LWRK2) 467 ENDIF 468C 469 DTIME = SECOND() - DTIME 470 TITRAN = TITRAN + DTIME 471 472C 473 DTIME = SECOND() 474C 475 DTIME = SECOND() - DTIME 476 TISORT = TISORT + DTIME 477C 478C----------------------------------------------------------- 479C Construct 2*C-E of the integrals. 480C Have integral for both (ij,k,a) and (a,k,j,i) 481C----------------------------------------------------------- 482C 483 CALL CCSDT_TCMEOCC(WORK(KTROC0),WORK(KTROC2),ISINT2) 484C 485 IF (CC3) THEN 486 IOFF = 1 487 IF (NTOTOC(ISINT2) .GT. 0) THEN 488 CALL GETWA2(LUTOC,FNTOC,WORK(KINTOC),IOFF,NTOTOC(ISINT2)) 489 ENDIF 490C 491 CALL CC3_TROCC(WORK(KINTOC),WORK(KTROC01),WORK(KLAMDH), 492 * WORK(KEND2),LWRK2,ISINT2) 493C 494 CALL CCSDT_TCMEOCC(WORK(KTROC01),WORK(KTROC21),ISINT2) 495C 496 CALL CCFOP_SORT(WORK(KTROC01),WORK(KTROC03),ISINT2,1) 497C 498 CALL CCFOP_SORT(WORK(KTROC21),WORK(KTROC23),ISINT2,1) 499 ENDIF 500C 501 CALL CCFOP_SORT(WORK(KTROC0),WORK(KTROC02),ISINT2,1) 502C 503 CALL CCFOP_SORT(WORK(KTROC2),WORK(KTROC22),ISINT2,1) 504C 505C------------------------------- 506C Write out norms of arrays. 507C------------------------------- 508C 509 IF (IPRINT .GT. 55) THEN 510 XINT = DDOT(NTOTOC(ISINT2),WORK(KINTOC),1, 511 * WORK(KINTOC),1) 512 WRITE(LUPRI,*) 'Norm of CKJDEL-INT ',XINT 513 ENDIF 514C 515 IF (IPRINT .GT. 55) THEN 516 XTROC0 = DDOT(NTRAOC(ISINT2),WORK(KTROC0),1, 517 * WORK(KTROC0),1) 518 WRITE(LUPRI,*) 'Norm of TROC0 ',XTROC0 519 ENDIF 520C 521 IF (IPRINT .GT. 55) THEN 522 XTROC0 = DDOT(NTRAOC(ISINT2),WORK(KTROC2),1, 523 * WORK(KTROC2),1) 524 WRITE(LUPRI,*) 'Norm of TROC2 ',XTROC0 525 ENDIF 526C 527C-------------------------------------------------------- 528C Open files to the one and two electron densities. 529C-------------------------------------------------------- 530C 531 LUPTIA = -1 532 FNDPTIA = 'DPTIA' 533C d_{ia} 534 CALL WOPEN2(LUPTIA,FNDPTIA,64,0) 535C 536 IF ((.NOT. CC3) .AND. (RELORB)) THEN 537 LUPTIAJB = -1 538 LUABI1 = -1 539 LUABI2 = -1 540 LUABI3 = -1 541 LUABI4 = -1 542 LUAIJK = -1 543 LUIAJK = -1 544 FNDIAJB = 'DPTIAJB' 545 FNDABI1 = 'DPTABIC1' 546 FNDABI2 = 'DPTABIC2' 547 FNDABI3 = 'DPTABCI1' 548 FNDABI4 = 'DPTABCI2' 549 FNDAIJK = 'DPTAIJK' 550 FNDIAJK = 'DPTIAJK' 551C 552C d_{iajb} 553 CALL WOPEN2(LUPTIAJB,FNDIAJB,64,0) 554C d_{abic_1} 555 CALL WOPEN2(LUABI1,FNDABI1,64,0) 556C d_{abic_2} 557 CALL WOPEN2(LUABI2,FNDABI2,64,0) 558C d_{abci_1} 559 CALL WOPEN2(LUABI3,FNDABI3,64,0) 560C d_{abci_2} 561 CALL WOPEN2(LUABI4,FNDABI4,64,0) 562C d_{aijk} 563 CALL WOPEN2(LUAIJK,FNDAIJK,64,0) 564C d_{iajk} 565 CALL WOPEN2(LUIAJK,FNDIAJK,64,0) 566 ELSE 567 LUPTIA2 = -1 568 LUPTAB = -1 569 LUPTIJ = -1 570 FNDPTIA2 = 'DPTIA2' 571 FNDPTAB = 'DPTAB' 572 FNDPTIJ = 'DPTIJ' 573C d_{ia} 574 CALL WOPEN2(LUPTIA2,FNDPTIA2,64,0) 575C d_{ab} 576 CALL WOPEN2(LUPTAB,FNDPTAB,64,0) 577C d_{ij} 578 CALL WOPEN2(LUPTIJ,FNDPTIJ,64,0) 579 ENDIF 580C 581C----------------------------------------- 582! Sonia: old code!!!! 583C Save workspace adress to regain workspace 584C in the end. 585C----------------------------------------- 586C 587 KENDSV = KEND1 588C 589C 590C---------------------------- 591C General loop structure. 592C---------------------------- 593C 594 DO ISYMD = 1,NSYM 595C 596 ISAIJ1 = MULD2H(ISYMD,ISYRES) 597 ISYCKB = MULD2H(ISYMD,ISYMOP) 598 ISCKB1 = MULD2H(ISINT1,ISYMD) 599 ISCKB2 = MULD2H(ISINT2,ISYMD) 600C 601 IF (IPRINT .GT. 55) THEN 602C 603 WRITE(LUPRI,*) 'In CCSDPT_DENS2: ISAIJ1 :',ISAIJ1 604 WRITE(LUPRI,*) 'In CCSDPT_DENS2: ISYCKB :',ISYCKB 605 WRITE(LUPRI,*) 'In CCSDPT_DENS2: ISCKB1 :',ISCKB1 606 WRITE(LUPRI,*) 'In CCSDPT_DENS2: ISCKB2 :',ISCKB2 607C 608 ENDIF 609C 610C-------------------------- 611C Memory allocation. 612C-------------------------- 613C 614 KTRVI1 = KEND1 615 KTRVI2 = KTRVI1 + NCKATR(ISCKB2) 616 KRMAT1 = KTRVI2 + NCKATR(ISCKB2) 617 KEND2 = KRMAT1 + NCKI(ISAIJ1) 618 LWRK2 = LWORK - KEND2 619C 620 KTRVI0 = KEND2 621 KTRVI3 = KTRVI0 + NCKATR(ISCKB2) 622 KTRVI4 = KTRVI3 + NCKATR(ISCKB2) 623 KTRVI5 = KTRVI4 + NCKATR(ISCKB2) 624 KTRVI6 = KTRVI5 + NCKATR(ISCKB2) 625 KTRVI7 = KTRVI6 + NCKATR(ISCKB2) 626 KVIR1 = KTRVI7 + NCKATR(ISCKB2) 627 KVIR2 = KVIR1 + NCKATR(ISAIJ1) 628 KVIR3 = KVIR2 + NCKATR(ISAIJ1) 629 KVIR4 = KVIR3 + NCKATR(ISAIJ1) 630 KEND3 = KVIR4 + NCKATR(ISAIJ1) 631 LWRK3 = LWORK - KEND3 632C 633 IF (CC3) THEN 634 KTRVI14 = KEND3 635 KTRVI15 = KTRVI14 + NCKATR(ISCKB2) 636 KTRVI18 = KTRVI15 + NCKATR(ISCKB2) 637 KTRVI19 = KTRVI18 + NCKATR(ISCKB2) 638 KEND3 = KTRVI19 + NCKATR(ISCKB2) 639 LWRK3 = LWORK - KEND3 640 ENDIF 641C 642 KINTVI = KEND3 643 KEND4 = KINTVI + MAX(NCKA(ISYMD),NCKA(ISCKB2)) 644 LWRK4 = LWORK - KEND4 645C 646 IF (LWRK4 .LT. 0) THEN 647 WRITE(LUPRI,*) 'Memory available : ',LWORK 648 WRITE(LUPRI,*) 'Memory needed : ',KEND4 649 CALL QUIT('Insufficient space in CCSDPT_DENS2') 650 END IF 651C 652C--------------------- 653C Sum over D 654C--------------------- 655C 656 DO D = 1,NVIR(ISYMD) 657C 658C------------------------------------ 659C Initialize the R1 matrix. 660C------------------------------------ 661C 662 CALL DZERO(WORK(KRMAT1),NCKI(ISAIJ1)) 663 CALL DZERO(WORK(KVIR1),NCKATR(ISAIJ1)) 664 CALL DZERO(WORK(KVIR2),NCKATR(ISAIJ1)) 665 CALL DZERO(WORK(KVIR3),NCKATR(ISAIJ1)) 666 CALL DZERO(WORK(KVIR4),NCKATR(ISAIJ1)) 667C 668C----------------------------------------------- 669C Integrals used in s3am. 670C----------------------------------------------- 671C 672 IF (CC3) THEN 673 IOFF = ICKBD(ISCKB2,ISYMD) + NCKATR(ISCKB2)*(D - 1) + 1 674 IF (NCKATR(ISCKB2) .GT. 0) THEN 675 CALL GETWA2(LUDKBC,FNDKBC,WORK(KTRVI0),IOFF, 676 & NCKATR(ISCKB2)) 677 ENDIF 678 ELSE 679 IOFF = ICKAD(ISYCKB,ISYMD) + NCKA(ISYCKB)*(D - 1) + 1 680 IF (NCKA(ISYCKB) .GT. 0) THEN 681 CALL GETWA2(LUDKBC3,FNDKBC3,WORK(KINTVI),IOFF, 682 & NCKA(ISYCKB)) 683 ENDIF 684C 685 CALL CCSDT_TRVIR(WORK(KINTVI),WORK(KTRVI0),WORK(KCMO), 686 * ISYMD,D,ISINT2,WORK(KEND4),LWRK4) 687 ENDIF 688C 689C------------------------------------------------------ 690C Read 2*C-E of integral used for t3-bar 691C------------------------------------------------------ 692C 693 IF (CC3) THEN 694 IOFF = ICKBD(ISCKB2,ISYMD) + NCKATR(ISCKB2)*(D - 1) + 1 695 IF (NCKATR(ISCKB2) .GT. 0) THEN 696 CALL GETWA2(LU3FOP2X,FN3FOP2X,WORK(KTRVI4),IOFF, 697 & NCKATR(ISCKB2)) 698 ENDIF 699 ELSE 700 IOFF = ICKAD(ISYCKB,ISYMD) + NCKA(ISYCKB)*(D - 1) + 1 701 IF (NCKA(ISYCKB) .GT. 0) THEN 702 CALL GETWA2(LU3FOP2,FN3FOP2,WORK(KINTVI),IOFF, 703 * NCKA(ISYCKB)) 704 ENDIF 705C 706 CALL CCSDT_TRVIR(WORK(KINTVI),WORK(KTRVI4),WORK(KCMO), 707 * ISYMD,D,ISYMOP,WORK(KEND4),LWRK4) 708 ENDIF 709C 710C------------------------------------------------------------ 711C Integrals used for t3-bar for cc3 712C------------------------------------------------------------ 713C 714 IF (CC3) THEN 715 IOFF = ICKBD(ISCKB2,ISYMD) + NCKATR(ISCKB2)*(D - 1) + 1 716 IF (NCKATR(ISCKB2) .GT. 0) THEN 717 CALL GETWA2(LUDKBC3,FNDKBC3,WORK(KTRVI14),IOFF, 718 & NCKATR(ISCKB2)) 719 ENDIF 720 CALL CCSDT_SRVIR3(WORK(KTRVI14),WORK(KEND4), 721 * ISYMD,D,ISINT2) 722 CALL CCSDT_SRTVIR(WORK(KTRVI14),WORK(KTRVI15),WORK(KEND4) 723 * ,LWRK4,ISYMD,ISINT2) 724 ENDIF 725C 726C----------------------------------------------------------- 727C Sort the integrals for s3am and for t3-bar 728C----------------------------------------------------------- 729C 730 DTIME = SECOND() 731 CALL CCSDT_SRTVIR(WORK(KTRVI0),WORK(KTRVI2),WORK(KEND4), 732 * LWRK4,ISYMD,ISINT2) 733C 734 CALL CCSDT_SRTVIR(WORK(KTRVI4),WORK(KTRVI5),WORK(KEND4), 735 * LWRK4,ISYMD,ISINT2) 736C 737 DTIME = SECOND() - DTIME 738 TISORT = TISORT + DTIME 739C 740 IF (IPRINT .GT. 55) THEN 741 XTRVI0= DDOT(NCKATR(ISCKB2),WORK(KTRVI0),1, 742 * WORK(KTRVI0),1) 743 WRITE(LUPRI,*) 'Norm of TRVI0 ',XTRVI0 744 ENDIF 745C 746 IF (IPRINT .GT. 55) THEN 747 XTRVI2= DDOT(NCKATR(ISCKB2),WORK(KTRVI2),1, 748 * WORK(KTRVI2),1) 749 WRITE(LUPRI,*) 'Norm of TRVI2 ',XTRVI2 750 ENDIF 751C 752 IF (IPRINT .GT. 55) THEN 753 XTRVI0= DDOT(NCKATR(ISCKB2),WORK(KTRVI4),1, 754 * WORK(KTRVI4),1) 755 WRITE(LUPRI,*) 'Norm of TRVI4 ',XTRVI0 756 ENDIF 757C 758 IF (IPRINT .GT. 55) THEN 759 XTRVI2= DDOT(NCKATR(ISCKB2),WORK(KTRVI5),1, 760 * WORK(KTRVI5),1) 761 WRITE(LUPRI,*) 'Norm of TRVI5 ',XTRVI2 762 ENDIF 763C 764C------------------------------------------------------ 765C Read virtual integrals used in contraction. 766C------------------------------------------------------ 767C 768 IF (CC3) THEN 769 IOFF = ICKAD(ISCKB2,ISYMD) + NCKA(ISCKB2)*(D - 1) + 1 770 IF (NCKA(ISCKB2) .GT. 0) THEN 771 CALL GETWA2(LUDELD,FNDELD,WORK(KINTVI),IOFF, 772 * NCKA(ISCKB2)) 773 ENDIF 774C 775 CALL CCSDT_TRVIR(WORK(KINTVI),WORK(KTRVI1),WORK(KLAMDH), 776 * ISYMD,D,ISINT2,WORK(KEND4),LWRK4) 777C 778 ELSE 779 IOFF = ICKAD(ISYCKB,ISYMD) + NCKA(ISYCKB)*(D - 1) + 1 780 IF (NCKA(ISYCKB) .GT. 0) THEN 781 CALL GETWA2(LU3VI,FN3VI,WORK(KINTVI),IOFF, 782 & NCKA(ISYCKB)) 783 ENDIF 784C 785 CALL CCSDT_TRVIR(WORK(KINTVI),WORK(KTRVI1),WORK(KCMO), 786 * ISYMD,D,ISINT2,WORK(KEND4),LWRK4) 787 ENDIF 788C 789C-------------------------------------------------------- 790C Calculate virtual integrals used in q3am. 791C-------------------------------------------------------- 792C 793 CALL DCOPY(NCKATR(ISCKB2),WORK(KTRVI1),1,WORK(KTRVI3),1) 794C 795 IF (LWRK3 .LT. NCKATR(ISYCKB)) THEN 796 CALL QUIT('Insufficient space for allocation in '// 797 & 'CCSDPT_DENS2 (1)') 798 END IF 799C 800 DTIME = SECOND() 801 CALL CCSDT_SRVIR3(WORK(KTRVI3),WORK(KEND4),ISYMD,D,ISINT2) 802C 803 DTIME = SECOND() - DTIME 804 TISORT = TISORT + DTIME 805C 806 IF (IPRINT .GT. 55) THEN 807 XTRVI3= DDOT(NCKATR(ISCKB2),WORK(KTRVI3),1, 808 * WORK(KTRVI3),1) 809 WRITE(LUPRI,*) 'Norm of TRVI3 ',XTRVI3 810 ENDIF 811C 812C--------------------------------------------------------------- 813C Read virtual integrals used in q3am/u3am for t3-bar. 814C--------------------------------------------------------------- 815C 816 IF (CC3) THEN 817 IOFF = ICKAD(ISYCKB,ISYMD) + NCKA(ISYCKB)*(D - 1) + 1 818 IF (NCKA(ISYCKB) .GT. 0) THEN 819 CALL GETWA2(LU3VI,FN3VI,WORK(KINTVI),IOFF, 820 * NCKA(ISYCKB)) 821 ENDIF 822C 823 CALL CCSDT_TRVIR(WORK(KINTVI),WORK(KTRVI19),WORK(KLAMDP), 824 * ISYMD,D,ISYMOP,WORK(KEND4),LWRK4) 825C 826 IF (LWRK4 .LT. NCKATR(ISYCKB)) THEN 827 CALL QUIT('Insufficient space for allocation in '// 828 * 'CCSDPT_DENS2 (CC3 TRVI)') 829 END IF 830C 831 CALL CCSDT_SRTVIR(WORK(KTRVI19),WORK(KTRVI18),WORK(KEND4) 832 * ,LWRK4,ISYMD,ISINT2) 833 ENDIF 834C 835 IF (CC3) THEN 836 IOFF = ICKBD(ISYCKB,ISYMD) + NCKATR(ISYCKB)*(D - 1) + 1 837 IF (NCKATR(ISYCKB) .GT. 0) THEN 838 CALL GETWA2(LU3FOPX,FN3FOPX,WORK(KTRVI6),IOFF, 839 * NCKATR(ISYCKB)) 840 ENDIF 841 ELSE 842 IOFF = ICKAD(ISYCKB,ISYMD) + NCKA(ISYCKB)*(D - 1) + 1 843 IF (NCKA(ISYCKB) .GT. 0) THEN 844 CALL GETWA2(LU3FOP,FN3FOP,WORK(KINTVI),IOFF, 845 * NCKA(ISYCKB)) 846 ENDIF 847C 848 CALL CCSDT_TRVIR(WORK(KINTVI),WORK(KTRVI6),WORK(KCMO), 849 * ISYMD,D,ISINT2,WORK(KEND4),LWRK4) 850 ENDIF 851C 852 IF (LWRK3 .LT. NCKATR(ISYCKB)) THEN 853 CALL QUIT('Insufficient space for allocation in '// 854 & 'CCSDPT_DENS2 (2)') 855 END IF 856C 857 CALL DCOPY(NCKATR(ISCKB2),WORK(KTRVI6),1,WORK(KTRVI7),1) 858C 859 DTIME = SECOND() 860 CALL CCSDT_SRVIR3(WORK(KTRVI6),WORK(KEND4),ISYMD,D,ISINT2) 861C 862 DTIME = SECOND() - DTIME 863 TISORT = TISORT + DTIME 864C 865 IF (IPRINT .GT. 55) THEN 866 XTRVI3= DDOT(NCKATR(ISCKB2),WORK(KTRVI6),1, 867 * WORK(KTRVI6),1) 868 WRITE(LUPRI,*) 'Norm of TRVI6 ',XTRVI3 869 XTRVI3= DDOT(NCKATR(ISCKB2),WORK(KTRVI7),1, 870 * WORK(KTRVI7),1) 871 WRITE(LUPRI,*) 'Norm of TRVI7 ',XTRVI3 872 ENDIF 873C 874 IF (IPRINT .GT. 55) THEN 875 XTRVI1= DDOT(NCKATR(ISCKB2),WORK(KTRVI1),1, 876 * WORK(KTRVI1),1) 877 WRITE(LUPRI,*) 'Norm of TRVI1 ',XTRVI1 878 ENDIF 879C 880C--------------------- 881C Calculate. 882C--------------------- 883C 884 DO ISYMB = 1,NSYM 885C 886 ISYALJ = MULD2H(ISYMB,ISYMT2) 887 ISYALJ2 = MULD2H(ISYMD,ISYMT2) 888 ISAIJ2 = MULD2H(ISYMB,ISYRES) 889 ISYMBD = MULD2H(ISYMB,ISYMD) 890 ISCKIJ = MULD2H(ISYMBD,ISYMIM) 891 ISYCKD = MULD2H(ISYMOP,ISYMB) 892 ISCKD2 = MULD2H(ISINT2,ISYMB) 893C 894 IF ((IPRINT .GT. 55)) THEN 895C 896 WRITE(LUPRI,*) 'In CCSDPT_DENS2: ISYMD :',ISYMD 897 WRITE(LUPRI,*) 'In CCSDPT_DENS2: ISYMB :',ISYMB 898 WRITE(LUPRI,*) 'In CCSDPT_DENS2: ISYALJ:',ISYALJ 899 WRITE(LUPRI,*) 'In CCSDPT_DENS2: ISAIJ2:',ISAIJ2 900 WRITE(LUPRI,*) 'In CCSDPT_DENS2: ISYMBD:',ISYMBD 901 WRITE(LUPRI,*) 'In CCSDPT_DENS2: ISCKIJ:',ISCKIJ 902C 903 ENDIF 904C 905C Can use kend3 since we do not need the integrals anymore. 906 KSMAT = KEND3 907 KQMAT = KSMAT + NCKIJ(ISCKIJ) 908 KSMAT2 = KQMAT + NCKIJ(ISCKIJ) 909 KSMAT3 = KSMAT2 + NCKIJ(ISCKIJ) 910 KQMAT2 = KSMAT3 + NCKIJ(ISCKIJ) 911 KUMAT = KQMAT2 + NCKIJ(ISCKIJ) 912 KUMAT2 = KUMAT + NCKIJ(ISCKIJ) 913 KUMAT3 = KUMAT2 + NCKIJ(ISCKIJ) 914 KDIAG = KUMAT3 + NCKIJ(ISCKIJ) 915 KINDSQ = KDIAG + NCKIJ(ISCKIJ) 916 KINDEX = KINDSQ + (6*NCKIJ(ISCKIJ) - 1)/IRAT + 1 917 KINDEX2 = KINDEX + (NCKI(ISYALJ) - 1)/IRAT + 1 918 KTMAT = KINDEX2 + (NCKI(ISYALJ2) - 1)/IRAT + 1 919 KRMAT2 = KTMAT + NCKIJ(ISCKIJ) 920 KTRVI8 = KRMAT2 + NCKI(ISAIJ2) 921 KTRVI9 = KTRVI8 + NCKATR(ISCKD2) 922 KTRVI10 = KTRVI9 + NCKATR(ISCKD2) 923 KEND4 = KTRVI10 + NCKATR(ISCKD2) 924 LWRK4 = LWORK - KEND4 925C 926 IF (CC3) THEN 927 KTRVI16 = KEND4 928 KTRVI17 = KTRVI16 + NCKATR(ISCKD2) 929 KTRVI20 = KTRVI17 + NCKATR(ISCKD2) 930 KEND4 = KTRVI20 + NCKATR(ISCKD2) 931 LWRK4 = LWORK - KEND4 932 ENDIF 933C 934 IF (.NOT. RELORB) THEN 935 KSMAT4 = KEND4 936 KUMAT4 = KSMAT4 + NCKIJ(ISCKIJ) 937 KTRVI11 = KUMAT4 + NCKIJ(ISCKIJ) 938 KTRVI12 = KTRVI11 + NCKATR(ISCKD2) 939 KTRVI13 = KTRVI12 + NCKATR(ISCKD2) 940 KEND4 = KTRVI13 + NCKATR(ISCKD2) 941 LWRK4 = LWORK-KEND4 942 ENDIF 943C 944 KINTVI = KEND4 945COMMENT COMMENT 946C KEND5 = KINTVI + NCKA(ISCKD2) 947 KEND5 = KINTVI + MAX(NCKA(ISYMB),NCKA(ISCKD2)) 948COMMENT COMMENT 949 LWRK5 = LWORK - KEND5 950C 951 IF (LWRK5 .LT. 0) THEN 952 WRITE(LUPRI,*) 'Memory available : ',LWORK 953 WRITE(LUPRI,*) 'Memory needed : ',KEND5 954 CALL QUIT('Insufficient space in CCSDPT_DENS2') 955 END IF 956C 957C--------------------------------------------- 958C Construct part of the diagonal. 959C--------------------------------------------- 960C 961 CALL CC3_DIAG(WORK(KDIAG),WORK(KFOCKD),ISCKIJ) 962C 963 IF ((IPRINT .GT. 55)) THEN 964 XDIA = DDOT(NCKIJ(ISCKIJ),WORK(KDIAG),1, 965 * WORK(KDIAG),1) 966 WRITE(LUPRI,*) 'Norm of DIA ',XDIA 967 ENDIF 968 969C 970C------------------------------------- 971C Construct index arrays. 972C------------------------------------- 973C 974 LENSQ = NCKIJ(ISCKIJ) 975 CALL CC3_INDSQ(WORK(KINDSQ),LENSQ,ISCKIJ) 976 CALL CC3_INDEX(WORK(KINDEX),ISYALJ) 977 CALL CC3_INDEX(WORK(KINDEX2),ISYALJ2) 978C 979 DO B = 1,NVIR(ISYMB) 980C 981C----------------------------------------- 982C Initialize the R2 matrix. 983C----------------------------------------- 984C 985 CALL DZERO(WORK(KRMAT2),NCKI(ISAIJ2)) 986C 987C------------------------------------------------------------- 988C Read and transform integrals used in second S 989C------------------------------------------------------------- 990C 991 IF (CC3) THEN 992 IOFF = ICKBD(ISYCKD,ISYMB) 993 * + NCKATR(ISYCKD)*(B - 1) + 1 994 IF (NCKATR(ISYCKD) .GT. 0) THEN 995 CALL GETWA2(LUDKBC,FNDKBC,WORK(KTRVI8),IOFF, 996 * NCKATR(ISYCKD)) 997 ENDIF 998 ELSE 999C 1000 IOFF = ICKAD(ISYCKD,ISYMB) 1001 * + NCKA(ISYCKD)*(B - 1) + 1 1002 IF (NCKA(ISYCKD) .GT. 0) THEN 1003 CALL GETWA2(LUDKBC3,FNDKBC3,WORK(KINTVI),IOFF, 1004 * NCKA(ISYCKD)) 1005 ENDIF 1006C 1007 CALL CCSDT_TRVIR(WORK(KINTVI),WORK(KTRVI8), 1008 * WORK(KCMO),ISYMB,B,ISINT2, 1009 * WORK(KEND5),LWRK5) 1010 ENDIF 1011C 1012 CALL CCSDT_SRTVIR(WORK(KTRVI8),WORK(KTRVI9), 1013 * WORK(KEND4),LWRK4,ISYMB,ISINT2) 1014C 1015 IF (CC3) THEN 1016 IOFF = ICKBD(ISYCKD,ISYMB) 1017 * + NCKATR(ISYCKD)*(B - 1) + 1 1018 IF (NCKATR(ISYCKD) .GT. 0) THEN 1019 CALL GETWA2(LUDKBC3,FNDKBC3,WORK(KTRVI16),IOFF, 1020 * NCKATR(ISYCKD)) 1021 ENDIF 1022 CALL CCSDT_SRVIR3(WORK(KTRVI16),WORK(KEND5), 1023 * ISYMB,B,ISINT2) 1024 CALL CCSDT_SRTVIR(WORK(KTRVI16),WORK(KTRVI17), 1025 * WORK(KEND4),LWRK4,ISYMB,ISINT2) 1026 ENDIF 1027C 1028C---------------------------------------------------------- 1029C Read virtual integrals used in second U 1030C---------------------------------------------------------- 1031C 1032C 1033 IF (CC3) THEN 1034 IOFF = ICKAD(ISCKD2,ISYMB) 1035 * + NCKA(ISCKD2)*(B - 1) + 1 1036 IF (NCKA(ISYCKD) .GT. 0) THEN 1037 CALL GETWA2(LUDELD,FNDELD,WORK(KINTVI),IOFF, 1038 * NCKA(ISCKD2)) 1039 ENDIF 1040C 1041 CALL CCSDT_TRVIR(WORK(KINTVI),WORK(KTRVI10), 1042 * WORK(KLAMDH),ISYMB,B,ISINT2, 1043 * WORK(KEND5),LWRK5) 1044C 1045 ELSE 1046C 1047 IOFF = ICKAD(ISYCKD,ISYMB) 1048 * + NCKA(ISYCKD)*(B - 1) + 1 1049 IF (NCKA(ISYCKD) .GT. 0) THEN 1050 CALL GETWA2(LU3VI,FN3VI,WORK(KINTVI),IOFF, 1051 * NCKA(ISYCKD)) 1052 ENDIF 1053C 1054 CALL CCSDT_TRVIR(WORK(KINTVI),WORK(KTRVI10), 1055 * WORK(KCMO),ISYMB,B,ISYMOP, 1056 * WORK(KEND5),LWRK5) 1057 ENDIF 1058C 1059C------------------------------------------------------------------------ 1060C Read and transform integrals used in second S-bar and U-bar 1061C NOT used for CC3 1062C------------------------------------------------------------------------ 1063C 1064 IF (.NOT. RELORB) THEN 1065C 1066 IF (CC3) THEN 1067 IOFF = ICKBD(ISYCKD,ISYMB) 1068 * + NCKATR(ISYCKD)*(B-1) + 1 1069 IF (NCKATR(ISYCKD) .GT. 0) THEN 1070 CALL GETWA2(LU3FOP2X,FN3FOP2X,WORK(KTRVI11), 1071 * IOFF,NCKATR(ISYCKD)) 1072 ENDIF 1073 ELSE 1074 IOFF = ICKAD(ISYCKD,ISYMB) 1075 * + NCKA(ISYCKD)*(B-1) + 1 1076 IF (NCKA(ISYCKD) .GT. 0) THEN 1077 CALL GETWA2(LU3FOP2,FN3FOP2,WORK(KINTVI), 1078 * IOFF,NCKA(ISYCKD)) 1079 ENDIF 1080C 1081 CALL CCSDT_TRVIR(WORK(KINTVI),WORK(KTRVI11), 1082 * WORK(KCMO),ISYMB,B,ISYMOP, 1083 * WORK(KEND5),LWRK5) 1084C 1085 ENDIF 1086C 1087 CALL CCSDT_SRTVIR(WORK(KTRVI11),WORK(KTRVI12), 1088 * WORK(KEND5),LWRK5,ISYMB, 1089 * ISINT2) 1090C 1091 IF (CC3) THEN 1092 IOFF = ICKBD(ISYCKD,ISYMB) 1093 * + NCKATR(ISYCKD)*(B - 1) + 1 1094 IF (NCKATR(ISYCKD) .GT. 0) THEN 1095 CALL GETWA2(LU3FOPX,FN3FOPX,WORK(KTRVI13), 1096 * IOFF,NCKATR(ISYCKD)) 1097 ENDIF 1098C 1099 IOFF = ICKAD(ISYCKD,ISYMB) 1100 * + NCKA(ISYCKD)*(B - 1) + 1 1101 IF (NCKA(ISYCKD) .GT. 0) THEN 1102 CALL GETWA2(LU3VI,FN3VI,WORK(KINTVI),IOFF, 1103 * NCKA(ISYCKD)) 1104 ENDIF 1105C 1106 CALL CCSDT_TRVIR(WORK(KINTVI),WORK(KTRVI20), 1107 * WORK(KLAMDP),ISYMB,B,ISYMOP, 1108 * WORK(KEND4),LWRK4) 1109 ELSE 1110 IOFF = ICKAD(ISYCKD,ISYMB) 1111 * + NCKA(ISYCKD)*(B-1) + 1 1112 IF (NCKA(ISYCKD) .GT. 0) THEN 1113 CALL GETWA2(LU3FOP,FN3FOP,WORK(KINTVI),IOFF, 1114 * NCKA(ISYCKD)) 1115 ENDIF 1116C 1117 CALL CCSDT_TRVIR(WORK(KINTVI),WORK(KTRVI13), 1118 * WORK(KCMO),ISYMB,B,ISINT2, 1119 * WORK(KEND5),LWRK5) 1120 ENDIF 1121 ENDIF 1122C 1123C------------------------------------------------------------------- 1124C Calculate the S(ci,bk,dj) matrix for T3 for B,D. 1125C------------------------------------------------------------------- 1126C 1127 DTIME = SECOND() 1128 CALL CC3_SMAT(0.0D0,T2TP,ISYMT2,WORK(KTMAT), 1129 * WORK(KTRVI0), 1130 * WORK(KTRVI2),WORK(KTROC0),ISINT2, 1131 * WORK(KFOCKD),WORK(KDIAG), 1132 * WORK(KSMAT),WORK(KEND4),LWRK4, 1133 * WORK(KINDEX),WORK(KINDSQ),LENSQ, 1134 * ISYMB,B,ISYMD,D) 1135C 1136 CALL T3_FORBIDDEN(WORK(KSMAT),ISYMIM,ISYMB,B,ISYMD,D) 1137C 1138 DTIME = SECOND() - DTIME 1139 TISMAT = TISMAT + DTIME 1140C 1141 IF (IPRINT .GT. 55) THEN 1142 XSMAT = DDOT(NCKIJ(ISCKIJ),WORK(KSMAT),1, 1143 * WORK(KSMAT),1) 1144 WRITE(LUPRI,*) 'Norm of SMAT ',XSMAT 1145 ENDIF 1146C 1147C------------------------------------------------------------------- 1148C Calculate the S(ci,bk,dj) matrix for T3 for D,B. 1149C------------------------------------------------------------------- 1150C 1151 DTIME = SECOND() 1152 CALL CC3_SMAT(0.0D0,T2TP,ISYMT2,WORK(KTMAT), 1153 * WORK(KTRVI8), 1154 * WORK(KTRVI9),WORK(KTROC0),ISINT2, 1155 * WORK(KFOCKD),WORK(KDIAG), 1156 * WORK(KSMAT3),WORK(KEND4),LWRK4, 1157 * WORK(KINDEX2),WORK(KINDSQ),LENSQ, 1158 * ISYMD,D,ISYMB,B) 1159C 1160 CALL T3_FORBIDDEN(WORK(KSMAT3),ISYMIM,ISYMD,D,ISYMB,B) 1161C 1162 DTIME = SECOND() - DTIME 1163 TISMAT = TISMAT + DTIME 1164C 1165 IF (IPRINT .GT. 55) THEN 1166 XSMAT = DDOT(NCKIJ(ISCKIJ),WORK(KSMAT3),1, 1167 * WORK(KSMAT3),1) 1168 WRITE(LUPRI,*) 'Norm of SMAT3 ',XSMAT 1169 ENDIF 1170C 1171C--------------------------------------------------------------------------- 1172C Calculate the S(ci,bk,dj) matrix for for B,D for T3-BAR. 1173C--------------------------------------------------------------------------- 1174C 1175 DTIME = SECOND() 1176C 1177 CALL DZERO(WORK(KSMAT2),NCKIJ(ISCKIJ)) 1178C 1179 IF (CC3) THEN 1180 CALL CCFOP_SMAT(0.0D0,L1AM,ISYML1,L2TP,ISYML2, 1181 * WORK(KTMAT), 1182 * WORK(KFCKBA),WORK(KXIAJB),ISINT1, 1183 * WORK(KTRVI14),WORK(KTRVI15), 1184 * WORK(KTRVI4),WORK(KTRVI5), 1185 * WORK(KTROC01),WORK(KTROC21), 1186 * ISINT2,WORK(KFOCKD),WORK(KDIAG), 1187 * WORK(KSMAT2),WORK(KEND4),LWRK4, 1188 * WORK(KINDEX),WORK(KINDSQ),LENSQ, 1189 * ISYMB,B,ISYMD,D) 1190C 1191 CALL DSCAL(NCKIJ(ISCKIJ),HALF,WORK(KSMAT2),1) 1192C 1193 ELSE 1194C 1195 CALL CCFOP_SMAT(0.0D0,T1AM,ISYMT1,WORK(KT2TCME), 1196 * ISYMT2,WORK(KTMAT), 1197 * WORK(KFCKBA),WORK(KXIAJB),ISINT1, 1198 * WORK(KTRVI0),WORK(KTRVI2), 1199 * WORK(KTRVI4),WORK(KTRVI5), 1200 * WORK(KTROC0),WORK(KTROC2), 1201 * ISINT2,WORK(KFOCKD),WORK(KDIAG), 1202 * WORK(KSMAT2),WORK(KEND4),LWRK4, 1203 * WORK(KINDEX),WORK(KINDSQ),LENSQ, 1204 * ISYMB,B,ISYMD,D) 1205 ENDIF 1206C 1207 CALL T3_FORBIDDEN(WORK(KSMAT2),ISYMIM,ISYMB,B,ISYMD,D) 1208C 1209 DTIME = SECOND() - DTIME 1210 TISMAT = TISMAT + DTIME 1211C 1212 IF (IPRINT .GT. 55) THEN 1213 XSMAT = DDOT(NCKIJ(ISCKIJ),WORK(KSMAT2),1, 1214 * WORK(KSMAT2),1) 1215 WRITE(LUPRI,*) 'Norm of SMAT-BAR-1 ',XSMAT 1216 ENDIF 1217C 1218C----------------------------------------------------------------------------- 1219C Calculate the S(ci,bk,dj) matrix for for D,B for T3-BAR. 1220C----------------------------------------------------------------------------- 1221C 1222 IF (.NOT. RELORB) THEN 1223C 1224 DTIME = SECOND() 1225C 1226 CALL DZERO(WORK(KSMAT4),NCKIJ(ISCKIJ)) 1227C 1228 IF (CC3) THEN 1229 CALL CCFOP_SMAT(0.0D0,L1AM,ISYML1,L2TP, 1230 * ISYML2,WORK(KTMAT),WORK(KFCKBA), 1231 * WORK(KXIAJB),ISINT1, 1232 * WORK(KTRVI16),WORK(KTRVI17), 1233 * WORK(KTRVI11),WORK(KTRVI12), 1234 * WORK(KTROC01),WORK(KTROC21), 1235 * ISINT2,WORK(KFOCKD),WORK(KDIAG), 1236 * WORK(KSMAT4),WORK(KEND4),LWRK4, 1237 * WORK(KINDEX2),WORK(KINDSQ), 1238 * LENSQ,ISYMD,D,ISYMB,B) 1239C 1240 CALL DSCAL(NCKIJ(ISCKIJ),HALF,WORK(KSMAT4),1) 1241C 1242 ELSE 1243C 1244 CALL CCFOP_SMAT(0.0D0,T1AM,ISYMT1,WORK(KT2TCME), 1245 * ISYMT2,WORK(KTMAT),WORK(KFCKBA), 1246 * WORK(KXIAJB),ISINT1, 1247 * WORK(KTRVI8),WORK(KTRVI9), 1248 * WORK(KTRVI11),WORK(KTRVI12), 1249 * WORK(KTROC0),WORK(KTROC2), 1250 * ISINT2,WORK(KFOCKD),WORK(KDIAG), 1251 * WORK(KSMAT4),WORK(KEND4),LWRK4, 1252 * WORK(KINDEX2),WORK(KINDSQ), 1253 * LENSQ,ISYMD,D,ISYMB,B) 1254 ENDIF 1255C 1256 CALL T3_FORBIDDEN(WORK(KSMAT4),ISYMIM, 1257 * ISYMD,D,ISYMB,B) 1258C 1259 DTIME = SECOND() - DTIME 1260 TISMAT = TISMAT + DTIME 1261C 1262 IF (IPRINT .GT. 55) THEN 1263 XSMAT = DDOT(NCKIJ(ISCKIJ),WORK(KSMAT4),1, 1264 * WORK(KSMAT4),1) 1265 WRITE(LUPRI,*) 'Norm of SMAT-BAR-2 ',XSMAT 1266 ENDIF 1267C 1268 ENDIF 1269C 1270C-------------------------------------------------- 1271C Calculate Q(ci,jk) for fixed b,d. 1272C-------------------------------------------------- 1273C 1274 DTIME = SECOND() 1275 CALL CC3_QMAT(0.0D0,T2TP,ISYMT2,WORK(KTRVI3), 1276 * WORK(KTROC0),ISINT2,WORK(KFOCKD), 1277 * WORK(KDIAG),WORK(KQMAT), 1278 * WORK(KEND4),LWRK4,WORK(KINDSQ),LENSQ, 1279 * ISYMB,B,ISYMD,D) 1280C 1281 CALL T3_FORBIDDEN(WORK(KQMAT),ISYMIM,ISYMB,B,ISYMD,D) 1282C 1283 DTIME = SECOND() - DTIME 1284 TIQMAT = TIQMAT + DTIME 1285C 1286 IF (IPRINT .GT. 55) THEN 1287 XQMAT = DDOT(NCKIJ(ISCKIJ),WORK(KQMAT),1, 1288 * WORK(KQMAT),1) 1289 WRITE(LUPRI,*) 'Norm of QMAT ',XQMAT 1290 ENDIF 1291C 1292C------------------------------------------------------------------- 1293C Calculate Q(ci,jk) for fixed b,d for t3-bar. 1294C------------------------------------------------------------------- 1295C 1296 DTIME = SECOND() 1297C 1298 CALL DZERO(WORK(KQMAT2),NCKIJ(ISCKIJ)) 1299C 1300 IF (CC3) THEN 1301 CALL CCFOP_QMAT(0.0D0,L1AM,ISYML1,L2TP,ISYML2, 1302 * WORK(KTMAT),WORK(KFCKBA), 1303 * WORK(KXIAJB),ISINT1,WORK(KTRVI18), 1304 * WORK(KTRVI6),WORK(KTROC01), 1305 * WORK(KTROC21),ISINT2,WORK(KFOCKD), 1306 * WORK(KDIAG),WORK(KQMAT2), 1307 * WORK(KEND4),LWRK4,WORK(KINDSQ), 1308 * LENSQ,ISYMB,B,ISYMD,D) 1309C 1310 CALL DSCAL(NCKIJ(ISCKIJ),HALF,WORK(KQMAT2),1) 1311C 1312 ELSE 1313 CALL CCFOP_QMAT(0.0D0,T1AM,ISYMT1, 1314 * WORK(KT2TCME),ISYMT2, 1315 * WORK(KTMAT),WORK(KFCKBA), 1316 * WORK(KXIAJB),ISINT1,WORK(KTRVI3), 1317 * WORK(KTRVI6),WORK(KTROC0), 1318 * WORK(KTROC2),ISINT2,WORK(KFOCKD), 1319 * WORK(KDIAG),WORK(KQMAT2), 1320 * WORK(KEND4),LWRK4,WORK(KINDSQ), 1321 * LENSQ,ISYMB,B,ISYMD,D) 1322 ENDIF 1323C 1324 CALL T3_FORBIDDEN(WORK(KQMAT2),ISYMIM,ISYMB,B,ISYMD,D) 1325C 1326 DTIME = SECOND() - DTIME 1327 TIQMAT = TIQMAT + DTIME 1328C 1329 IF (IPRINT .GT. 55) THEN 1330 XQMAT = DDOT(NCKIJ(ISCKIJ),WORK(KQMAT2),1, 1331 * WORK(KQMAT2),1) 1332 WRITE(LUPRI,*) 'Norm of QMAT-BAR ',XQMAT 1333 ENDIF 1334C 1335C-------------------------------------------------- 1336C Calculate U(ci,jk) for fixed b,d. 1337C-------------------------------------------------- 1338C 1339 DTIME = SECOND() 1340 CALL CC3_UMAT(0.0D0,T2TP,ISYMT2,WORK(KTRVI1), 1341 * WORK(KTROC02),ISINT2,WORK(KFOCKD), 1342 * WORK(KDIAG),WORK(KUMAT),WORK(KTMAT), 1343 * WORK(KEND4),LWRK4,WORK(KINDSQ),LENSQ, 1344 * ISYMB,B,ISYMD,D) 1345C 1346 CALL T3_FORBIDDEN(WORK(KUMAT),ISYMIM,ISYMB,B,ISYMD,D) 1347C 1348 DTIME = SECOND() - DTIME 1349 TIQMAT = TIQMAT + DTIME 1350C 1351 IF (IPRINT .GT. 55) THEN 1352 XQMAT = DDOT(NCKIJ(ISCKIJ),WORK(KUMAT),1, 1353 * WORK(KUMAT),1) 1354 WRITE(LUPRI,*) 'Norm of UMAT ',XQMAT 1355 ENDIF 1356C 1357C-------------------------------------------------- 1358C Calculate U(ci,jk) for fixed d,b. 1359C-------------------------------------------------- 1360C 1361 DTIME = SECOND() 1362 CALL CC3_UMAT(0.0D0,T2TP,ISYMT2,WORK(KTRVI10), 1363 * WORK(KTROC02),ISINT2,WORK(KFOCKD), 1364 * WORK(KDIAG),WORK(KUMAT3),WORK(KTMAT), 1365 * WORK(KEND4),LWRK4,WORK(KINDSQ),LENSQ, 1366 * ISYMD,D,ISYMB,B) 1367C 1368 CALL T3_FORBIDDEN(WORK(KUMAT3),ISYMIM,ISYMD,D,ISYMB,B) 1369C 1370 DTIME = SECOND() - DTIME 1371 TIQMAT = TIQMAT + DTIME 1372C 1373 IF (IPRINT .GT. 55) THEN 1374 XQMAT = DDOT(NCKIJ(ISCKIJ),WORK(KUMAT),1, 1375 * WORK(KUMAT),1) 1376 WRITE(LUPRI,*) 'Norm of UMAT3 ',XQMAT 1377 ENDIF 1378C 1379C----------------------------------------------------------------- 1380C Calculate U(ci,jk) for fixed b,d for t3-bar. 1381C----------------------------------------------------------------- 1382C 1383 DTIME = SECOND() 1384C 1385 CALL DZERO(WORK(KUMAT2),NCKIJ(ISCKIJ)) 1386C 1387 IF (CC3) THEN 1388 CALL CCFOP_UMAT(0.0D0,L1AM,ISYML1,L2TP,ISYML2, 1389 * WORK(KXIAJB),ISINT1,WORK(KFCKBA), 1390 * WORK(KTRVI19),WORK(KTRVI7), 1391 * WORK(KTROC03),WORK(KTROC23),ISINT2, 1392 * WORK(KFOCKD),WORK(KDIAG), 1393 * WORK(KUMAT2), 1394 * WORK(KTMAT),WORK(KEND4),LWRK4, 1395 * WORK(KINDSQ),LENSQ,ISYMB,B,ISYMD,D) 1396C 1397 CALL DSCAL(NCKIJ(ISCKIJ),HALF,WORK(KUMAT2),1) 1398C 1399 ELSE 1400 CALL CCFOP_UMAT(0.0D0,T1AM,ISYMT1,WORK(KT2TCME), 1401 * ISYMT2, 1402 * WORK(KXIAJB),ISINT1,WORK(KFCKBA), 1403 * WORK(KTRVI1),WORK(KTRVI7), 1404 * WORK(KTROC02),WORK(KTROC22),ISINT2, 1405 * WORK(KFOCKD),WORK(KDIAG), 1406 * WORK(KUMAT2), 1407 * WORK(KTMAT),WORK(KEND4),LWRK4, 1408 * WORK(KINDSQ),LENSQ,ISYMB,B,ISYMD,D) 1409 ENDIF 1410C 1411 CALL T3_FORBIDDEN(WORK(KUMAT2),ISYMIM,ISYMB,B,ISYMD,D) 1412C 1413 DTIME = SECOND() - DTIME 1414 TIQMAT = TIQMAT + DTIME 1415C 1416 IF (IPRINT .GT. 55) THEN 1417 XQMAT = DDOT(NCKIJ(ISCKIJ),WORK(KUMAT2),1, 1418 * WORK(KUMAT2),1) 1419 WRITE(LUPRI,*) 'Norm of UMAT-BAR-1 ',XQMAT 1420 ENDIF 1421C 1422C----------------------------------------------------------------- 1423C Calculate U(ci,jk) for fixed d,b for t3-bar. 1424C----------------------------------------------------------------- 1425C 1426 IF (.NOT. RELORB) THEN 1427C 1428 DTIME = SECOND() 1429C 1430 CALL DZERO(WORK(KUMAT4),NCKIJ(ISCKIJ)) 1431C 1432 IF (CC3) THEN 1433 CALL CCFOP_UMAT(0.0D0,L1AM,ISYML1,L2TP,ISYML2, 1434 * WORK(KXIAJB),ISINT1, 1435 * WORK(KFCKBA),WORK(KTRVI20), 1436 * WORK(KTRVI13),WORK(KTROC03), 1437 * WORK(KTROC23),ISINT2, 1438 * WORK(KFOCKD),WORK(KDIAG), 1439 * WORK(KUMAT4),WORK(KTMAT), 1440 * WORK(KEND4),LWRK4,WORK(KINDSQ), 1441 * LENSQ,ISYMD,D,ISYMB,B) 1442C 1443 CALL DSCAL(NCKIJ(ISCKIJ),HALF,WORK(KUMAT4),1) 1444C 1445 ELSE 1446 CALL CCFOP_UMAT(0.0D0,T1AM,ISYMT1,WORK(KT2TCME), 1447 * ISYMT2,WORK(KXIAJB),ISINT1, 1448 * WORK(KFCKBA),WORK(KTRVI10), 1449 * WORK(KTRVI13),WORK(KTROC02), 1450 * WORK(KTROC22),ISINT2, 1451 * WORK(KFOCKD),WORK(KDIAG), 1452 * WORK(KUMAT4),WORK(KTMAT), 1453 * WORK(KEND4),LWRK4,WORK(KINDSQ), 1454 * LENSQ,ISYMD,D,ISYMB,B) 1455 ENDIF 1456C 1457 CALL T3_FORBIDDEN(WORK(KUMAT4),ISYMIM, 1458 * ISYMD,D,ISYMB,B) 1459C 1460 DTIME = SECOND() - DTIME 1461 TIQMAT = TIQMAT + DTIME 1462C 1463 IF (IPRINT .GT. 55) THEN 1464 XQMAT = DDOT(NCKIJ(ISCKIJ),WORK(KUMAT4),1, 1465 * WORK(KUMAT4),1) 1466 WRITE(LUPRI,*) 'Norm of UMAT-BAR-2 ',XQMAT 1467 ENDIF 1468C 1469 ENDIF 1470C 1471C----------------------------------------------------------- 1472C Construct Kappabar_{aa} and Kappabar_{ii} 1473C----------------------------------------------------------- 1474C 1475 IF ((.NOT. CC3) .AND. (RELORB)) THEN 1476C 1477 CALL CCSDT_KAPPADIAG(WORK(KKAPAA),WORK(KKAPII), 1478 * WORK(KSMAT2),WORK(KSMAT), 1479 * WORK(KSMAT3),WORK(KUMAT2), 1480 * WORK(KUMAT),WORK(KUMAT3), 1481 * WORK(KTMAT),WORK(KINDSQ), 1482 * LENSQ,ISCKIJ, 1483 * WORK(KEND4),LWRK4) 1484C 1485 ENDIF 1486C 1487C---------------------------------------------------------------- 1488C Calculate the three extra contributions to the 1489C one-electron density if nonrelaxed 1490C---------------------------------------------------------------- 1491C 1492 IF (.NOT. RELORB) THEN 1493 CALL CCFOP_NONREL(WORK(KOMG12),WORK(KDENSAB), 1494 * WORK(KDENSIJ),ISCKIJ, 1495 * WORK(KSMAT),WORK(KSMAT3), 1496 * WORK(KSMAT2),WORK(KSMAT4), 1497 * WORK(KUMAT),WORK(KUMAT3), 1498 * WORK(KUMAT2),WORK(KUMAT4), 1499 * WORK(KTMAT),T2TP,ISYMT2, 1500 * WORK(KINDSQ),LENSQ, 1501 * ISYMB,B,ISYMD,D, 1502 * WORK(KEND4),LWRK4) 1503 ENDIF 1504C 1505C--------------------------------------------- 1506C Contract with integrals. 1507C--------------------------------------------- 1508C 1509 DTIME = SECOND() 1510C 1511 IF ((.NOT. CC3) .AND. (RELORB)) THEN 1512C 1513 CALL CCFOP_DENVIR_SC(WORK(KVIR1),WORK(KVIR2), 1514 * WORK(KSMAT),WORK(KQMAT), 1515 * WORK(KTMAT),ISYMIM, 1516 * WORK(KT2TCME),ISYMT2, 1517 & WORK(KEND4),LWRK4, 1518 & WORK(KINDSQ),LENSQ, 1519 * ISYMB,B,ISYMD,D,1) 1520C 1521 IF ((IPRINT .GT. 55)) THEN 1522 XRMAT = DDOT(NCKATR(ISAIJ1),WORK(KVIR1),1, 1523 * WORK(KVIR1),1) 1524 WRITE(LUPRI,*) 'Norm DENS1 - CCFOP_DENVIR',XRMAT 1525 XRMAT = DDOT(NCKATR(ISAIJ1),WORK(KVIR2),1, 1526 * WORK(KVIR2),1) 1527 WRITE(LUPRI,*) 'Norm DENS2 - CCFOP_CONVIR',XRMAT 1528 ENDIF 1529C 1530 CALL CCFOP_DENVIR_SC(WORK(KVIR3),WORK(KVIR4), 1531 * WORK(KSMAT2),WORK(KQMAT2), 1532 * WORK(KTMAT),ISYMIM, 1533 * T2TP,ISYMT2,WORK(KEND4), 1534 * LWRK4,WORK(KINDSQ),LENSQ, 1535 * ISYMB,B,ISYMD,D,2) 1536C 1537 IF ((IPRINT .GT. 55)) THEN 1538 XRMAT = DDOT(NCKATR(ISAIJ1),WORK(KVIR1),1, 1539 * WORK(KVIR1),1) 1540 WRITE(LUPRI,*) 'Norm DENS1 - CCFOP_DENVIR',XRMAT 1541 XRMAT = DDOT(NCKATR(ISAIJ1),WORK(KVIR2),1, 1542 * WORK(KVIR2),1) 1543 WRITE(LUPRI,*) 'Norm DENS2 - CCFOP_CONVIR',XRMAT 1544 ENDIF 1545C 1546C 1547 CALL CCFOP_DENOCC_SC(WORK(KOCC1),WORK(KSMAT), 1548 * WORK(KQMAT),WORK(KTMAT),ISYMIM, 1549 * WORK(KT2TCME),ISYMT2,WORK(KEND4), 1550 * LWRK4,WORK(KINDSQ),LENSQ, 1551 * ISYMB,B,ISYMD,D,1) 1552 1553C 1554 CALL CCFOP_DENOCC_SC(WORK(KOCC2),WORK(KSMAT2), 1555 * WORK(KQMAT2),WORK(KTMAT),ISYMIM, 1556 * T2TP,ISYMT2,WORK(KEND4), 1557 * LWRK4,WORK(KINDSQ),LENSQ, 1558 * ISYMB,B,ISYMD,D,2) 1559C 1560 1561 IF ((IPRINT .GT. 55)) THEN 1562 XRMAT = DDOT(NCKIJ(ISYRES),WORK(KOCC1),1, 1563 * WORK(KOCC1),1) 1564 WRITE(LUPRI,*) 'Norm DENS1 - CCFOP_DENOCC',XRMAT 1565 XRMAT = DDOT(NCKIJ(ISYRES),WORK(KOCC2),1, 1566 * WORK(KOCC2),1) 1567 WRITE(LUPRI,*) 'Norm DENS2 - CCFOP_DENOCC',XRMAT 1568 END IF 1569C--------------------------------------- 1570C Calculate Omega22. 1571C--------------------------------------- 1572C 1573 DTIME = SECOND() 1574C 1575 CALL CCFOP_ONEL(WORK(KOMG22),WORK(KRMAT1), 1576 * WORK(KRMAT2),T1AM,WORK(KSMAT), 1577 * WORK(KTMAT),ISYMIM,ISINT1, 1578 * WORK(KINDSQ),LENSQ, 1579 * WORK(KEND4),LWRK4,ISYMB,B,ISYMD,D) 1580C 1581 IF ((IPRINT .GT. 55)) THEN 1582 RHO2N = DDOT(NT2AM(ISYRES),WORK(KOMG22),1, 1583 * WORK(KOMG22),1) 1584 WRITE(LUPRI,*) 'Norm of Rho22 after CC3_ONEL', 1585 * RHO2N 1586 ENDIF 1587C 1588 IF (IPRINT .GT. 220) THEN 1589 CALL AROUND('After CC3_ONEL: ') 1590 CALL CC_PRP(DUMMY,WORK(KOMG22),ISYRES,0,1) 1591 ENDIF 1592C 1593 IF (IPRINT .GT. 55) THEN 1594 XRMAT = DDOT(NCKI(ISAIJ1),WORK(KRMAT1),1, 1595 * WORK(KRMAT1),1) 1596 WRITE(LUPRI,*) 'Norm of RMAT1 -after ONEL',XRMAT 1597 XRMAT = DDOT(NCKI(ISAIJ2),WORK(KRMAT2),1, 1598 * WORK(KRMAT2),1) 1599 WRITE(LUPRI,*) 'Norm of RMAT2 -after ONEL',XRMAT 1600 XTMAT = DDOT(NCKIJ(ISCKIJ),WORK(KSMAT),1, 1601 * WORK(KSMAT),1) 1602 WRITE(LUPRI,*) 'Norm of SMAT -after ONEL',XSMAT 1603 XTMAT = DDOT(NCKIJ(ISCKIJ),WORK(KTMAT),1, 1604 * WORK(KTMAT),1) 1605 WRITE(LUPRI,*) 'Norm of TMAT -after ONEL',XTMAT 1606 ENDIF 1607C 1608 ENDIF ! RELORB 1609C 1610C--------------------------------------------------- 1611C Calculate Omega1. 1612C--------------------------------------------------- 1613C 1614 DTIME = SECOND() - DTIME 1615 TIOME1 = TIOME1 + DTIME 1616C 1617 IF (CC3) THEN 1618 CALL CCFOP_ONED(WORK(KOMG1),L2TP,ISYML2, 1619 * WORK(KSMAT),WORK(KTMAT),ISYMIM, 1620 * WORK(KINDSQ),LENSQ, 1621 * WORK(KEND4),LWRK4,ISYMB,B,ISYMD,D) 1622 ELSE 1623 CALL CCFOP_ONED(WORK(KOMG1),WORK(KT2TCME),ISYMT2, 1624 * WORK(KSMAT),WORK(KTMAT),ISYMIM, 1625 * WORK(KINDSQ),LENSQ, 1626 * WORK(KEND4),LWRK4,ISYMB,B,ISYMD,D) 1627 ENDIF 1628C 1629 IF ((IPRINT .GT. 55)) THEN 1630 XT2TP = DDOT(NT1AM(ISYMOP),WORK(KOMG1),1, 1631 * WORK(KOMG1),1) 1632 WRITE(LUPRI,*) 'Norm of 1 e- density : ',XT2TP 1633 ENDIF 1634C 1635 IF (IPRINT .GT. 220) THEN 1636 CALL AROUND('After CCFOP_ONED: ') 1637 CALL CC_PRP(WORK(KOMG1),DUMMY,ISYRES,1,0) 1638 ENDIF 1639C 1640 DTIME = SECOND() - DTIME 1641 TIOME1 = TIOME1 + DTIME 1642C 1643C--------------------------------------------------------- 1644C Accumulate the R2 matrix in Omega22 1645C--------------------------------------------------------- 1646C 1647 IF ((.NOT. CC3) .AND. (RELORB)) THEN 1648C 1649 CALL CC3_RACC(WORK(KOMG22),WORK(KRMAT2),ISYMB,B, 1650 * ISYRES) 1651C 1652 IF ((IPRINT .GT. 55)) THEN 1653 RHO2N = DDOT(NT2AM(ISYRES),WORK(KOMG22),1, 1654 * WORK(KOMG22),1) 1655 WRITE(LUPRI,*) 'Norm of Rho22-after CC3_RACC', 1656 * RHO2N 1657 ENDIF 1658C 1659 IF (IPRINT .GT. 220) THEN 1660 CALL AROUND('After CC3_RACC: ') 1661 CALL CC_PRP(DUMMY,WORK(KOMG22),ISYRES,0,1) 1662 ENDIF 1663C 1664 ENDIF 1665C 1666 ENDDO ! B 1667 ENDDO ! ISYMB 1668C 1669C--------------------------------------------------- 1670C Accumulate the R1 matrix in Omega22. 1671C--------------------------------------------------- 1672C 1673 IF ((.NOT. CC3) .AND. (RELORB)) THEN 1674C 1675 CALL CC3_RACC(WORK(KOMG22),WORK(KRMAT1),ISYMD,D,ISYRES) 1676C 1677 IF (IPRINT .GT. 55) THEN 1678 RHO2N = DDOT(NT2AM(ISYRES),WORK(KOMG22),1, 1679 * WORK(KOMG22),1) 1680 WRITE(LUPRI,*) 'Norm of Rho22-after CC3_RACC-2',RHO2N 1681 ENDIF 1682C 1683 IF (IPRINT .GT. 220) THEN 1684 CALL AROUND('After CC3_RACC-2: ') 1685 CALL CC_PRP(DUMMY,WORK(KOMG22),ISYRES,0,1) 1686 ENDIF 1687C 1688C-------------------------------------------------------------- 1689C Sort the two electron densities from T3 for a constant 1690C d and write them to file. 1691C-------------------------------------------------------------- 1692C 1693 IF (LWRK4 .LT. NCKATR(ISAIJ1)) THEN 1694 CALL QUIT('Exceeded memory in CCSDPT_DENS2 (sort)') 1695 ENDIF 1696! 1697!Sonia: following not working with symmetry 1698C 1699! CALL DEN_AIBSORT(WORK(KVIR1),WORK(KEND4),ISAIJ1) 1700C 1701! CALL DEN_AIBSORT(WORK(KVIR2),WORK(KEND4),ISAIJ1) 1702C 1703 CALL DEN_BIASORT(WORK(KVIR1),WORK(KEND4),ISAIJ1) 1704C 1705 CALL DEN_BIASORT(WORK(KVIR2),WORK(KEND4),ISAIJ1) 1706C 1707 IOFF = ICKBD(ISAIJ1,ISYMD) 1708 * + NCKATR(ISAIJ1)*(D-1) 1709 * + 1 1710 CALL PUTWA2(LUABI2,FNDABI2,WORK(KVIR1),IOFF, 1711 * NCKATR(ISAIJ1)) 1712C 1713 CALL PUTWA2(LUABI1,FNDABI1,WORK(KVIR2),IOFF, 1714 * NCKATR(ISAIJ1)) 1715C 1716 IF (IPRINT .GT. 55) THEN 1717 RHO1N = DDOT(NCKATR(ISAIJ1),WORK(KVIR1),1, 1718 * WORK(KVIR1),1) 1719 WRITE(LUPRI,*) 'Norm of VIR1 : ',RHO1N 1720 RHO1N = DDOT(NCKATR(ISAIJ1),WORK(KVIR2),1, 1721 * WORK(KVIR2),1) 1722 WRITE(LUPRI,*) 'Norm of VIR2 : ',RHO1N 1723 ENDIF 1724C 1725C---------------------------------------------------------------------- 1726C Sort the two electron densities from T3-bar for a constant 1727C d and write them to file. 1728C---------------------------------------------------------------------- 1729C 1730 IF (LWRK4 .LT. NCKATR(ISAIJ1)) THEN 1731 CALL QUIT('Exceeded memory in CCSDPT_DENS2 (sort)') 1732 ENDIF 1733C 1734! 1735!Sonia: following not working with symmetry 1736C 1737! CALL DEN_AIBSORT(WORK(KVIR3),WORK(KEND4),ISAIJ1) 1738C 1739! CALL DEN_AIBSORT(WORK(KVIR4),WORK(KEND4),ISAIJ1) 1740C 1741 CALL DEN_BIASORT(WORK(KVIR3),WORK(KEND4),ISAIJ1) 1742C 1743 CALL DEN_BIASORT(WORK(KVIR4),WORK(KEND4),ISAIJ1) 1744C 1745 IOFF = ICKBD(ISAIJ1,ISYMD) 1746 * + NCKATR(ISAIJ1)*(D-1) 1747 * + 1 1748 CALL PUTWA2(LUABI4,FNDABI4,WORK(KVIR3),IOFF, 1749 * NCKATR(ISAIJ1)) 1750C 1751 CALL PUTWA2(LUABI3,FNDABI3,WORK(KVIR4),IOFF, 1752 * NCKATR(ISAIJ1)) 1753C 1754 IF (IPRINT .GT. 55) THEN 1755 RHO1N = DDOT(NCKATR(ISAIJ1),WORK(KVIR3),1, 1756 * WORK(KVIR3),1) 1757 WRITE(LUPRI,*) 'Norm of VIR3 : ',RHO1N 1758 RHO1N = DDOT(NCKATR(ISAIJ1),WORK(KVIR4),1, 1759 * WORK(KVIR4),1) 1760 WRITE(LUPRI,*) 'Norm of VIR4 : ',RHO1N 1761 ENDIF 1762C 1763 ENDIF ! RELORB 1764C 1765 ENDDO ! D 1766 ENDDO ! ISYMD 1767C 1768C ------------------------------------------------------------- 1769C The following section is taken from old code, as the new does 1770C not seem to work with symmetry 1771C Sonia, Dec 2002 1772C ------------------------------------------------------------- 1773C 1774 IF ((.NOT. CC3)) THEN 1775 1776 KEND4 = KENDSV 1777 LWRK4 = LWORK - KEND4 1778 1779 IF (RELORB) THEN 1780C 1781 NTOT = 0 1782 DO ISYMD = 1, NSYM 1783 IF (NT1AM(ISYMD) .GT. NTOT) THEN 1784 NTOT = NT1AM(ISYMD) 1785 ENDIF 1786 ENDDO 1787C 1788 KVIR1 = KEND4 1789 KVIR2 = KVIR1 + NTOT 1790 KEND4 = KVIR2 + NTOT 1791 LWRK4 = LWORK - KEND4 1792 1793C-------------------------------------------- 1794C Resort the densities occ1 and occ2 1795C-------------------------------------------- 1796 IF ((IPRINT.GT.55).or.(locdbg)) THEN 1797 RHO2N = ddot(nckij(isyres),work(kocc1),1,work(kocc1),1) 1798 WRITE(lupri,*) 'Norm of Occ1 (before cc3_lsort1) : ',rho2n 1799 RHO2N = ddot(nckij(isyres),work(kocc2),1,work(kocc2),1) 1800 WRITE(lupri,*) 'Norm of Occ2 (before cc3_lsort1) : ',rho2n 1801 END IF 1802C 1803 CALL CC3_LSORT1(WORK(KOCC1),ISYRES,WORK(KEND4),LWRK4,3) 1804C 1805 CALL CC3_LSORT1(WORK(KOCC2),ISYRES,WORK(KEND4),LWRK4,3) 1806C 1807 IF ((IPRINT.GT.55).or.(locdbg)) THEN 1808 rho2n = ddot(nckij(isyres),work(kocc1),1,work(kocc1),1) 1809 write(lupri,*) 'Norm of Occ1 (after cc3_lsort1) : ',rho2n 1810 rho2n = ddot(nckij(isyres),work(kocc2),1,work(kocc2),1) 1811 write(lupri,*) 'Norm of Occ2 (after cc3_lsort1) : ',rho2n 1812 END IF 1813 1814 ENDIF 1815 END IF 1816C 1817C-------------------------------------- End of old part.... 1818C--------------------------------------------------------- 1819C Construct 2*C-E of work(komg22) and write to file. 1820C--------------------------------------------------------- 1821C 1822 IF ((.NOT. CC3) .AND. (RELORB)) THEN 1823 IOPTTCME = 1 1824 ISYOPE = ISYRES 1825 CALL CCSD_TCMEPK(WORK(KOMG22),1.0D0,ISYOPE,IOPTTCME) 1826C 1827 IF ((IPRINT .GT. 55)) THEN 1828 RHO2N = DDOT(NT2AM(ISYRES),WORK(KOMG22),1,WORK(KOMG22),1) 1829 WRITE(LUPRI,*) 'Norm of Rho22 at the end ',RHO2N 1830 ENDIF 1831C 1832 IF (IPRINT .GT. 100) THEN 1833 CALL AROUND('TWO ELECTRON DENSITY : D_{IAJB}') 1834 CALL CC_PRP(T1AM,WORK(KOMG22),ISYRES,0,1) 1835 ENDIF 1836C 1837 IF (NT2AM(ISYRES) .GT. 0) THEN 1838 IOFF = 1 1839 CALL PUTWA2(LUPTIAJB,FNDIAJB,WORK(KOMG22),IOFF, 1840 * NT2AM(ISYRES)) 1841 ENDIF 1842C 1843 IF (LDEBUG .AND. (.NOT. CC3)) THEN 1844 LENGTH = IRAT*NT2AM(ISYRES) 1845C 1846 REWIND(LUIAJB) 1847 CALL READI(LUIAJB,LENGTH,WORK(KXIAJB)) 1848 CALL CCLR_DIASCL(WORK(KXIAJB),0.5D0,ISYMop) 1849 CALL DSCAL(NT2AM(ISYMOP),2.0D0,WORK(KOMG22),1) 1850C 1851 XQMAT = DDOT(NT2AM(ISYRES),WORK(KOMG22),1,WORK(KXIAJB),1) 1852 WRITE(LUPRI,*) 'DEBUGGING CCSD(T) : E5 = ',XQMAT 1853 ENDIF 1854C 1855 ENDIF 1856C 1857C--------------------------------------- 1858C Scale and store the 1 e- density : 1859C--------------------------------------- 1860C 1861 IF (CC3) THEN 1862 CALL DSCAL(NT1AM(ISYRES),-ONE,WORK(KOMG1),1) 1863 ELSE 1864 CALL DSCAL(NT1AM(ISYRES),-TWO,WORK(KOMG1),1) 1865 ENDIF 1866C 1867 IF (IPRINT .GT. 55) THEN 1868 RHO1N = DDOT(NT1AM(ISYRES),WORK(KOMG1),1,WORK(KOMG1),1) 1869 WRITE(LUPRI,*) 'Norm of OMEG1 at the end : ',RHO1N 1870 ENDIF 1871C 1872 IF (IPRINT .GT. 100) THEN 1873 CALL AROUND('1 e- in CCSDPT_DENS2 : ') 1874 CALL CC_PRP(WORK(KOMG1),DUMMY,ISYRES,1,0) 1875 ENDIF 1876C 1877 IF (NT1AM(ISYRES) .GT. 0) THEN 1878 IOFF = 1 1879 CALL PUTWA2(LUPTIA,FNDPTIA,WORK(KOMG1),IOFF,NT1AM(ISYRES)) 1880 ENDIF 1881C 1882C---------------------------------------------------------- 1883C Add the 1e- to the d_{iajk} for j=k and for i=k 1884C---------------------------------------------------------- 1885C 1886 IF ((IPRINT .GT. 55).or.locdbg) THEN 1887 RHO1N = DDOT(NCKIJ(ISYRES),WORK(KOCC1),1,WORK(KOCC1),1) 1888 WRITE(LUPRI,*) 'Norm of OCC1 (iajk) (before dens1to2) =',RHO1N 1889 ENDIF 1890C 1891 IF ((.NOT. CC3) .AND. (RELORB)) THEN 1892 CALL DENS1TO2(WORK(KOMG1),WORK(KOCC1),ISYRES) 1893 ENDIF 1894C 1895 IF ((IPRINT .GT. 55).or.locdbg) THEN 1896 RHO1N = DDOT(NCKIJ(ISYRES),WORK(KOCC1),1,WORK(KOCC1),1) 1897 WRITE(LUPRI,*) 'Norm of OCC1 (iajk) (after dens1to2) =',RHO1N 1898 ENDIF 1899C 1900C----------------------------------------------------------- 1901C If nonrel store the three extra terms on disc 1902C----------------------------------------------------------- 1903C 1904 IF (.NOT. RELORB) THEN 1905C 1906 IF (NMATAB(ISYRES) .GT. 0) THEN 1907 IOFF = 1 1908 CALL PUTWA2(LUPTAB,FNDPTAB,WORK(KDENSAB),IOFF,NMATAB(ISYRES)) 1909 ENDIF 1910C 1911 IF (NMATIJ(ISYRES) .GT. 0) THEN 1912 IOFF = 1 1913 CALL PUTWA2(LUPTIJ,FNDPTIJ,WORK(KDENSIJ),IOFF,NMATIJ(ISYRES)) 1914 ENDIF 1915C 1916 CALL DSCAL(NT1AM(ISYRES),-TWO,WORK(KOMG12),1) 1917 IF (NT1AM(ISYRES) .GT. 0) THEN 1918 IOFF = 1 1919 CALL PUTWA2(LUPTIA2,FNDPTIA2,WORK(KOMG12),IOFF,NT1AM(ISYRES)) 1920 ENDIF 1921C 1922 ENDIF 1923C 1924C--------------------------------------------------------------- 1925C Construct the total d(ab,ic) density stored as (ai,b,c) 1926C from the T3 amplitudes. 1927C--------------------------------------------------------------- 1928C 1929 IF ((.NOT. CC3) .AND. (RELORB)) THEN 1930C 1931 CALL DENSTORE(WORK(KVIR2),LUABI1,FNDABI1, 1932 * WORK(KVIR1),LUABI2,FNDABI2,ISYRES) 1933C 1934C 1935C--------------------------------------------------------------- 1936C Construct the total d(ab,ci) density stored as (bi,a,c) 1937C from the T3-bar amplitudes. 1938C--------------------------------------------------------------- 1939C 1940 CALL DENSTORE(WORK(KVIR4),LUABI3,FNDABI3, 1941 * WORK(KVIR3),LUABI4,FNDABI4,ISYRES) 1942C 1943C----------------------------------------- 1944C Store the d_{iajk} as kjia 1945C----------------------------------------- 1946C 1947 IF (NCKIJ(ISYRES) .GT. 0) THEN 1948 IOFF = 1 1949 CALL PUTWA2(LUIAJK,FNDIAJK,WORK(KOCC1),IOFF,NCKIJ(ISYRES)) 1950 ENDIF 1951C 1952C----------------------------------------- 1953C Store the d_{aijk} as jkia 1954C----------------------------------------- 1955C 1956 IF ((IPRINT .GT. 55).or.locdbg) THEN 1957 RHO1N = DDOT(NCKIJ(ISYRES),WORK(KOCC2),1,WORK(KOCC2),1) 1958 WRITE(LUPRI,*) 'Norm of OCC2 (aijk) = ',RHO1N 1959 ENDIF 1960C 1961 IF (NCKIJ(ISYRES) .GT. 0) THEN 1962 IOFF = 1 1963 CALL PUTWA2(LUAIJK,FNDAIJK,WORK(KOCC2),IOFF,NCKIJ(ISYRES)) 1964 ENDIF 1965C 1966C------------------------------------------ 1967C Store kappabar_{aa} and kappabar_{ii} 1968C------------------------------------------ 1969C 1970 IF ((IPRINT .GT. 55).or.locdbg) THEN 1971 RHO1N = DDOT(NRHFT,WORK(KKAPII),1,WORK(KKAPII),1) 1972 WRITE(LUPRI,*) 'Norm of KAPII : ',RHO1N 1973 RHO1N = DDOT(NVIRT,WORK(KKAPAA),1,WORK(KKAPAA),1) 1974 WRITE(LUPRI,*) 'Norm of KAPAA : ',RHO1N 1975 ENDIF 1976C 1977 LUKAPAB = -1 1978 LUKAPIJ = -1 1979 FNKAPAB = 'KAPAB' 1980 FNKAPIJ = 'KAPIJ' 1981 CALL WOPEN2(LUKAPAB,FNKAPAB,64,0) 1982 CALL WOPEN2(LUKAPIJ,FNKAPIJ,64,0) 1983C 1984 IF (NVIRT .GT. 0) THEN 1985 IOFF = 1 1986 CALL PUTWA2(LUKAPAB,FNKAPAB,WORK(KKAPAA),IOFF,NVIRT) 1987 ENDIF 1988C 1989 IF (NRHFT .GT. 0) THEN 1990 IOFF = 1 1991 CALL PUTWA2(LUKAPIJ,FNKAPIJ,WORK(KKAPII),IOFF,NRHFT) 1992 ENDIF 1993C 1994 CALL WCLOSE2(LUKAPAB,FNKAPAB,'KEEP') 1995 CALL WCLOSE2(LUKAPIJ,FNKAPIJ,'KEEP') 1996C 1997 ENDIF ! RELORB 1998C 1999C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 2000C SONIA: symmetrize/reorder some two-electron densities 2001C before closing, and generate backtransformed ones!!! 2002C Symmetrize the vir.vir.occ.vir and vir.vir.vir.occ 2003C Symmetrize the occ.occ.occ.vir and occ.occ.vir.occ 2004C and backtransform last index to delta 2005C---------------------------------------------------------- 2006C 2007 IF ((.NOT. CC3) .AND. RELORB) THEN 2008 if (.false.) then 2009 CALL SYMMBACK(MODEL,LUIAJK,FNDIAJK,LUAIJK,FNDAIJK, 2010 * LUABI1,FNDABI1,LUABI3,FNDABI3, 2011 * LUPTIAJB,FNDIAJB, 2012 * ISYRES,WORK(KEND4),LWRK4) 2013 else 2014 CALL SYMMBACK_1(MODEL,LUIAJK,FNDIAJK,LUAIJK,FNDAIJK, 2015 * LUABI1,FNDABI1,LUABI3,FNDABI3, 2016 * LUPTIAJB,FNDIAJB, 2017 * ISYRES,WORK(KEND4),LWRK4) 2018 end if 2019 ENDIF 2020C 2021C--------------------------------------- 2022C Close files. 2023C--------------------------------------- 2024C 2025 CALL WCLOSE2(LUPTIA,FNDPTIA,'KEEP') 2026C 2027 IF ((.NOT. CC3) .AND. (RELORB)) THEN 2028 CALL WCLOSE2(LUPTIAJB,FNDIAJB,'KEEP') 2029 CALL WCLOSE2(LUABI1,FNDABI1,'KEEP') 2030 CALL WCLOSE2(LUABI2,FNDABI2,'DELETE') 2031 CALL WCLOSE2(LUABI3,FNDABI3,'KEEP') 2032 CALL WCLOSE2(LUABI4,FNDABI4,'DELETE') 2033 CALL WCLOSE2(LUAIJK,FNDAIJK,'KEEP') 2034 CALL WCLOSE2(LUIAJK,FNDIAJK,'KEEP') 2035 ELSE 2036 CALL WCLOSE2(LUPTIA2,FNDPTIA2,'KEEP') 2037 CALL WCLOSE2(LUPTAB,FNDPTAB,'KEEP') 2038 CALL WCLOSE2(LUPTIJ,FNDPTIJ,'KEEP') 2039 ENDIF 2040C 2041C------------------- 2042C Print timings. 2043C------------------- 2044C 2045 IF (IPRINT .GT. 9) THEN 2046 WRITE(LUPRI,*) 2047 WRITE(LUPRI,*) 2048 WRITE(LUPRI,1) 'CC3_TRAN : ',TITRAN 2049 WRITE(LUPRI,1) 'CC3_SORT : ',TISORT 2050 WRITE(LUPRI,1) 'CC3_SMAT : ',TISMAT 2051 WRITE(LUPRI,1) 'CC3_QMAT : ',TIQMAT 2052 WRITE(LUPRI,1) 'CC3_OME1 : ',TIOME1 2053 WRITE(LUPRI,*) 2054 END IF 2055C 2056C------------- 2057C End 2058C------------- 2059C 2060 CALL QEXIT('CCSDPT_DENS2_SC') 2061C 2062 RETURN 2063C 2064 1 FORMAT(7X,'Time used in',2X,A12,F12.2,' seconds') 2065C 2066 END 2067C /* Deck ccfop_denvir_sc */ 2068 SUBROUTINE CCFOP_DENVIR_SC(RINTE1,RINTE2,SMAT,QMAT,TMAT,ISYMIM, 2069 * T2TCME,ISYMT2,WORK,LWORK,INDSQ,LENSQ, 2070 * ISYMB,B,ISYMD,D,IOPT) 2071C 2072C Kasper Hald, Fall 2001. 2073C OLd version working for CCSD(T) gradient, Sonia Coriani 2074C 2075C Calculate the two electron density (abic) for a constant index D, 2076C and add to the density RINTE1 and RINTE2. 2077C 2078C ISYMIM is the symmetry of the SMAT and TMAT intermdiates. 2079C ISYMT2 is the symmetry of the T2 amplitudes. 2080C 2081C IOPT = 1. Calculate the terms from T3AM. 2082C IOPT = 2. Calculate the terms from T3BAR. 2083C 2084 IMPLICIT NONE 2085C 2086#include "priunit.h" 2087#include "ccorb.h" 2088#include "ccsdinp.h" 2089#include "ccsdsym.h" 2090C 2091 INTEGER ISYMIM, ISYMT2, LWORK, LENSQ, ISYMB, ISYMD, IOPT 2092 INTEGER INDSQ(LENSQ,6) 2093 INTEGER INDEX, ISYRES, ISYMBD, ISCKIJ, LENGTH, ISYAIJ, ISYMAI 2094 INTEGER ISYMA, ISYMIJ, ISYMI, ISYMJ, NAI, KOFF1, KOFF2, KOFF3 2095 INTEGER ISYMK, ISYCIJ, ISYMC, ISYMAB, ISYMCK, NTOTCK, NTOTIJ 2096 INTEGER ISYBIJ, ISYMBJ, NTOTA, KEND1, LWRK1 2097C 2098#if defined (SYS_CRAY) 2099 REAL RINTE1(*), RINTE2(*), SMAT(*), QMAT(*) 2100 REAL TMAT(*), T2TCME(*), WORK(LWORK) 2101 REAL ZERO, ONE, TWO, HALF 2102#else 2103 DOUBLE PRECISION RINTE1(*), RINTE2(*), SMAT(*), QMAT(*) 2104 DOUBLE PRECISION TMAT(*), T2TCME(*), WORK(LWORK) 2105 DOUBLE PRECISION ZERO, ONE, TWO, HALF 2106#endif 2107 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0, HALF = 0.5D0) 2108C 2109C 2110C INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 2111C 2112 CALL QENTER('CCFOP_DENVIR_SC') 2113C 2114C----------------------------------------------- 2115C Sanity check and symmetry calculation. 2116C----------------------------------------------- 2117C 2118 IF (IOPT .NE. 1 .AND. IOPT .NE. 2) THEN 2119 CALL QUIT('Wrong IOPT in CCFOP_DENVIR_SC') 2120 ENDIF 2121C 2122 ISYRES = MULD2H(ISYMIM,ISYMT2) 2123C 2124 ISYMBD = MULD2H(ISYMB,ISYMD) 2125 ISCKIJ = MULD2H(ISYMBD,ISYMIM) 2126C 2127 LENGTH = NCKIJ(ISCKIJ) 2128C 2129C 2130C---------------------------------------- 2131C Sort the T2 for a constant B 2132C---------------------------------------- 2133C 2134 ISYAIJ = MULD2H(ISYMT2,ISYMB) 2135C 2136 IF (LWORK .LT. NCKATR(ISYAIJ)) THEN 2137 CALL QUIT('Exceeded work memory in CCFOP_DENVIR') 2138 ENDIF 2139C 2140 DO ISYMA = 1, NSYM 2141 ISYMIJ = MULD2H(ISYAIJ,ISYMA) 2142 ISYBIJ = MULD2H(ISYMIJ,ISYMB) 2143 DO ISYMI = 1, NSYM 2144C 2145 ISYMJ = MULD2H(ISYMIJ,ISYMI) 2146 ISYMBJ = MULD2H(ISYMJ,ISYMB) 2147 ISYMAI = MULD2H(ISYMA,ISYMI) 2148C 2149 DO A = 1, NVIR(ISYMA) 2150 DO I = 1, NRHF(ISYMI) 2151C 2152!Sonia: NAI is NOT used anywhere.... 2153! 2154! NAI = IT1AM(ISYMA,ISYMI) 2155! * + NVIR(ISYMA)*(I-1) + A 2156C 2157 KOFF1 = IT2SP(ISYBIJ,ISYMA) 2158 * + NCKI(ISYBIJ)*(A - 1) 2159 * + ICKI(ISYMBJ,ISYMI) 2160 * + NT1AM(ISYMBJ)*(I-1) 2161 * + IT1AM(ISYMB,ISYMJ) 2162 * + B 2163C 2164! KOFF2 = ISAIK(ISYMAI,ISYMJ) 2165! * + IT1AM(ISYMA,ISYMI) 2166! * + NVIR(ISYMA)*(I-1) 2167! * + A 2168C 2169! CALL DCOPY(NRHF(ISYMJ),T2TCME(KOFF1),NVIR(ISYMB), 2170! * WORK(KOFF2),NT1AM(ISYMAI)) 2171!Sonia: replace what above with the following: 2172! 2173 KOFF2 = IMAIJA(ISYMIJ,ISYMA) 2174 * + NMATIJ(ISYMIJ)*(A-1) 2175 * + IMATIJ(ISYMI,ISYMJ) 2176C * + NRHF(ISYMI)*(J-1) 2177 * + I 2178C 2179 CALL DCOPY(NRHF(ISYMJ),T2TCME(KOFF1),NVIR(ISYMB), 2180 * WORK(KOFF2),NRHF(ISYMI)) 2181C 2182 ENDDO ! I 2183 ENDDO ! A 2184 ENDDO ! ISYMI 2185 ENDDO ! ISYMA 2186C 2187C------------------------ 2188C First term. 2189C------------------------ 2190C 2191 DO I = 1,LENGTH 2192C 2193 IF (IOPT .EQ. 1) THEN 2194C 2195 TMAT(I) = TWO*SMAT(I) 2196 * - SMAT(INDSQ(I,1)) 2197 * - SMAT(INDSQ(I,5)) 2198 * + TWO*QMAT(INDSQ(I,3)) 2199 * - QMAT(INDSQ(I,2)) 2200 * - QMAT(INDSQ(I,4)) 2201C 2202 ELSE 2203 TMAT(I) =-HALF*SMAT(I) 2204 * -HALF*QMAT(INDSQ(I,3)) 2205 ENDIF 2206C 2207 ENDDO 2208! 2209!Sonia: this IF-block is also from old routine.... 2210! 2211 IF (NSYM .GT. 1) THEN 2212 KEND1 = 1 + NCKI(ISYAIJ) 2213 LWRK1 = LWORK - KEND1 2214C 2215 IF (LWRK1 .LT. LENGTH) THEN 2216 CALL QUIT('Out of memory in CCFOP_DENVIR (GATHER-1)') 2217 ENDIF 2218C 2219 CALL CC_GATHER(LENGTH,WORK(KEND1),TMAT,INDSQ(1,6)) 2220 CALL DCOPY(LENGTH,WORK(KEND1),1,TMAT,1) 2221 ENDIF 2222! 2223C------------------------------ 2224C Contract with T2 2225C------------------------------ 2226C 2227! 2228!Sonia: replace following with old version..... 2229! 2230! DO ISYMK = 1,NSYM 2231C 2232! ISYCIJ = MULD2H(ISCKIJ,ISYMK) 2233C 2234! DO ISYMC = 1, NSYM 2235C 2236! ISYMIJ = MULD2H(ISYCIJ,ISYMC) 2237! ISYMAB = MULD2H(ISYMT2,ISYMIJ) 2238! ISYMA = MULD2H(ISYMAB,ISYMB) 2239! ISYMCK = MULD2H(ISYMC,ISYMK) 2240C 2241! NTOTCK = MAX(NT1AM(ISYMCK),1) 2242! NTOTIJ = MAX(NMATIJ(ISYMIJ),1) 2243! NTOTA = MAX(NVIR(ISYMA),1) 2244C 2245! KOFF1 = ISAIK(ISYMAI,ISYMJ) + 1 2246! KOFF2 = ISAIKL(ISYMCK,ISYMIJ) + 1 2247! KOFF3 = ICKATR(ISYMCK,ISYMA) + 1 2248C 2249! CALL DGEMM('N','T',NVIR(ISYMA),NT1AM(ISYMCK), 2250! * NMATIJ(ISYMIJ),TWO,WORK(KOFF1),NTOTA, 2251! * TMAT(KOFF2),NTOTCK,ONE,RINTE1(KOFF3), 2252! * NTOTA) 2253C 2254! ENDDO ! ISYMC 2255! ENDDO ! ISYMK 2256 2257 DO ISYMCK = 1,NSYM 2258C 2259 ISYMIJ = MULD2H(ISCKIJ,ISYMCK) 2260 ISYMAB = MULD2H(ISYMT2,ISYMIJ) 2261 ISYMA = MULD2H(ISYMAB,ISYMB) 2262C 2263 NTOTCK = MAX(NT1AM(ISYMCK),1) 2264 NTOTIJ = MAX(NMATIJ(ISYMIJ),1) 2265 NTOTA = MAX(NVIR(ISYMA),1) 2266C 2267 KOFF1 = ISAIKL(ISYMCK,ISYMIJ) + 1 2268 KOFF2 = IMAIJA(ISYMIJ,ISYMA) + 1 2269 KOFF3 = ICKATR(ISYMCK,ISYMA) + 1 2270C 2271 CALL DGEMM('N','N',NT1AM(ISYMCK),NVIR(ISYMA), 2272 * NMATIJ(ISYMIJ),TWO,TMAT(KOFF1),NTOTCK, 2273 * WORK(KOFF2),NTOTIJ,ONE,RINTE1(KOFF3), 2274 * NTOTCK) 2275C 2276 ENDDO 2277C 2278C------------------------- 2279C Second term. 2280C------------------------- 2281C 2282 DO I = 1,LENGTH 2283C 2284 IF (IOPT .EQ. 1) THEN 2285C 2286 TMAT(I) = TWO*SMAT(INDSQ(I,1)) 2287 * - SMAT(I) 2288 * - SMAT(INDSQ(I,2)) 2289 * + TWO*QMAT(INDSQ(I,2)) 2290 * - QMAT(INDSQ(I,3)) 2291 * - QMAT(INDSQ(I,1)) 2292C 2293 ELSE 2294 TMAT(I) =-HALF*SMAT(INDSQ(I,1)) 2295 * -HALF*QMAT(INDSQ(I,2)) 2296 ENDIF 2297 ENDDO 2298! 2299!Sonia: This IF-block is also from old version.... 2300! 2301 IF (NSYM .GT. 1) THEN 2302 KEND1 = 1 + NCKI(ISYAIJ) 2303 LWRK1 = LWORK - KEND1 2304C 2305 IF (LWRK1 .LT. LENGTH) THEN 2306 CALL QUIT('Out of memory in CCFOP_DENVIR (GATHER-1)') 2307 ENDIF 2308C 2309 CALL CC_GATHER(LENGTH,WORK(KEND1),TMAT,INDSQ(1,6)) 2310 CALL DCOPY(LENGTH,WORK(KEND1),1,TMAT,1) 2311 ENDIF 2312C 2313C------------------------------ 2314C Contract with T2 2315C------------------------------ 2316C 2317! 2318!Sonia: As before, replace following with old version..... 2319! 2320! DO ISYMK = 1,NSYM 2321C 2322! ISYCIJ = MULD2H(ISCKIJ,ISYMK) 2323C 2324! DO ISYMC = 1, NSYM 2325C 2326! ISYMIJ = MULD2H(ISYCIJ,ISYMC) 2327! ISYMAB = MULD2H(ISYMT2,ISYMIJ) 2328! ISYMA = MULD2H(ISYMAB,ISYMB) 2329! ISYMCK = MULD2H(ISYMC,ISYMK) 2330C 2331! NTOTCK = MAX(NT1AM(ISYMCK),1) 2332! NTOTIJ = MAX(NMATIJ(ISYMIJ),1) 2333! NTOTA = MAX(NVIR(ISYMA),1) 2334C 2335! KOFF1 = ISAIK(ISYMAI,ISYMJ) + 1 2336! KOFF2 = ISAIKL(ISYMCK,ISYMIJ) + 1 2337! KOFF3 = ICKATR(ISYMCK,ISYMA) + 1 2338C 2339! CALL DGEMM('N','T',NVIR(ISYMA),NT1AM(ISYMCK), 2340! * NMATIJ(ISYMIJ),TWO,WORK(KOFF1),NTOTA, 2341! * TMAT(KOFF2),NTOTCK,ONE,RINTE2(KOFF3), 2342! * NTOTA) 2343C 2344! ENDDO ! ISYMC 2345! ENDDO ! ISYMK 2346C 2347 DO ISYMCK = 1,NSYM 2348C 2349 ISYMIJ = MULD2H(ISCKIJ,ISYMCK) 2350 ISYMAB = MULD2H(ISYMT2,ISYMIJ) 2351 ISYMA = MULD2H(ISYMAB,ISYMB) 2352C 2353 NTOTCK = MAX(NT1AM(ISYMCK),1) 2354 NTOTIJ = MAX(NMATIJ(ISYMIJ),1) 2355 NTOTA = MAX(NVIR(ISYMA),1) 2356C 2357 KOFF1 = ISAIKL(ISYMCK,ISYMIJ) + 1 2358 KOFF2 = IMAIJA(ISYMIJ,ISYMA) + 1 2359 KOFF3 = ICKATR(ISYMCK,ISYMA) + 1 2360C 2361 CALL DGEMM('N','N',NT1AM(ISYMCK),NVIR(ISYMA), 2362 * NMATIJ(ISYMIJ),TWO,TMAT(KOFF1),NTOTCK, 2363 * WORK(KOFF2),NTOTIJ,ONE,RINTE2(KOFF3), 2364 * NTOTCK) 2365C 2366 ENDDO 2367C 2368 CALL QEXIT('CCFOP_DENVIR_SC') 2369C 2370 RETURN 2371 END 2372C /* Deck ccfop_denocc_sc */ 2373 SUBROUTINE CCFOP_DENOCC_SC(OCC,SMAT,QMAT,TMAT,ISYMIM,T2AM,ISYMT2, 2374 * WORK,LWORK,INDSQ,LENSQ,ISYMB,B, 2375 * ISYMD,D,IOPT) 2376C 2377C Written by Kasper Hald, Fall 2001. 2378C Original version working with symmetry, Sonia Coriani 2002 2379C 2380C Purpose : Calculate the contributions to the t3 and t3-bar 2381C densities d_{iajk} and d_{aijk} respectively. 2382C 2383 IMPLICIT NONE 2384C 2385#include "priunit.h" 2386#include "ccsdsym.h" 2387#include "ccorb.h" 2388C 2389 INTEGER ISYMIM, ISYMT2, LWORK,LENSQ, ISYMB, ISYMD, IOPT 2390 INTEGER INDSQ(LENSQ,6) 2391 INTEGER ISYMBD, ISELJI, ISYELK, ISYIJK, ISYMK, ISYMEL, ISYMIJ 2392 INTEGER NTOTEL, NTOTK, KOFF1, KOFF2, KOFF3, ISYML, ISYME 2393 INTEGER ISYMEK, NTOTIJ, isyres 2394C 2395#if defined (SYS_CRAY) 2396 REAL OCC(*), SMAT(*), QMAT(*), TMAT(*), T2AM(*) 2397 REAL WORK(LWORK), TWO, ONE, HALF 2398#else 2399 DOUBLE PRECISION OCC(*), SMAT(*), QMAT(*), TMAT(*), T2AM(*) 2400 DOUBLE PRECISION WORK(LWORK), TWO, ONE, HALF 2401 double precision xnorm, ddot 2402#endif 2403C 2404 PARAMETER (ONE = 1.0D0, TWO = 2.0D0, HALF = 0.5D0) 2405C 2406 CALL QENTER('CCFOP_DENOCC_SC') 2407C 2408C-------------------------------------- 2409C Symmetries and sanity check. 2410C-------------------------------------- 2411C 2412 IF (IOPT .NE. 1 .AND. IOPT .NE. 2) THEN 2413 CALL QUIT('Wrong IOPT in CCFOP_DENOCC_SC') 2414 ENDIF 2415C 2416! ISYMBD = MULD2H(ISYMB,ISYMD) 2417! ISYELK = MULD2H(ISYMT2,ISYMD) 2418! ISELJI = MULD2H(ISYMIM,ISYMBD) 2419! ISYIJK = MULD2H(ISYELK,ISELJI) 2420! 2421 2422 ISYRES = MULD2H(ISYMT2,ISYMIM) 2423 ISYMBD = MULD2H(ISYMB,ISYMD) 2424 ISYELK = MULD2H(ISYMT2,ISYMD) 2425 ISELJI = MULD2H(ISYMIM,ISYMBD) 2426 ISYIJK = MULD2H(ISYRES,ISYMB) 2427C 2428 IF (LWORK .LT. NCKI(ISYELK)) THEN 2429 CALL QUIT('Not enough memory in CCFOP_DENOCC_SC') 2430 ENDIF 2431C 2432C------------------------------------ 2433C Sort T2 to first term. 2434C------------------------------------ 2435C 2436C DO ISYMK = 1, NSYM 2437C ISYMEL = MULD2H(ISYELK,ISYMK) 2438C DO ISYML = 1, NSYM 2439C ISYME = MULD2H(ISYMEL,ISYML) 2440CC 2441C DO K = 1, NRHF(ISYMK) 2442C DO L = 1, NRHF(ISYML) 2443CC 2444C KOFF1 = IT2SP(ISYELK,ISYMD) 2445C * + NCKI(ISYELK)*(D - 1) 2446C * + ISAIK(ISYMEL,ISYMK) 2447C * + NT1AM(ISYMEL)*(K - 1) 2448C * + IT1AM(ISYME,ISYML) 2449C * + NVIR(ISYME)*(L-1) 2450C * + E 2451CC 2452C KOFF2 = ISAIK(ISYMEL,ISYMK) 2453C * + NT1AM(ISYMEL)*(K - 1) 2454C * + IT1AM(ISYME,ISYML) 2455C * + NVIR(ISYME)*(L-1) 2456C * + E 2457CC 2458C CALL DCOPY(NVIR(ISYME),T2AM(KOFF1),1,WORK(KOFF2),1) 2459CC 2460C ENDDO 2461C ENDDO 2462C ENDDO 2463C ENDDO 2464C 2465C-------------------------------------------- 2466C Contract with S and Q intermediates. 2467C-------------------------------------------- 2468C 2469 DO I = 1, NCKIJ(ISELJI) 2470C 2471 IF (IOPT .EQ. 1) THEN 2472 TMAT(I) = SMAT(INDSQ(I,5)) 2473 * - TWO*SMAT(I) 2474 * + SMAT(INDSQ(I,3)) 2475 * + QMAT(INDSQ(I,4)) 2476 * - TWO*QMAT(INDSQ(I,3)) 2477 * + QMAT(I) 2478 ELSE 2479 TMAT(I) = -HALF*SMAT(I) 2480 * -HALF*QMAT(INDSQ(I,3)) 2481 ENDIF 2482C 2483 ENDDO 2484! 2485!Sonia: following bit comes from old code 2486! 2487 IF (NSYM .GT. 1) THEN 2488 IF (LWORK .LT. NCKIJ(ISELJI)) THEN 2489 CALL QUIT('Out of memory in CCFOP_DENOCC_SC (Gather-1)') 2490 ENDIF 2491C 2492 CALL CC_GATHER(NCKIJ(ISELJI),WORK,TMAT,INDSQ(1,6)) 2493 CALL DCOPY(NCKIJ(ISELJI),WORK,1,TMAT,1) 2494 ENDIF 2495 2496! xnorm = ddot(NCKIJ(ISELJI),TMAT,1,TMAT,1) 2497! write(lupri,*) ' TMAT-1 in FOP_DENOC ', xnorm 2498! 2499!Sonia: end of it 2500! 2501C 2502 DO ISYMK = 1, NSYM 2503 ISYMEL = MULD2H(ISYELK,ISYMK) 2504 ISYMIJ = MULD2H(ISYIJK,ISYMK) 2505C 2506 NTOTEL = MAX(NT1AM(ISYMEL),1) 2507 NTOTK = MAX(NRHF(ISYMK),1) 2508C 2509! 2510!Sonia: removed following and change to old version. 2511! It does not work with symmetry! 2512! 2513! KOFF1 = IT2SP(ISYELK,ISYMD) 2514! * + NCKI(ISYELK)*(D-1) 2515! * + ISAIK(ISYMEL,ISYMK) + 1 2516! KOFF2 = ISAIKL(ISYMEL,ISYMIJ) + 1 2517! KOFF3 = I3OVIR(ISYIJK,ISYMB) 2518! * + NMAIJK(ISYIJK)*(B-1) 2519! * + IMAIJK(ISYMIJ,ISYMK) 2520! * + 1 2521C 2522! CALL DGEMM('T','N',NRHF(ISYMK),NMATIJ(ISYMIJ), 2523! * NT1AM(ISYMEL),TWO,T2AM(KOFF1),NTOTEL, 2524! * TMAT(KOFF2),NTOTEL,ONE,OCC(KOFF3), 2525! * NTOTK) 2526 2527 NTOTIJ = MAX(NMATIJ(ISYMIJ),1) 2528 2529 KOFF1 = ISAIKL(ISYMEL,ISYMIJ) + 1 2530 KOFF2 = IT2SP(ISYELK,ISYMD) 2531 * + NCKI(ISYELK)*(D-1) 2532 * + ISAIK(ISYMEL,ISYMK) + 1 2533 KOFF3 = I3OVIR(ISYIJK,ISYMB) 2534 * + NMAIJK(ISYIJK)*(B-1) 2535 * + IMAIJK(ISYMIJ,ISYMK) 2536 * + 1 2537C 2538 CALL DGEMM('T','N',NMATIJ(ISYMIJ),NRHF(ISYMK), 2539 * NT1AM(ISYMEL),TWO,TMAT(KOFF1),NTOTEL, 2540 * T2AM(KOFF2),NTOTEL,ONE,OCC(KOFF3), 2541 * NTOTIJ) 2542C 2543 ENDDO 2544! xnorm = ddot(NCKIJ(isyres),OCC,1,OCC,1) 2545! write(lupri,*) ' OCC in FOP_DENOC ', xnorm 2546! xnorm = ddot(NT2SQ(isyres),T2AM,1,T2AM,1) 2547! write(lupri,*) ' T2AM in FOP_DENOC ', xnorm 2548C 2549C-------------------------------------------- 2550C Contract with S and Q intermediates. 2551C 1. Build second TMAT first 2552C 2. Contract TMAT with sorted T2 2553C-------------------------------------------- 2554C 2555 DO I = 1, NCKIJ(ISELJI) 2556C 2557 IF (IOPT .EQ. 1) THEN 2558 TMAT(I) = SMAT(INDSQ(I,2)) 2559 * - TWO*SMAT(INDSQ(I,1)) 2560 * + SMAT(INDSQ(I,4)) 2561 * + QMAT(INDSQ(I,1)) 2562 * - TWO*QMAT(INDSQ(I,2)) 2563 * + QMAT(INDSQ(I,5)) 2564 ELSE 2565 TMAT(I) = -HALF*SMAT(INDSQ(I,1)) 2566 * -HALF*QMAT(INDSQ(I,2)) 2567 ENDIF 2568C 2569 ENDDO 2570! 2571!Sonia: following bit from old code 2572! 2573 IF (NSYM .GT. 1) THEN 2574 IF (LWORK .LT. NCKIJ(ISELJI)) THEN 2575 CALL QUIT('Out of memory in CCFOP_DENOCC (Gather-1)') 2576 ENDIF 2577C 2578 CALL CC_GATHER(NCKIJ(ISELJI),WORK,TMAT,INDSQ(1,6)) 2579 CALL DCOPY(NCKIJ(ISELJI),WORK,1,TMAT,1) 2580 ENDIF 2581! 2582!Sonia: end of it 2583! 2584C 2585C------------------------------------ 2586C Sort T2 to second term. 2587C------------------------------------ 2588C 2589 DO ISYMK = 1, NSYM 2590 ISYMEL = MULD2H(ISYELK,ISYMK) 2591 DO ISYML = 1, NSYM 2592 ISYME = MULD2H(ISYMEL,ISYML) 2593 ISYMEK = MULD2H(ISYME,ISYMK) 2594C 2595 DO K = 1, NRHF(ISYMK) 2596 DO L = 1, NRHF(ISYML) 2597C 2598 KOFF1 = IT2SP(ISYELK,ISYMD) 2599 * + NCKI(ISYELK)*(D - 1) 2600 * + ISAIK(ISYMEL,ISYMK) 2601 * + NT1AM(ISYMEL)*(K - 1) 2602 * + IT1AM(ISYME,ISYML) 2603 * + NVIR(ISYME)*(L-1) 2604 * + 1 2605C 2606 KOFF2 = ISAIK(ISYMEK,ISYML) 2607 * + NT1AM(ISYMEK)*(L - 1) 2608 * + IT1AM(ISYME,ISYMK) 2609 * + NVIR(ISYME)*(K-1) 2610 * + 1 2611C 2612 CALL DCOPY(NVIR(ISYME),T2AM(KOFF1),1,WORK(KOFF2),1) 2613C 2614 ENDDO 2615 ENDDO 2616 ENDDO 2617 ENDDO 2618C 2619C 2620 DO ISYMK = 1, NSYM 2621 ISYMEL = MULD2H(ISYELK,ISYMK) 2622 ISYMIJ = MULD2H(ISYIJK,ISYMK) 2623C 2624 NTOTEL = MAX(NT1AM(ISYMEL),1) 2625 NTOTK = MAX(NRHF(ISYMK),1) 2626C 2627! 2628!Sonia: following bit does not work with symmetry! 2629! Replaced with older version 2630! 2631! KOFF1 = ISAIK(ISYMEL,ISYMK) + 1 2632! KOFF2 = ISAIKL(ISYMEL,ISYMIJ) + 1 2633! KOFF3 = I3OVIR(ISYIJK,ISYMB) 2634! * + NMAIJK(ISYIJK)*(B-1) 2635! * + IMAIJK(ISYMIJ,ISYMK) 2636! * + 1 2637C 2638! CALL DGEMM('T','N',NRHF(ISYMK),NMATIJ(ISYMIJ), 2639! * NT1AM(ISYMEL),TWO,WORK(KOFF1),NTOTEL, 2640! * TMAT(KOFF2),NTOTEL,ONE,OCC(KOFF3), 2641! * NTOTK) 2642! 2643 NTOTIJ = MAX(NMATIJ(ISYMIJ),1) 2644 KOFF1 = ISAIKL(ISYMEL,ISYMIJ) + 1 2645 KOFF2 = ISAIK(ISYMEL,ISYMK) + 1 2646 KOFF3 = I3OVIR(ISYIJK,ISYMB) 2647 * + NMAIJK(ISYIJK)*(B-1) 2648 * + IMAIJK(ISYMIJ,ISYMK) 2649 * + 1 2650C 2651 CALL DGEMM('T','N',NMATIJ(ISYMIJ),NRHF(ISYMK), 2652 * NT1AM(ISYMEL),TWO,TMAT(KOFF1),NTOTEL, 2653 * WORK(KOFF2),NTOTEL,ONE,OCC(KOFF3), 2654 * NTOTIJ) 2655 ENDDO 2656! xnorm = ddot(NCKIJ(isyres),OCC,1,OCC,1) 2657! write(lupri,*) ' OCC in FOP_DENOC ', xnorm 2658C 2659C----------------------- 2660C End. 2661C----------------------- 2662C 2663 CALL QEXIT('CCFOP_DENOCC_SC') 2664C 2665 RETURN 2666 END 2667C /* Deck den_biasort */ 2668 SUBROUTINE DEN_BIASORT(VIRREAL,VIRTMP,ISYVIR) 2669C 2670C Written by Kasper Hald, 2001. 2671C 2672C Purpose : Sort the two electron densities d(bia) -> d(aib) 2673C where the densities have a constant C. 2674C 2675 IMPLICIT NONE 2676C 2677#include "priunit.h" 2678#include "ccsdsym.h" 2679#include "ccorb.h" 2680C 2681 INTEGER ISYVIR, ISYMB, ISYMA, ISYMI, ISYMAI 2682 INTEGER KOFF1, KOFF2, ISYMBI 2683COMMENT COMMENT 2684 integer khcount 2685COMMENT COMMENT 2686C 2687#if defined (SYS_CRAY) 2688 REAL VIRREAL(*), VIRTMP(*) 2689#else 2690 DOUBLE PRECISION VIRREAL(*), VIRTMP(*), tmp 2691#endif 2692C 2693 CALL QENTER('DEN_BIASORT') 2694C 2695C---------------------------------------- 2696C Sort matrix. 2697C---------------------------------------- 2698C 2699 DO ISYMB = 1, NSYM 2700 ISYMAI = MULD2H(ISYMB,ISYVIR) 2701 DO ISYMA = 1, NSYM 2702 ISYMI = MULD2H(ISYMAI,ISYMA) 2703 ISYMBI = MULD2H(ISYMI,ISYMB) 2704 DO B = 1, NVIR(ISYMB) 2705 DO I = 1, NRHF(ISYMI) 2706C 2707 KOFF1 = ICKATR(ISYMAI,ISYMB) 2708 * + NT1AM(ISYMAI)*(B-1) 2709 * + IT1AM(ISYMA,ISYMI) 2710 * + NVIR(ISYMA)*(I-1) 2711 * + 1 2712C 2713 KOFF2 = ICKATR(ISYMBI,ISYMA) 2714 * + IT1AM(ISYMB,ISYMI) 2715 * + NVIR(ISYMB)*(I-1) 2716 * + B 2717C 2718 CALL DCOPY(NVIR(ISYMA),VIRREAL(KOFF1),1, 2719 * VIRTMP(KOFF2),NT1AM(ISYMBI)) 2720C 2721 ENDDO 2722 ENDDO 2723 ENDDO 2724 ENDDO 2725C 2726C------------------------------------ 2727C Copy back to original matrix 2728C------------------------------------ 2729C 2730 CALL DCOPY(NCKATR(ISYVIR),VIRTMP(1),1,VIRREAL(1),1) 2731C 2732C----------------------- 2733C End. 2734C----------------------- 2735C 2736 CALL QEXIT('DEN_BIASORT') 2737C 2738 RETURN 2739 END 2740 2741