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_TRDRV 20 SUBROUTINE CC_TRDRV(ECURR,FRHO1,LUFR1,FRHO2,LUFR2,FRHO12,LUFR12, 21 * FC1AM,LUFC1,FC2AM,LUFC2,FC12AM,LUFC12, 22 * TRIPLET,CC3DIIS,ITRAN,FREQ,FS12AM,LUFS12, 23 * FS2AM,LUFS2, 24 * ISIDE,IST,NSIM,WORK,LWORK,APROXR12) 25C 26C------------------------------------------------------------------- 27C 28C "transformation drive" 29C 30C Directs the transformation of trial vecors multiplyed onto the 31C Jacobian from the right and left. 32C 33C Performs NSIM linear transformations. 34C 35C The trial vectors are taken from files: 36C 37C (name,unit nr.) 38C c1am: (FC1AM,LUFC1) singles 39C c2am: (FC2AM,LUFC2) doubles 40C cr12am: (FC12AM,LUFC12) R12 doubles 41C Rho1: (FRHO1,LUFR1) singles 42C Rho2: (FRHO2,LUFR2) doubles 43C RhoR12: (FRHO12,LUFR12) R12 doubles 44C 45C Cholesky CC2 by Thomas Bondo Pedersen, Feb. 2003. 46C - apart from a different trf. routine, the major difference from 47C the viewpoint of the solver (i.e. this driver) is that no doubles 48C are ever stored in core but are assembled on-the-fly. Furthermore, 49C we have to attach a frequency to each trial vector in order to 50C calculate the correct doubles contribution. Thus, it is assumed 51C that all frequencies passed in FREQ are identical (should be checked 52C by calling routine CCEQ_SOL)! 53C 54C It is assumed that on disc are the global intermediates 55C For singlet: 56C CC2: E1, E2, F 57C CCSD & CC3: E1, E2, F, Gamma, BF, C, D intermedias. 58C RCCD, DRCCD: we take a detour/noddy, no global used (FRAN) 59C For triplet: 60C CCSD : E1,E2,F,Gamma,BF,C,D,CD intermediates 61C 62C Local intermediates are put to scratch as well: 63C 64C For right hand transformation: 65C C-tilde,D-tilde intermediates. 66C 67C For left hand side: 68C 69C O3V integrals. 70C 71C For Triplet : (1)C-tilde, (3)C-tilde, (1)D-tilde 72C (P)D-tilde, (M)D-tilde. 73C 74C 75C The variable MINSCR controls how the vectors are to be 76C transformed: one by one or all in one integral calculation. 77C 78C Variables IVEC is the start vector number to be transformed 79C and ITR is the vector number for the result vector. 80C (Used slightly different in CCS and CC2 relative to CCSD,.. 81C since in these cases) 82C 83C Biorthonormal basis: 84C 85C Singlet: 86C (L: 1/2*Eia,(2*EiaEjb + EjaEib)/[6*(1 + delta(ab)delta(ij))] ) 87C (R: Eai, EaiEbj) 88C 89C Triplet: See CC_RHTR3.F 90C 91C Version 1.0 2-11-1996 (base on cclr_trr) 92C Written by Ove Christiansen 2-11-1996. 93C Changes for triplet transformation by Christof Haettig, 1999. 94C Changes for CC-R12 by Christof Haettig, Jun 2003. 95C RCCD and DRCCD handled on their own. MFIozzi (FRAN), Nov 2009 96C 97C------------------------------------------------------------------- 98C 99#include "implicit.h" 100#include "priunit.h" 101#include "maxash.h" 102#include "maxorb.h" 103#include "ccorb.h" 104#include "ccisao.h" 105#include "ccsdsym.h" 106#include "ccsdinp.h" 107#include "ccfield.h" 108#include "cclr.h" 109#include "ccsdio.h" 110#include "leinf.h" 111#include "aovec.h" 112#include "r12int.h" 113Cholesky 114#include "ccdeco.h" 115Cholesky 116 117 PARAMETER ( TWO = 2.0D00,XHALF=0.5D00 ) 118 DIMENSION WORK(LWORK), FREQ(*) 119C 120 LOGICAL ORSAVE,T2TSAV,MSCRS,TRIPLET,CC3DIIS,DEBUGV 121 PARAMETER (DEBUGV = .FALSE.) 122 CHARACTER*8 FRHO1,FRHO2,FC1AM,FC2AM,FR2SD,FRHO12,FC12AM,FS12AM, 123 & FS2AM 124 CHARACTER*3 APROXR12 125 INTEGER ITRAN(NSIM) 126 LOGICAL LOCDBG 127 PARAMETER (LOCDBG = .FALSE.) 128C 129C------------------------- 130C Save and initialize. 131C------------------------- 132C 133 CALL QENTER('CC_TRDRV') 134C 135 MSCRS = MINSCR 136 IF ( CCSDT .OR. (ISIDE.EQ.2) .OR. FDJAC .OR. JACEXP 137 * .OR. MLCC3) THEN 138 MINSCR = .TRUE. 139 ENDIF 140C 141 IF (IPRINT .GT. 5) THEN 142 CALL AROUND(' START OF CC_TRDRV ') 143Chol IF (DIRECT) WRITE(LUPRI,*) ' Atomic direct calculation' 144 IF (CHOINT) THEN 145 WRITE(LUPRI,*) ' Cholesky calculation' 146 ELSE IF (DIRECT) THEN 147 WRITE(LUPRI,*) ' Atomic direct calculation' 148 END IF 149C 150 IF (CCR12) WRITE(LUPRI,*) ' CC-R12 calculation' 151 WRITE(LUPRI,*) ' Workspace input: ',LWORK 152 WRITE(LUPRI,*) ' Triplet flag : ',TRIPLET 153 WRITE(LUPRI,'(/A48,I5)') 154 * ' CC_TRDRV: The total number of trial vectors is: ',NSIM 155 ENDIF 156C 157Cholesky 158 IF (.NOT. (CHOINT .AND. CC2)) THEN 159 T2TSAV = T2TCOR 160 IF (CCS .OR. CC2) T2TCOR = .FALSE. 161 IF (CC2 .AND.(ISIDE .EQ. 2)) T2TCOR = T2TSAV 162 IF (TRIPLET) T2TCOR = .FALSE. 163 END IF 164C 165 NC2 = 1 166 IF ( CCS ) NC2 = 0 167C 168 IF (CHOINT .AND. CC2) THEN 169 NCAMP = NT1AM(ISYMTR) 170 NC2 = 0 171 ELSE 172 NCAMP = NT1AM(ISYMTR) + NT2AM(ISYMTR) 173 IF (CCR12) NCAMP = NCAMP + NTR12AM(ISYMTR) 174 IF (TRIPLET) NCAMP = NCAMP + NT2AMA(ISYMTR) 175 IF (CCS) NCAMP = NT1AM(ISYMTR) 176 END IF 177C 178 IF (CHOINT .AND. JACEXP ) CALL QUIT( 179 * 'CC_TRDRV: JACEXP test flag does not work in CHOINT calc.') 180 IF (DIRECT .AND. JACEXP ) CALL QUIT( 181 * 'CC_TRDRV: JACEXP test flag does not work in DIRECT calc.') 182C 183C------------------------------------------------------------------ 184C Make rho2 file name. 185C For CCSD rho2 has to be stored on different file due to 186C different length. 187C------------------------------------------------------------------ 188C 189 IF (.NOT. (CCS.OR.CC2)) THEN 190 LUFSD = -1 191 FR2SD = 'CC_TRA2_' 192 ELSE 193 LUFSD = LUFR2 194 FR2SD = FRHO2 195 ENDIF 196C 197C------------------------------------------------------------------ 198C Cheat and do CCS left transformation by right hand 199C transformation. 200C------------------------------------------------------------------ 201C 202 NSIDSA = ISIDE 203 IF ((ISIDE .EQ. -1 ) .AND. CCS ) ISIDE = 1 204C 205C-------------------------------------------------------------------- 206C Make total transformation. 207C Orthonormal basis with omegor packed rho in 208C delta loop in ccsd. 209C NB: It is assumed that also omegor for intermediate BF 210C for ccsd. 211C-------------------------------------------------------------------- 212C 213 ONLY21 = .FALSE. 214C 215 ORSAVE = OMEGOR 216C IF ( CCS .OR. CC2) THEN 217C OMEGOR = .FALSE. 218C ELSE 219 OMEGOR = .TRUE. 220C ENDIF 221 222C 223C------------------------------------------------------ 224C Transform vectors: ONE or ALL. 225C If minscr then one at a time, else all vectors in one 226C integral calculation. 227C 228C Keep C1 in core and C2 on disk. 229C Keep rho1 in core and rho2 on disk. ????????????????????? 230C------------------------------------------------------ 231C 232 IF ( MINSCR ) THEN 233 NSIMTR = 1 234 ELSE 235C NSIMTR = NSIM 236 NSIMTR = MAX(NSIM,1) 237 ENDIF 238C 239 IF (NSIMTR .GT. 100 ) CALL QUIT( 240 * ' Maximum nr of simultaneously trial vectors') 241C 242 IF (CC3DIIS .AND. (.NOT.MINSCR)) 243 * CALL QUIT('CC3DIIS option requires MINSCR=.TRUE.') 244C 245C-------------------------------------------------------- 246C Max is 100 due to limitations in the array IT2DLR 247C set in ccsd.cdk. But irrelavant for Cholesky CC2 248C-------------------------------------------------------- 249C 250 NL = (NSIM -1 )/NSIMTR + 1 251C 252 DO 100 I = 1, NL 253C 254 IF ( .NOT.MINSCR) THEN 255C 256 K1 = IST 257C 258 IF (IPRINT .GT. 5 ) THEN 259 IF (CHOINT .AND. CC2) THEN 260 WRITE(LUPRI,'(A24,I3,A37)') 261 * ' CC_TRDRV: Transforming ', 262 * NSIMTR,' vectors in one call to CC_CHOATR.' 263 ELSE 264 WRITE(LUPRI,'(A24,I3,A37)') 265 * ' CC_TRDRV: Transforming ', 266 * NSIMTR,' vectors in one call to CC_RHTR.' 267 END IF 268 ENDIF 269C 270 ELSE 271C 272 IF (CC3DIIS) THEN 273 K1 = ITRAN(I) 274 ECURR = FREQ(K1) 275 ELSE 276 K1 = I + IST - 1 277 END IF 278 279C 280 IF (IPRINT .GT. 5 ) THEN 281 WRITE(LUPRI,'(A35,I4)') 282 * ' CC_TRDRV: Transforming vector nr.:',K1 283 IF (CCSDT) WRITE(LUPRI,'(A,F8.4)') 284 * ' CC_TRDRV: Frequency in R3 denominators:',ECURR 285 ENDIF 286C 287 ENDIF 288C 289 CALL FLSHFO(LUPRI) 290 291C ----------------------------------------------- 292C Prepare memory depending on which case we have: 293C ----------------------------------------------- 294 295 IF (.NOT. TRIPLET) THEN 296Cholesky 297C 298 IF (CHOINT .AND. CC2) THEN 299 300C Cholesky CC2 section. 301C --------------------- 302 303CTODO: Do we actually want to batch over trial vectors for Cholesky CC2 ? 304C..... All it does is take up perfectly good work space for the extensive 305C..... batchings in the transformation routines! 306 307 NRHO1 = NT1AM(ISYMTR)*NSIMTR 308 309 KRHO1 = 1 310 KC1AM = KRHO1 + NRHO1 311 KEND1 = KC1AM + NRHO1 312 LWRK1 = LWORK - KEND1 + 1 313 314 IF (LWRK1 .LT. 0) THEN 315 WRITE(LUPRI,'(//,5X,A)') 316 & 'Insufficient memory in CC_TRDRV' 317 WRITE(LUPRI,'(5X,A,I10)') 318 & 'Number of trial/trf. vectors in core: ',NSIMTR 319 WRITE(LUPRI,'(5X,A,I10,/,5X,A,I10,/)') 320 & 'Need (more than): ',KEND1-1, 321 & 'Available : ',LWORK 322 CALL QUIT('Insufficient memory in CC_TRDRV') 323 ENDIF 324 325C Read trial vectors. 326C ------------------- 327 328 DO IV = 1,NSIMTR 329 KOFF1 = KC1AM + NT1AM(ISYMTR)*(IV - 1) 330 NR1 = IV + K1 - 1 331 CALL CC_RVEC(LUFC1,FC1AM,NT1AM(ISYMTR),NT1AM(ISYMTR), 332 & NR1,WORK(KOFF1)) 333 ENDDO 334 335C Print trial vectors. 336C -------------------- 337 338 IF (IPRINT .GT. 45) THEN 339 KRHO2 = KEND1 ! just a dummy... 340 DO IV = 1, NSIMTR 341 KOFF1 = KC1AM + NT1AM(ISYMTR)*(IV - 1) 342 NR1 = IV + K1 - 1 343 CALL AROUND('CC_TRDRV: C1') 344 WRITE(LUPRI,*) 'Vector nr is = ',NR1 345 CALL CC_PRP(WORK(KOFF1),WORK(KRHO2),ISYMTR,1,NC2) 346 ENDDO 347 ENDIF 348 349 ELSE 350C 351C Conventional 352C 353C ------------ 354C Allocations: 355C ------------ 356 NRHO2 = MAX(NT2AM(ISYMTR),NT2AM(1),2*NT2ORT(ISYMTR)) 357 IF (ISIDE .EQ. -1) NRHO2 = MAX(NRHO2,2*NT2ORT(1)) 358 IF ( CC2 ) NRHO2 = MAX(NT2AM(ISYMTR),NT2AM(1)) 359 IF ( CCS ) NRHO2 = 2 360C 361 NC2AM = MAX(NT2SQ(ISYMTR),NT2SQ(1), 362 * NT2AM(ISYMTR)+2*NT2ORT(1),NT2R12(1)) 363 IF ( CC2 ) NC2AM = MAX(NT2SQ(ISYMTR),NT2SQ(1)) 364 IF ( CCS ) NC2AM = 2 365C 366 NRHO1 = NT1AM(ISYMTR)*NSIMTR 367C 368 KRHO1 = 1 369 KRHO2 = KRHO1 + NRHO1 370 KC1AM = KRHO2 + NRHO2 371 KC2AM = KC1AM + NT1AM(ISYMTR)*NSIMTR 372 KEND1 = KC2AM + NC2AM 373 LWRK1 = LWORK - KEND1 374 IF (LWRK1 .LE. 0 ) 375 * CALL QUIT('Too little workspace in cclr_trr') 376C 377C --------------------------------- 378C Read the CC trial vectors from disk. 379C --------------------------------- 380 IF (IPRINT .GT. 45 .OR. LOCDBG) CALL AROUND('CC_TRDRV: C1') 381 DO 150 IV = 1, NSIMTR 382 KOFF1 = KC1AM + NT1AM(ISYMTR)*(IV - 1) 383 CALL CC_RVEC(LUFC1,FC1AM,NT1AM(ISYMTR),NT1AM(ISYMTR), 384 * IV+K1-1,WORK(KOFF1)) 385 IF (IPRINT.GT.45 .AND. (CCS.OR.(.NOT.MINSCR)) 386 * .OR. LOCDBG) THEN 387 WRITE(LUPRI,*) 'Vector nr.',IV+K1-1 388 CALL CC_PRP(WORK(KOFF1),WORK(KRHO2),ISYMTR,1,0) 389 ENDIF 390 150 CONTINUE 391C 392 IF ((.NOT.CCS).AND.( MINSCR)) THEN 393 CALL CC_RVEC(LUFC2,FC2AM,NT2AM(ISYMTR),NT2AM(ISYMTR), 394 * K1,WORK(KRHO2)) 395 IF (IPRINT .GT. 45 .OR. LOCDBG) THEN 396 CALL AROUND('CC_TRDRV: (C1,C2) packed ') 397 WRITE(LUPRI,*) 'Vector nr.',K1 398 CALL CC_PRP(WORK(KC1AM),WORK(KRHO2),ISYMTR,1,1) 399 ENDIF 400 ENDIF 401C 402C IF (CCR12 .AND. (IPRINT .GT. 45 .OR. LOCDBG)) THEN 403C KRHO12 = KEND1 404C KEND2 = KRHO12 + NTR12AM(ISYMTR) 405C CALL CC_RVEC(LUFC12,FC12AM,NTR12AM(ISYMTR),NTR12AM(ISYMTR), 406C * K1,WORK(KRHO12)) 407C CALL AROUND('CC_TRDRV: R12 trial vector packed ') 408C WRITE(lupri,*) 'Vector nr.',K1 409C CALL CC_PRPR12(WORK(KRHO12),ISYMTR,1,.FALSE.) 410C END IF 411C 412C ----------------------------------------------- 413C Test with finited difference CC jacobian if FDJAC. 414C ----------------------------------------------- 415 IF (FDJAC .OR. JACEXP) THEN 416 IF (CCR12) THEN 417 KC12AM = KEND1 418 KS12AM = KC12AM + NTR12AM(ISYMTR) 419 KRHO12 = KS12AM + NTR12AM(ISYMTR) 420 KEND1 = KRHO12 + NTR12AM(ISYMTR) 421 IF (IANR12.EQ.2) THEN 422 KS2AM = KEND1 423 KEND1 = KS2AM + NT2AM(ISYMTR) 424 END IF 425 LWRK1 = LWORK - KEND1 426 IF (LWRK1 .LE. 0 ) 427 * CALL QUIT('Too little workspace in cclr_trr (2)') 428 429 CALL CC_RVEC(LUFC12,FC12AM,NTR12AM(ISYMTR), 430 * NTR12AM(ISYMTR),K1,WORK(KC12AM)) 431 IF (IPRINT .GT. 45 .OR. LOCDBG) THEN 432 CALL AROUND('R12 double excitation part of vector') 433 CALL OUTPUT(WORK(KRHO12),1,NTR12AM(ISYMTR),1,1, 434 * NTR12AM(ISYMTR),1,1,LUPRI) 435 END IF 436 END IF 437 438 CALL DCOPY(NT2AM(ISYMTR),WORK(KRHO2),1,WORK(KC2AM),1) 439 CALL CCLR_DUMTRR(WORK(KC1AM),WORK(KRHO1),WORK(KRHO2), 440 * WORK(KRHO12),WORK(KS12AM),WORK(KS2AM), 441 * WORK(KEND1),LWRK1) 442C 443C ------------------------------ 444C Print the transformed vectors. 445C ------------------------------ 446 IF (IPRINT .GT. 50 .OR. LOCDBG) THEN 447 CALL AROUND('CC_TRDRV: RHO trans. by DUMTRR .') 448 CALL CC_PRP(WORK(KRHO1),WORK(KRHO2),ISYMTR,1,1) 449 IF (CCR12) THEN 450 WRITE(LUPRI,*) 'A * C_12:' 451 CALL OUTPUT(WORK(KRHO12),1,NTR12AM(ISYMTR),1,1, 452 * NTR12AM(ISYMTR),1,1,LUPRI) 453 IF (IANR12.EQ.2) THEN 454 WRITE(LUPRI,*) 'S * C2:' 455 CALL OUTPUT(WORK(KS2AM),1,NT2AM(ISYMTR),1,1, 456 * NT2AM(ISYMTR),1,1,LUPRI) 457 END IF 458 WRITE(LUPRI,*) 'S * C_12:' 459 CALL OUTPUT(WORK(KS12AM),1,NTR12AM(ISYMTR),1,1, 460 * NTR12AM(ISYMTR),1,1,LUPRI) 461 END IF 462 ENDIF 463 464 CALL CC_WVEC(LUFR1,FRHO1,NT1AM(ISYMTR), 465 * NT1AM(ISYMTR),K1,WORK(KRHO1)) 466 CALL CC_WVEC(LUFR2,FRHO2,NT2AM(ISYMTR), 467 * NT2AM(ISYMTR),K1,WORK(KRHO2)) 468 IF (CCR12) THEN 469 CALL CC_WVEC(LUFR12,FRHO12,NTR12AM(ISYMTR), 470 * NTR12AM(ISYMTR),K1,WORK(KRHO12)) 471 CALL CC_WVEC(LUFS12,FS12AM,NTR12AM(ISYMTR), 472 * NTR12AM(ISYMTR),K1,WORK(KS12AM)) 473 IF (IANR12.EQ.2) THEN 474 CALL CC_WVEC(LUFS2,FS2AM,NT2AM(ISYMTR),NT2AM(ISYMTR), 475 & K1,WORK(KS2AM)) 476 END IF 477 END IF 478 ELSE 479C 480C ---------------------------------------- 481C Prepare the C-amplitudes. 482C ---------------------------------------- 483 IF ( .NOT. CCS ) THEN 484 IF ( ISIDE .GE. 1) THEN 485 CALL CCLR_DIASCL(WORK(KRHO2),TWO,ISYMTR) 486 ENDIF 487 CALL CC_T2SQ(WORK(KRHO2),WORK(KC2AM),ISYMTR) 488 IF (IPRINT.GT.50 .OR. LOCDBG) THEN 489 CALL AROUND('CC_TRDRV: (C1,C2) squared ') 490 CALL CC_PRSQ(WORK(KC1AM),WORK(KC2AM),ISYMTR,1,1) 491 ENDIF 492 ENDIF 493C 494C ---------------- 495C Zero rho vector. 496C ---------------- 497 CALL DZERO(WORK(KRHO1),NT1AM(ISYMTR)*NSIMTR) 498 CALL DZERO(WORK(KRHO2),NRHO2) 499C 500C ------------------- 501C File read and save. 502C ------------------ 503 IF (.NOT. (CCS.OR.CC2)) THEN 504 CALL WOPEN2(LUFSD,FR2SD,64,0) 505 ENDIF 506C 507 NRHO2 = MAX(NT2AM(ISYMTR),2*NT2ORT(ISYMTR)) 508 IF (CC2 ) NRHO2 = NT2AM(ISYMTR) 509 DO 80 IV = 1, NSIMTR 510 NR1 = IV + K1 - 1 511 IF (.NOT. (CCS.OR.CC2)) THEN 512 NR2 = IV 513 ELSE 514 NR2 = NR1 515 ENDIF 516 CALL CC_WVEC(LUFR1,FRHO1,NT1AM(ISYMTR), 517 * NT1AM(ISYMTR),NR1,WORK(KRHO1)) 518 IF (.NOT.CCS) THEN 519 CALL CC_WVEC(LUFSD,FR2SD,NRHO2,NRHO2,NR2, 520 * WORK(KRHO2)) 521 ENDIF 522 80 CONTINUE 523C 524 ENDIF 525C 526 END IF 527C 528 ELSE 529C 530C ----------------------------- 531C For Triplet start at the 532C beginning of the workspace. 533C ----------------------------- 534C 535 NRHO1 = NT1AM(ISYMTR)*NSIMTR 536 NRHO2 = NT2AM(ISYMTR)+NT2AMA(ISYMTR) 537 IF (ISIDE .EQ. -1) THEN 538 NRHO2 = MAX(NT2SQ(ISYMTR),NT2ORT(ISYMTR)+NT2ORT3(ISYMTR), 539 & 2*NT2ORT(1)) 540 LRHO1 = NT1AM(ISYMTR) 541 NC2AM = NT2SQ(ISYMTR) 542 ELSE 543 LRHO1 = 0 544 NC2AM = 0 545 END IF 546C 547 KRHO1 = 1 548 KRHO2 = KRHO1 + NRHO1 549 KC1AM = KRHO2 + NRHO2 550 KC2AM = KC1AM + LRHO1*NSIMTR 551 KEND1 = KC2AM + NC2AM 552 LWRK1 = LWORK - KEND1 553 IF (LWRK1 .LE. 0) THEN 554 CALL QUIT('Too little workspace in CC_TRDRV ') 555 ENDIF 556C----------------------------------------------- 557C The left transform need C1 in memory 558C----------------------------------------------- 559C 560 IF (ISIDE .EQ.-1) THEN 561 DO IV = 1, NSIMTR 562 KOFF1 = KC1AM + NT1AM(ISYMTR)*(IV - 1) 563 CALL CC_RVEC(LUFC1,FC1AM,NT1AM(ISYMTR),NT1AM(ISYMTR), 564 * IV+K1-1,WORK(KOFF1)) 565 IF (IPRINT.GT.45 .AND. (CCS.OR.(.NOT.MINSCR)) 566 * .OR. LOCDBG) THEN 567 WRITE(LUPRI,*) 'Vector nr.',IV+K1-1 568 CALL CC_PRP(WORK(KOFF1),WORK(KRHO2),ISYMTR,1,0) 569 ENDIF 570 END DO 571 END IF 572C 573C ------------------------------------- 574C File opening for the triplet case 575C ------------------------------------- 576 IF (.NOT. (CCS.OR.CC2)) THEN 577 CALL WOPEN2(LUFSD,FR2SD,64,0) 578 ENDIF 579 ENDIF 580C 581C---------------------------------------- 582C Calculate transformed vectors. 583C----------------------------------------- 584C 585 IF (.NOT. (FDJAC .OR. JACEXP)) THEN 586 IVEC = K1 587 IF ( .NOT.(CCS.OR.CC2) ) THEN 588 ITR = 1 589 ELSE 590 ITR = K1 591 ENDIF 592C 593 LRHO1 = NT1AM(ISYMTR) 594Cholesky 595C 596C 597 IF (CHOINT .AND. CC2) THEN 598 599C Cholesky CC2 section. 600C --------------------- 601 602 IF (TRIPLET) THEN 603 CALL QUIT('CC_TRDRV: Triplet Cholesky not implemented') 604 ELSE 605 IF (CHEXDI) THEN 606c DO III = 1,NSIMTR 607c write(lupri,*) 'cc_trdrv. freq :',iii,freq(iii) 608c FREQ(III) = ECURR 609c END DO 610c CALL CC_CHOATR(WORK(KRHO1),WORK(KC1AM),FREQ, 611 CALL CC_CHOATR(WORK(KRHO1),WORK(KC1AM),ecurr, 612 & WORK(KEND1),LWRK1,ISYMTR,NSIMTR,ISIDE) 613 ELSE 614 CALL CC_CHOATR(WORK(KRHO1),WORK(KC1AM),FREQ, 615 & WORK(KEND1),LWRK1,ISYMTR,NSIMTR,ISIDE) 616 END IF 617 ENDIF 618 619 GOTO 1234 620 621 END IF 622C 623C Conventional section. 624C --------------------- 625C 626 IF (TRIPLET) THEN 627 IF (CCR12) CALL QUIT('No triplet yet for CCR12') 628 IF (ISIDE .EQ. 1) THEN 629 CALL CC_RHTR3(ECURR,FRHO1,LUFR1,FR2SD,LUFSD, 630 * FC1AM,LUFC1,FC2AM,LUFC2, 631 * WORK,LWORK,NSIMTR, 632 * IVEC,ITR) 633 ELSE IF (ISIDE .EQ. -1) THEN 634 CALL CC_LHTR3(ECURR, 635 * FRHO1,LUFR1,FR2SD,LUFSD, 636 * FC1AM,LUFC1,FC2AM,LUFC2, 637 * WORK(KRHO1),WORK(KRHO2), 638 * WORK(KC1AM),WORK(KC2AM), 639 * WORK(KEND1),LWRK1,NSIMTR, 640 * IVEC,ITR,LRHO1) 641 ELSE 642 CALL QUIT(' ISIDE should be -1 or +1 ') 643 END IF 644C 645 ELSE 646 IF (ISIDE .EQ. 1) THEN 647C 648 CALL CC_RHTR(ECURR, 649 * FRHO1,LUFR1,FR2SD,LUFSD,FRHO12,LUFR12, 650 * FC1AM,LUFC1,FC2AM,LUFC2,FC12AM,LUFC12, 651 * WORK(KRHO1),WORK(KRHO2), 652 * WORK(KC1AM),WORK(KC2AM), 653 * WORK(KEND1),LWRK1,NSIMTR, 654 * IVEC,ITR,LRHO1,.FALSE.,DUMMY,APROXR12) 655C 656 657 IF (DEBUGV) THEN 658 WRITE(LUPRI,*)'analytical RHO1 and RHO2:' 659 CALL CC_PRP(WORK(KRHO1),WORK(KRHO2),1,1,1) 660 ! calculate V bar numerically 661 CALL CC_R12FDVINT(WORK(KC1AM),WORK(KC2AM), 662 & WORK(KEND1),LWRK1,APROXR12, 663 & FC12AM,LUFC12,IVEC) 664c CALL QUIT('DEBUG V BAR') 665 END IF 666 667 ELSE IF (ISIDE .EQ. -1) THEN 668 if ((RCCD).or.(DRCCD)) then 669 !write(lupri,*)'FRAN: call noddy for LHT in RCCD' 670 call flshfo(lupri) 671 call cc_lhtr_rccd(ECURR, 672 & FRHO1,LUFR1,FR2SD,LUFSD, 673 & FC1AM,LUFC1,FC2AM,LUFC2, 674 & WORK(KRHO1),WORK(KRHO2), 675 & WORK(KC1AM),WORK(KC2AM), 676 & WORK(KEND1),LWRK1,NSIMTR, 677 & IVEC,ITR,LRHO1) 678 ! write(lupri,*)'FRAN: out of noddy for LHT in RCCD' 679 call flshfo(lupri) 680 else 681 682 CALL CC_LHTR(ECURR, 683 * FRHO1,LUFR1,FR2SD,LUFSD,FRHO12,LUFR12, 684 * FC1AM,LUFC1,FC2AM,LUFC2,FC12AM,LUFC12, 685 * WORK(KRHO1),WORK(KRHO2), 686 * WORK(KC1AM),WORK(KC2AM), 687 * WORK(KEND1),LWRK1,NSIMTR, 688 * IVEC,ITR,LRHO1,APROXR12) 689 end if 690 ELSE 691 CALL QUIT(' ISIDE should be -1 or +1 ') 692 ENDIF 693 694C IF ( CCR12 ) THEN 695C DO IV = 1, NSIMTR 696C IF (ISIDE .EQ. 1) THEN 697C CALL CC_R12METRIC(ISYMTR,BRASCL,KETSCL, 698C * WORK(KEND1),LWRK1,FC2AM,LUFC2, 699C * FC12AM,LUFC12,FS12AM,LUFS12, 700C * FS2AM,LUFS2,IVEC-1+IV,.FALSE.,DUMMY) 701C ELSE IF (ISIDE .EQ. -1) THEN 702C CALL CC_R12METRIC(ISYMTR,0.5D0*KETSCL,2.0D0*BRASCL, 703C * WORK(KEND1),LWRK1,FC2AM,LUFC2, 704C * FC12AM,LUFC12,FS12AM,LUFS12, 705C * FS2AM,LUFS2,IVEC-1+IV,.FALSE.,DUMMY) 706C END IF 707C END DO 708C END IF 709 710C 711 ENDIF 712C 713 1234 CONTINUE ! From Cholesky section 714C 715C ------------------------------ 716C Print the transformed vectors. 717C ------------------------------ 718Cholesky 719C 720 IF (CHOINT .AND. CC2) THEN 721 722C Cholesky section: print and save trf. vectors. 723C ---------------------------------------------- 724 725 IF (IPRINT .GT. 45) THEN 726 KRHO2 = KEND1 ! just a dummy... 727 DO IV = 1,NSIMTR 728 KOFF1 = KRHO1 + NT1AM(ISYMTR)*(IV - 1) 729 NR1 = IV + K1 - 1 730 CALL AROUND('CC_TRDRV: RHO = trans. Vector ') 731 WRITE (LUPRI,*) 'number of vector on file:',NR1 732 CALL CC_PRP(WORK(KOFF1),WORK(KRHO2),ISYMTR,1,NC2) 733 ENDDO 734 ENDIF 735 736 DO IV = 1,NSIMTR 737 KOFF1 = KRHO1 + NT1AM(ISYMTR)*(IV - 1) 738 NR1 = IV + K1 - 1 739 CALL CC_WVEC(LUFR1,FRHO1,NT1AM(ISYMTR),NT1AM(ISYMTR), 740 & NR1,WORK(KOFF1)) 741 ENDDO 742 743 GOTO 999 744 745 END IF 746C 747C Conventional 748C 749 IF (( IPRINT .GT. 15).OR.(.NOT.(CCS.OR.CC2)).OR.LOCDBG) THEN 750 IF (TRIPLET) THEN 751 NRHO2 = NT2AM(ISYMTR) + NT2AMA(ISYMTR) 752 ELSE 753 NRHO2 = MAX(NT2AM(ISYMTR),2*NT2ORT(ISYMTR)) 754 IF (CC2 ) NRHO2 = NT2AM(ISYMTR) 755 END IF 756 DO 90 IV = 1, NSIMTR 757 NR1 = IV + K1 - 1 758 IF (.NOT. (CCS.OR.CC2)) THEN 759 NR2 = IV 760 ELSE 761 NR2 = NR1 762 ENDIF 763 CALL CC_RVEC(LUFR1,FRHO1,NT1AM(ISYMTR),NT1AM(ISYMTR), 764 * NR1,WORK(KRHO1)) 765 IF (.NOT.CCS) THEN 766 CALL CC_RVEC(LUFSD,FR2SD,NRHO2,NRHO2,NR2,WORK(KRHO2)) 767 END IF 768 IF (CCR12) THEN 769 KRHO12 = KEND1 770 KEND2 = KRHO12 + NTR12AM(ISYMTR) 771 CALL CC_RVEC(LUFR12,FRHO12,NTR12AM(ISYMTR), 772 * NTR12AM(ISYMTR),NR2,WORK(KRHO12)) 773 END IF 774 775 IF (IPRINT .GT. 45 .OR. LOCDBG) THEN 776 CALL AROUND('CC_TRDRV: RHO = trans. Vector ') 777 WRITE (LUPRI,*) 'number of vector on file:',NR2 778 IF (TRIPLET.AND. (.NOT.CCS)) NC2 = 2 779 CALL CC_PRP(WORK(KRHO1),WORK(KRHO2),ISYMTR,1,NC2) 780 IF (CCR12) THEN 781 CALL CC_PRPR12(WORK(KRHO12),ISYMTR,1,.TRUE.) 782 ENDIF 783 ENDIF 784 785 IF (.NOT.(CCS.OR.CC2) .AND. (.NOT. TRIPLET)) THEN 786 CALL CC_WVEC(LUFR2,FRHO2,NT2AM(ISYMTR), 787 * NT2AM(ISYMTR),NR1,WORK(KRHO2)) 788 ENDIF 789 790 IF ((TRIPLET) .AND. (.NOT. (CCS.OR.CC2))) THEN 791 CALL CC_WVEC3(LUFR2,FRHO2, 792 * NT2AM(ISYMTR)+NT2AMA(ISYMTR), 793 * NT2AM(ISYMTR)+NT2AMA(ISYMTR), 794 * NR1,0,WORK(KRHO2)) 795 ENDIF 796 CALL FLSHFO(LUPRI) 797 90 CONTINUE 798 ENDIF 799C 800 999 CONTINUE ! From Cholesky section 801C 802C----------------------- 803C Close files. 804C----------------------- 805C 806 IF ( .NOT.(CCS.OR.CC2)) THEN 807 CALL WCLOSE2(LUFSD,FR2SD,'DELETE') 808 ENDIF 809C 810 ENDIF 811C 812 100 CONTINUE 813C 814C------------- 815C Restore. 816C------------- 817C 818 T2TCOR = T2TSAV 819 ISIDE = NSIDSA 820 OMEGOR = ORSAVE 821 MINSCR = MSCRS 822C 823 IF (IPRINT .GT. 10 .OR. LOCDBG) THEN 824 CALL AROUND(' END OF CC_TRDRV ') 825 CALL FLSHFO(LUPRI) 826 ENDIF 827C 828 1 FORMAT(1x,A35,1X,E20.10) 829 CALL QEXIT('CC_TRDRV') 830 RETURN 831 END 832