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 19*=====================================================================* 20 SUBROUTINE CC_FDFBMAT(CMOQ,LISTL,IDLSTL,LISTR,IDLSTR, 21 & RHODIF,WORK,LWORK) 22C 23C--------------------------------------------------------------------- 24C Test routine for calculating the orbital relaxation contribution 25C to the CC F^B T^A transformation by finite difference on the usual 26C F matrix transformation wrt the CMO coefficients. 27C We pass the derivative C^(1) MO = CMOQ 28C 29C S. Coriani & Ch. Haettig, march 1999 30C--------------------------------------------------------------------- 31C 32#include "implicit.h" 33#include "priunit.h" 34#include "maxorb.h" 35#include "iratdef.h" 36#include "ccorb.h" 37#include "aovec.h" 38#include "ccsdinp.h" 39#include "cclr.h" 40#include "ccsdsym.h" 41#include "ccsdio.h" 42#include "leinf.h" 43#include "ccinftap.h" 44#include "dummy.h" 45C 46C------------------------------------------------------------ 47C the displacement for the finite difference calculation: 48C------------------------------------------------------------ 49 PARAMETER (DELTA = 1.0D-07, DELINV = 1.0D+7) 50C------------------------------------------------------------ 51 DIMENSION WORK(LWORK), RHODIF(*), CMOQ(*) 52 PARAMETER (ONE = 1.0d0, ZERO = 0.0d0, TWO = 2.0d0, XHALF = 0.5D0) 53 CHARACTER MODEL*10, LISTL*(*), LISTR*(*) 54 LOGICAL LHTF 55C 56 IF (IPRINT.GT.10) THEN 57 CALL AROUND('IN CC_FDFB : MAKING FIN. DIFF. CC F*T^A Vector') 58 END IF 59C 60 IF (CCR12) CALL QUIT('Finite-difference of F*T^A Vector for '// 61 & 'CCR12 not adapted') 62C 63C---------------------------- 64C Work space allocations. 65C---------------------------- 66C 67 ISYMTR = 1 68 ISYMOP = 1 69C 70 NTAMP = NT1AMX + NT2AMX 71 IF (CCS) NTAMP = NT1AMX 72C 73 KCMOREF = 1 74 KCMO = KCMOREF + NLAMDS 75 KRHREF1 = KCMO + NLAMDS 76 KRHREF2 = KRHREF1 + NT1AMX 77 KEND0 = KRHREF2 + NT2AMX 78 LWRK0 = LWORK - KEND0 79C 80 KT1AMP0 = KEND0 81 KOMEGA1 = KT1AMP0 + NT1AMX 82 KOMEGA2 = KOMEGA1 + NT1AMX 83 KT2AMP0 = KOMEGA2 + MAX(NT2AMX,2*NT2ORT(1),NT2AO(1)) 84 KSCR2 = KT2AMP0 + NT2AMX 85 KEND1 = KSCR2 + NT2AMX + NT1AMX 86 LWRK1 = LWORK - KEND1 87C 88 KRHO1 = KEND0 89 KRHO2 = KRHO1 + NT1AMX 90 KEND1A = KRHO2 + NT2AMX 91 LWRK1A = LWORK - KEND1A 92C 93 IF ( LWRK1.LT.0 .OR. LWRK1A.LT.0) THEN 94 WRITE(LUPRI,*) 'Too little work space in CC_FDFBMAT ' 95 WRITE(LUPRI,*) 'AVAILABLE: LWORK = ',LWORK 96 WRITE(LUPRI,*) 'NEEDED (AT LEAST) = ',MAX(KEND1,KEND1A) 97 CALL QUIT('TOO LITTLE WORKSPACE IN CC_FDFBMAT ') 98 ENDIF 99C 100C--------------------- 101C Initializations. 102C--------------------- 103C 104 CALL DZERO(RHODIF,NTAMP) 105C 106 IPRSAVE = IPRINT 107C 108C--------------------------------------------- 109C Read the reference CMO vector from disk: 110C--------------------------------------------- 111C 112 LUSIFC = -1 113 CALL GPOPEN(LUSIFC,'SIRIFC','UNKNOWN',' ','UNFORMATTED', 114 * IDUMMY,.FALSE.) 115 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 116 READ(LUSIFC) 117 READ(LUSIFC) 118 READ(LUSIFC) (WORK(KCMOREF-1+I),I=1,NLAMDS) 119 CALL GPCLOSE(LUSIFC,'KEEP') 120C 121C------------------------------------------------------------------- 122C make sure that the correct response intermediates are on disc: 123C------------------------------------------------------------------- 124C 125 CALL DZERO(WORK(KT1AMP0),NT1AMX) 126 LHTF = .FALSE. 127 CALL CCSD_IAJB(WORK(KT2AMP0),WORK(KT1AMP0),LHTF, 128 & .FALSE.,.FALSE.,WORK(KEND1),LWRK1) 129 REWIND(LUIAJB) 130 CALL WRITI(LUIAJB,IRAT*NT2AMX,WORK(KT2AMP0)) 131C 132 IOPT = 3 133 CALL CC_RDRSP('R0 ',0,1,IOPT,MODEL,WORK(KT1AMP0),WORK(KT2AMP0)) 134C 135 RSPIM = .TRUE. 136 CALL CCRHSN(WORK(KOMEGA1),WORK(KOMEGA2), 137 & WORK(KT1AMP0),WORK(KT2AMP0), 138 & WORK(KEND1),LWRK1,'XXX') 139C 140C------------------------------------------------------------------- 141C calculate the reference vector: 142C------------------------------------------------------------------- 143C use short-cut routine (for one single transformation) 144C 145 IPRINT = 0 146 CALL CC_FTRAN(LISTL, IDLSTL, LISTR, IDLSTR, 147 & WORK(KRHO1), LWRK0) 148C 149 IF (CCS) CALL DZERO(WORK(KRHO2),NT2AMX) 150C 151 RHO1N = DDOT(NT1AMX,WORK(KRHO1),1,WORK(KRHO1),1) 152 IF (.NOT.CCS) RHO2N = DDOT(NT2AMX,WORK(KRHO2),1,WORK(KRHO2),1) 153C 154 IF (IPRSAVE.GT.10) THEN 155 WRITE (LUPRI,*) 'CC_FDFB: norm of reference RHO1:',RHO1N 156 IF (.NOT.CCS) 157 * WRITE (LUPRI,*) 'CC_FDFB: norm of reference RHO2:',RHO2N 158 END IF 159C 160c save result on Reference vector 161C 162 CALL DCOPY(NT1AMX,WORK(KRHO1),1,WORK(KRHREF1),1) 163 IF (.NOT.CCS) CALL DCOPY(NT2AMX,WORK(KRHO2),1,WORK(KRHREF2),1) 164 165*----------------------------------------------------------------------* 166* Calculate the derivative of the vector function with respect 167* to the CMO vector by finite difference: 168*----------------------------------------------------------------------* 169 170 DO IDXDIF = 1, NLAMDS 171C 172C ------------------------------- 173C add finite displacement to CMO: 174C ------------------------------- 175C 176 CALL DCOPY(NLAMDS,WORK(KCMOREF),1,WORK(KCMO),1) 177 WORK(KCMO-1 + IDXDIF) = WORK(KCMO-1 + IDXDIF) + DELTA 178C 179 CALL GPOPEN(LUSIFC,'SIRIFC','UNKNOWN',' ','UNFORMATTED', 180 * IDUMMY,.FALSE.) 181 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 182 READ(LUSIFC) 183 READ(LUSIFC) 184 WRITE(LUSIFC) (WORK(KCMO-1+I),I=1,NLAMDS) 185 CALL GPCLOSE(LUSIFC,'KEEP') 186C 187C ------------------------------------- 188C calculate new response (global) intermediates: 189C ------------------------------------- 190C 191 CALL DZERO(WORK(KT1AMP0),NT1AMX) 192 LHTF = .FALSE. 193 CALL CCSD_IAJB(WORK(KT2AMP0),WORK(KT1AMP0), 194 & LHTF,.FALSE.,.FALSE.,WORK(KEND1),LWRK1) 195 REWIND(LUIAJB) 196 CALL WRITI(LUIAJB,IRAT*NT2AMX,WORK(KT2AMP0)) 197C 198 IOPT = 3 199 CALL CC_RDRSP('R0 ',0,1,IOPT,MODEL,WORK(KT1AMP0),WORK(KT2AMP0)) 200C 201 RSPIM = .TRUE. 202 CALL CCRHSN(WORK(KRHO1),WORK(KRHO2), 203 & WORK(KT1AMP0),WORK(KT2AMP0), 204 & WORK(KEND1),LWRK1,'XXX') 205C 206C --------------------------------- 207C calculate the transformed vector: 208C --------------------------------- 209C 210 IPRINT = 0 211 CALL CC_FTRAN(LISTL, IDLSTL, LISTR, IDLSTR, 212 & WORK(KRHO1), LWRK0) 213C 214 IF (CCS) CALL DZERO(WORK(KRHO2),NT2AMX) 215C 216C -------------------------------------------------- 217C construct the row nb. IDXDIF of the result matrix: 218C -------------------------------------------------- 219C 220 RHO1N = DDOT(NT1AMX,WORK(KRHO1),1,WORK(KRHO1),1) 221 IF (.NOT.CCS) RHO2N = DDOT(NT2AMX,WORK(KRHO2),1,WORK(KRHO2),1) 222C 223 IF (IPRSAVE.GT.10) THEN 224 WRITE (LUPRI,*) 'CMO index:',IDXDIF 225 WRITE (LUPRI,*) 'Norm of RHO1: ',RHO1N 226 WRITE (LUPRI,*) 'Norm of RHO2: ',RHO2N 227 END IF 228C 229c Construct {[RES(delta)-RES(ref)]/delta}*C(1) 230C 231 CALL DAXPY(NT1AMX,-1.0D0,WORK(KRHREF1),1,WORK(KRHO1),1) ! 232 IF (.NOT.CCS) 233 & CALL DAXPY(NT2AMX,-1.0D0,WORK(KRHREF2),1,WORK(KRHO2),1) 234 CALL DSCAL(NT1AMX,DELINV,WORK(KRHO1),1) 235 IF (.NOT.CCS) 236 & CALL DSCAL(NT2AMX,DELINV,WORK(KRHO2),1) 237c 238 CALL DAXPY(NT1AMX,CMOQ(IDXDIF),WORK(KRHO1),1, 239 & RHODIF(1),1) 240 IF (.NOT.CCS) CALL DAXPY(NT2AMX,CMOQ(IDXDIF),WORK(KRHO2),1, 241 & RHODIF(1+NT1AMX),1) 242C 243 IF (IPRSAVE.GT.100) THEN 244 CALL CC_PRP(WORK(KRHO1),WORK(KRHO2),1,1,1) 245 WRITE (LUPRI,*) 'CMOQ(index) = ', CMOQ(IDXDIF) 246 WRITE (LUPRI,*) 'accumulated result vector:' 247 CALL CC_PRP(RHODIF,RHODIF(NT1AMX+1),1,1,1) 248 ENDIF 249C 250 END DO 251C 252C---------------------------------------------- 253C restore the reference CMO vector on disk: 254C---------------------------------------------- 255C 256 CALL GPOPEN(LUSIFC,'SIRIFC','UNKNOWN',' ','UNFORMATTED', 257 * IDUMMY,.FALSE.) 258 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 259 READ(LUSIFC) 260 READ(LUSIFC) 261 WRITE(LUSIFC) (WORK(KCMOREF-1+I),I=1,NLAMDS) 262 CALL GPCLOSE(LUSIFC,'KEEP') 263C 264C------------------------------------------------------------------ 265C make sure that all intermediates on file are calculated with 266C the reference CMO vector: 267C------------------------------------------------------------------ 268C 269 CALL DZERO(WORK(KT1AMP0),NT1AMX) 270 LHTF = .FALSE. 271 CALL CCSD_IAJB(WORK(KT2AMP0),WORK(KT1AMP0),LHTF, 272 & .FALSE.,.FALSE.,WORK(KEND1),LWRK1) 273 REWIND(LUIAJB) 274 CALL WRITI(LUIAJB,IRAT*NT2AMX,WORK(KT2AMP0)) 275C 276 IOPT = 3 277 CALL CC_RDRSP('R0 ',0,1,IOPT,MODEL,WORK(KT1AMP0),WORK(KT2AMP0)) 278C 279 RSPIM = .TRUE. 280 CALL CCRHSN(WORK(KOMEGA1),WORK(KOMEGA2), 281 & WORK(KT1AMP0),WORK(KT2AMP0), 282 & WORK(KEND1),LWRK1,'XXX') 283C 284 IPRINT = IPRSAVE 285C 286 IF (IPRINT.GT.10) THEN 287 CALL AROUND(' END OF CC_FDFBMAT ') 288 END IF 289C 290 RETURN 291 END 292*=====================================================================* 293 SUBROUTINE CC_FDFBMAT2(LISTR,IDLSTR,RHODIF1,RHODIF2, 294 & LABEL,IRELAX,WORK,LWORK) 295C 296C--------------------------------------------------------------------- 297C Test routine for calculating a generalized CC F{O} matrix transformed 298C vector by finite difference on the Eta{O} vector 299C 300C Christof Haettig, march 1999 301C--------------------------------------------------------------------- 302C 303#include "implicit.h" 304#include "priunit.h" 305#include "maxorb.h" 306#include "iratdef.h" 307#include "ccorb.h" 308#include "aovec.h" 309#include "ccsdinp.h" 310#include "cclr.h" 311#include "ccsdsym.h" 312#include "ccsdio.h" 313#include "leinf.h" 314#include "dummy.h" 315#include "ccfro.h" 316#include "ccroper.h" 317#include "cclists.h" 318C 319C------------------------------------------------------------ 320C the displacement for the finite difference calculation: 321C------------------------------------------------------------ 322 PARAMETER (DELTA = 1.0D-07, DELINV = 1.0D+07) 323C------------------------------------------------------------ 324C 325 DIMENSION WORK(LWORK) 326 DIMENSION RHODIF1(*), RHODIF2(*) 327 CHARACTER*(3) LISTL, LISTR 328 CHARACTER*(8) FILXI, FILETA 329 CHARACTER MODEL*(10) 330 INTEGER IXETRAN(MXDIM_XEVEC,3) 331C LOGICAL LORX 332C 333 INTEGER IR1TAMP, IRHSR1, IETA1, IROPER, ILSTSYM 334C 335 PARAMETER (ONE=1.0d0, ZERO=0.0d0, TWO=2.0d0, HALF=0.5d0) 336C 337 IF (IPRINT.GT.5) THEN 338 CALL AROUND('in CC_FDFBMAT2: making fini. diff. FB Matrix') 339 ENDIF 340C 341C---------------------------------------------- 342C set up IXETRAN array: 343C---------------------------------------------- 344C 345C IRELAX = 0 346C IF (LORX) THEN 347C IRELAX = IR1TAMP(LABEL,LORX,FREQ,ISYHOP) 348C END IF 349 350 ! set the IXETRAN array for one (XI,ETA) pair 351 IXETRAN(1,1) = IROPER(LABEL,ISYHOP) 352 IXETRAN(2,1) = 0 353 IXETRAN(3,1) = -1 354 IXETRAN(4,1) = 0 355 IXETRAN(5,1) = IRELAX 356C 357 LISTL = 'L0 ' 358 FILXI = 'CCFDFBO1' 359 FILETA = 'CCFDFBX1' 360 IOPTRES = 0 361 NXETRAN = 1 362C 363C---------------------------------------------- 364C initializations & work space allocations. 365C---------------------------------------------- 366C 367 CALL DZERO(RHODIF1,NT1AMX) 368 CALL DZERO(RHODIF2,NT2AMX) 369C 370 IPRSAV = IPRINT 371 IPRINT = 0 372C 373 MODEL = 'UNKNOWN' 374 IF (CCS) MODEL = 'CCS' 375 IF (CC2) MODEL = 'CC2' 376 IF (CCSD) MODEL = 'CCSD' 377C 378 KT1AMPSAV = 1 379 KT2AMPSAV = KT1AMPSAV + NT1AMX 380 KT1AMPA = KT2AMPSAV + NT2AMX 381 KT2AMPA = KT1AMPA + NT1AMX 382 KEND0 = KT2AMPA + NT2AMX 383 LWRK0 = LWORK - KEND0 384C 385 KT0AMP0 = KEND0 386 KT1AMP0 = KT0AMP0 + 2*NALLAI(1) 387 KOMEGA1 = KT1AMP0 + NT1AMX 388 KOMEGA2 = KOMEGA1 + NT1AMX 389 KT2AMP0 = KOMEGA2 + MAX(NT2AMX,2*NT2ORT(1),NT2AO(1)) 390 KSCR2 = KT2AMP0 + NT2AMX 391 KEND1A = KSCR2 + NT2AMX + NT1AMX 392 LWRK1A = LWORK - KEND1A 393C 394 KRHO1 = KEND0 395 KRHO2 = KRHO1 + NT1AMX 396 KEND1B = KRHO2 + NT2AMX 397 LWRK1B = LWORK - KEND1B 398C 399 IF (LWRK1A.LT.0 .OR. LWRK1B.LT.0) THEN 400 WRITE(LUPRI,*) 'TOO LITTLE WORK SPACE IN CC_FDFBMAT2::' 401 WRITE(LUPRI,*) 'AVAILABLE: LWORK = ',LWORK 402 WRITE(LUPRI,*) 'NEEDED (AT LEAST) = ',MAX(KEND1A,KEND1B) 403 CALL QUIT('TOO LITTLE WORKSPACE IN CC_FDFBMAT2: ') 404 ENDIF 405 406C ------------------------------------------- 407C Read the CC reference amplitudes from disk: 408C ------------------------------------------- 409 IOPT = 3 410 CALL CC_RDRSP('R0',0,1,IOPT,MODEL,WORK(KT1AMPSAV), 411 * WORK(KT2AMPSAV)) 412 413C -------------------------------------------- 414C Read the T^A reference amplitudes from disk: 415C -------------------------------------------- 416 IOPT = 3 417 ISYMR = ILSTSYM(LISTR,IDLSTR) 418 CALL CC_RDRSP(LISTR,IDLSTR,ISYMR,IOPT,MODEL, 419 & WORK(KT1AMPA),WORK(KT2AMPA)) 420 IF (.NOT.CCS) CALL CCLR_DIASCL(WORK(KT2AMPA),TWO,1) 421 422*---------------------------------------------------------------------* 423* Add delta x t^A to cluster amplitudes and recalculate the 424* response intermediates and the Eta^B vector: 425* Eta^B{t^0 + delta x t^A} 426*---------------------------------------------------------------------* 427 428C ------------------------------------------------------------- 429C add finite displadement to t^0 and recalculate intermediates: 430C ------------------------------------------------------------- 431 CALL DCOPY(NT1AMX,WORK(KT1AMPSAV),1,WORK(KT1AMP0),1) 432 IF (.NOT.CCS) 433 & CALL DCOPY(NT2AMX,WORK(KT2AMPSAV),1,WORK(KT2AMP0),1) 434 CALL DAXPY(NT1AMX,DELTA,WORK(KT1AMPA),1,WORK(KT1AMP0),1) 435 IF (.NOT.CCS) 436 & CALL DAXPY(NT2AMX,DELTA,WORK(KT2AMPA),1,WORK(KT2AMP0),1) 437 438 IOPT = 3 439 CALL CC_WRRSP('R0 ',0,1,IOPT,MODEL,WORK(KT0AMP0), 440 & WORK(KT1AMP0),WORK(KT2AMP0),WORK(KEND1A),LWRK1A) 441 442 RSPIM = .TRUE. 443 CALL CCRHSN(WORK(KOMEGA1),WORK(KOMEGA2),WORK(KT1AMP0), 444 * WORK(KT2AMP0),WORK(KEND1A),LWRK1A,'XXX') 445 446C --------------------------------- 447C calculate the transformed vector: 448C --------------------------------- 449 IORDER = 1 450 CALL CC_XIETA(IXETRAN, NXETRAN, IOPTRES, IORDER, LISTL, 451 & FILXI, IDUM, RDUM, 452 & FILETA, IDUM, RDUM, 453 & .FALSE.,0, WORK(KEND0), LWRK0 ) 454 455C IOPT = 3 456C IVEC = IXETRAN(4,1) 457C ISYETA = ISYHOP 458C CALL CC_RDRSP('X1 ',IVEC,ISYETA,IOPT,MODEL, 459C & WORK(KRHO1),WORK(KRHO2)) 460 461 LEN = NT1AMX + NT2AMX 462 IADRF_ETA = IXETRAN(4,1) 463 LUETA = -1 464 CALL WOPEN2(LUETA,FILETA,64,0) 465 CALL GETWA2(LUETA,FILETA,WORK(KRHO1),IADRF_ETA,LEN) 466 CALL WCLOSE2(LUETA,FILETA,'KEEP') 467 468 RHO1N = DDOT(NT1AMX,WORK(KRHO1),1,WORK(KRHO1),1) 469 IF (.NOT.CCS) RHO2N = DDOT(NT2AMX,WORK(KRHO2),1,WORK(KRHO2),1) 470 471 IF (IPRSAV.GT.10) THEN 472 WRITE (LUPRI,*) 'Norm of RHO1(t^0 + delta x t^C): ',RHO1N 473 WRITE (LUPRI,*) 'Norm of RHO2(t^0 + delta x t^C): ',RHO2N 474 END IF 475 476 ! divide by 2*delta and copy to result vector: 477 CALL DSCAL(NT1AMX,HALF*DELINV,WORK(KRHO1),1) 478 IF (.NOT.CCS) CALL DSCAL(NT2AMX,HALF*DELINV,WORK(KRHO2),1) 479 480 CALL DCOPY(NT1AMX,WORK(KRHO1),1,RHODIF1,1) 481 IF (.NOT.CCS) CALL DCOPY(NT2AMX,WORK(KRHO2),1,RHODIF2,1) 482 483*---------------------------------------------------------------------* 484* Substract delta x t^C to cluster amplitudes and recalculate the 485* response intermediates and the F matrix transformed T^B vector: 486* F{t^0 - delta x t^C} t^B 487*---------------------------------------------------------------------* 488 489C ------------------------------------------------------------- 490C add finite displadement to t^0 and recalculate intermediates: 491C ------------------------------------------------------------- 492 CALL DCOPY(NT1AMX,WORK(KT1AMPSAV),1,WORK(KT1AMP0),1) 493 IF (.NOT.CCS) 494 & CALL DCOPY(NT2AMX,WORK(KT2AMPSAV),1,WORK(KT2AMP0),1) 495 CALL DAXPY(NT1AMX,-DELTA,WORK(KT1AMPA),1,WORK(KT1AMP0),1) 496 IF (.NOT.CCS) 497 & CALL DAXPY(NT2AMX,-DELTA,WORK(KT2AMPA),1,WORK(KT2AMP0),1) 498 499 IOPT = 3 500 CALL CC_WRRSP('R0 ',0,1,IOPT,MODEL,WORK(KT0AMP0), 501 & WORK(KT1AMP0),WORK(KT2AMP0),WORK(KEND1A),LWRK1A) 502 503 RSPIM = .TRUE. 504 CALL CCRHSN(WORK(KOMEGA1),WORK(KOMEGA2), 505 & WORK(KT1AMP0),WORK(KT2AMP0),WORK(KEND1A),LWRK1A,'XXX') 506 507C --------------------------------- 508C calculate the transformed vector: 509C --------------------------------- 510 IORDER = 1 511 CALL CC_XIETA(IXETRAN, NXETRAN, IOPTRES, IORDER, LISTL, 512 & FILXI, IDUM, RDUM, 513 & FILETA, IDUM, RDUM, 514 & .FALSE.,0, WORK(KEND0), LWRK0 ) 515 516C IOPT = 3 517C IVEC = IXETRAN(4,1) 518C ISYETA = ISYHOP 519C CALL CC_RDRSP('X1 ',IVEC,ISYETA,IOPT,MODEL, 520C & WORK(KRHO1),WORK(KRHO2)) 521 522 LEN = NT1AMX + NT2AMX 523 IADRF_ETA = IXETRAN(4,1) 524 CALL WOPEN2(LUETA,FILETA,64,0) 525 CALL GETWA2(LUETA,FILETA,WORK(KRHO1),IADRF_ETA,LEN) 526 CALL WCLOSE2(LUETA,FILETA,'DELETE') 527 528 RHO1N = DDOT(NT1AMX,WORK(KRHO1),1,WORK(KRHO1),1) 529 IF (.NOT.CCS) RHO2N = DDOT(NT2AMX,WORK(KRHO2),1,WORK(KRHO2),1) 530 531 IF (IPRSAV.GT.10) THEN 532 WRITE (LUPRI,*) 'Norm of RHO1(t^0 + delta x t^C): ',RHO1N 533 WRITE (LUPRI,*) 'Norm of RHO2(t^0 + delta x t^C): ',RHO2N 534 END IF 535 536 ! divide by 2*delta and substract from final result: 537 CALL DAXPY(NT1AMX,-HALF*DELINV,WORK(KRHO1),1,RHODIF1,1) 538 IF (.NOT.CCS) 539 & CALL DAXPY(NT2AMX,-HALF*DELINV,WORK(KRHO2),1,RHODIF2,1) 540 541*---------------------------------------------------------------------* 542* fix the scale factor of the diagonal, print some output and 543* restore t^0 amplitudes and response intermediates on file: 544*---------------------------------------------------------------------* 545 546C ----------------------------------------------------------------- 547C scale diagonal with 1/2: (only for right vectors --> A,B,C,D mat) 548C ----------------------------------------------------------------- 549C CALL CCLR_DIASCL(RHODIF2,TWO,1) 550 551 IF (IPRSAV.GT.10) THEN 552 WRITE (LUPRI,*) 'RESULT VECTOR FROM CC_FDFBMAT2:' 553 CALL CC_PRP(RHODIF1,RHODIF2,1,1,1) 554 ENDIF 555 556C -------------------------------------------- 557C Restore the CC reference amplitudes on disk: 558C -------------------------------------------- 559 CALL DCOPY(NT1AMX,WORK(KT1AMPSAV),1,WORK(KT1AMP0),1) 560 IF (.NOT.CCS) 561 & CALL DCOPY(NT2AMX,WORK(KT2AMPSAV),1,WORK(KT2AMP0),1) 562 563 IOPT = 3 564 CALL CC_WRRSP('R0 ',0,1,IOPT,MODEL,WORK(KT0AMP0),WORK(KT1AMP0), 565 & WORK(KT2AMP0),WORK(KEND1A),LWRK1A) 566 567 RSPIM = .TRUE. 568 CALL CCRHSN(WORK(KOMEGA1),WORK(KOMEGA2), 569 & WORK(KT1AMP0),WORK(KT2AMP0),WORK(KEND1A),LWRK1A,'XXX') 570 571 IF (IPRSAV .GT. 5) THEN 572 CALL AROUND(' END OF CC_FDFBMAT2:') 573 ENDIF 574 575 IPRINT = IPRSAV 576 577 RETURN 578 END 579*=====================================================================* 580* END OF SUBROUTINE CC_FDFBMAT2: * 581*=====================================================================* 582 583*=====================================================================* 584 SUBROUTINE CC_FDFBMAT3(LISTR,IDLSTR,RHODIF1,RHODIF2, 585 & LABEL,IRELAX,WORK,LWORK) 586C 587C--------------------------------------------------------------------- 588C Test routine for calculating the Eta{O} vector to test F{O} 589C Sonia, dec 1999 590C--------------------------------------------------------------------- 591C 592#include "implicit.h" 593#include "priunit.h" 594#include "maxorb.h" 595#include "iratdef.h" 596#include "ccorb.h" 597#include "aovec.h" 598#include "ccsdinp.h" 599#include "cclr.h" 600#include "ccsdsym.h" 601#include "ccsdio.h" 602#include "leinf.h" 603#include "ccfro.h" 604#include "ccroper.h" 605#include "cclists.h" 606C 607 DIMENSION WORK(LWORK) 608 DIMENSION RHODIF1(*), RHODIF2(*) 609 CHARACTER*(3) LISTL, LISTR 610 CHARACTER*(8) FILXI, FILETA 611 CHARACTER MODEL*(10) 612 INTEGER IXETRAN(MXDIM_XEVEC,3) 613C LOGICAL LORX 614C 615 INTEGER IR1TAMP, IRHSR1, IETA1, IROPER, ILSTSYM 616C 617 PARAMETER (ONE=1.0d0, ZERO=0.0d0, TWO=2.0d0, HALF=0.5d0) 618 PARAMETER (DELTA = 1.0D-07, DELINV = 1.0D+7) 619C 620 IF (IPRINT.GT.5) THEN 621 CALL AROUND('in CC_FDFBMAT3: calculate XIETA for ITEST 5') 622 ENDIF 623C 624C---------------------------------------------- 625C set up IXETRAN array: 626C---------------------------------------------- 627C 628C IRELAX = 0 629C IF (LORX) THEN 630C IRELAX = IR1TAMP(LABEL,LORX,FREQ,ISYHOP) 631C END IF 632 633 ! set the IXETRAN array for one (XI,ETA) pair 634 IXETRAN(1,1) = IROPER(LABEL,ISYHOP) 635 IXETRAN(2,1) = 0 636 IXETRAN(3,1) = -1 637 IXETRAN(4,1) = 0 638 IXETRAN(5,1) = IRELAX 639C 640 LISTL = 'L0 ' 641 FILXI = 'CCFDFBO1' 642 FILETA = 'CCFDFBX1' 643 IOPTRES = 0 644 NXETRAN = 1 645C 646C---------------------------------------------- 647C initializations & work space allocations. 648C---------------------------------------------- 649C 650 CALL DZERO(RHO1,NT1AMX) 651 CALL DZERO(RHO2,NT2AMX) 652C 653 IPRSAV = IPRINT 654 IPRINT = 0 655C 656 MODEL = 'UNKNOWN' 657 IF (CCS) MODEL = 'CCS' 658 IF (CC2) MODEL = 'CC2' 659 IF (CCSD) MODEL = 'CCSD' 660C 661 KT1AMPSAV = 1 662 KT2AMPSAV = KT1AMPSAV + NT1AMX 663 KT1AMPA = KT2AMPSAV + NT2AMX 664 KT2AMPA = KT1AMPA + NT1AMX 665 KEND0 = KT2AMPA + NT2AMX 666 LWRK0 = LWORK - KEND0 667C 668 KT0AMP0 = KEND0 669 KT1AMP0 = KT0AMP0 + 2*NALLAI(1) 670 KOMEGA1 = KT1AMP0 + NT1AMX 671 KOMEGA2 = KOMEGA1 + NT1AMX 672 KT2AMP0 = KOMEGA2 + MAX(NT2AMX,2*NT2ORT(1),NT2AO(1)) 673 KSCR2 = KT2AMP0 + NT2AMX 674 KEND1A = KSCR2 + NT2AMX + NT1AMX 675 LWRK1A = LWORK - KEND1A 676C 677 KRHO1 = KEND0 678 KRHO2 = KRHO1 + NT1AMX 679 KEND1B = KRHO2 + NT2AMX 680 LWRK1B = LWORK - KEND1B 681C 682 IF (LWRK1A.LT.0 .OR. LWRK1B.LT.0) THEN 683 WRITE(LUPRI,*) 'TOO LITTLE WORK SPACE IN CC_FDFBMAT3::' 684 WRITE(LUPRI,*) 'AVAILABLE: LWORK = ',LWORK 685 WRITE(LUPRI,*) 'NEEDED (AT LEAST) = ',MAX(KEND1A,KEND1B) 686 CALL QUIT('TOO LITTLE WORKSPACE IN CC_FDFBMAT2: ') 687 ENDIF 688 689C ------------------------------------------- 690C Read the CC reference amplitudes from disk: 691C ------------------------------------------- 692 IOPT = 3 693 CALL CC_RDRSP('R0',0,1,IOPT,MODEL,WORK(KT1AMPSAV), 694 * WORK(KT2AMPSAV)) 695 696C -------------------------------------------- 697C Read the T^A reference amplitudes from disk: 698C -------------------------------------------- 699 IOPT = 3 700 ISYMR = ILSTSYM(LISTR,IDLSTR) 701 CALL CC_RDRSP(LISTR,IDLSTR,ISYMR,IOPT,MODEL, 702 & WORK(KT1AMPA),WORK(KT2AMPA)) 703 IF (.NOT.CCS) CALL CCLR_DIASCL(WORK(KT2AMPA),TWO,1) 704 705*---------------------------------------------------------------------* 706* Add delta x t^A to cluster amplitudes and recalculate the 707* response intermediates and the Eta^B vector: 708* Eta^B{t^0 + delta x t^A} 709*---------------------------------------------------------------------* 710C 711C --------------------------------- 712C calculate the transformed vector: 713C --------------------------------- 714 IORDER = 1 715 CALL CC_XIETA(IXETRAN, NXETRAN, IOPTRES, IORDER, LISTL, 716 & FILXI, IDUM, RDUM, 717 & FILETA, IDUM, RDUM, 718 & .FALSE.,0, WORK(KEND0), LWRK0 ) 719 720C IOPT = 3 721C IVEC = IXETRAN(4,1) 722C ISYETA = ISYHOP 723C CALL CC_RDRSP('X1 ',IVEC,ISYETA,IOPT,MODEL, 724C & WORK(KRHO1),WORK(KRHO2)) 725 726 LEN = NT1AMX + NT2AMX 727 IADRF_ETA = IXETRAN(4,1) 728 LUETA = -1 729 CALL WOPEN2(LUETA,FILETA,64,0) 730 CALL GETWA2(LUETA,FILETA,WORK(KRHO1),IADRF_ETA,LEN) 731 CALL WCLOSE2(LUETA,FILETA,'KEEP') 732 733 RHO1N = DDOT(NT1AMX,WORK(KRHO1),1,WORK(KRHO1),1) 734 IF (.NOT.CCS) RHO2N = DDOT(NT2AMX,WORK(KRHO2),1,WORK(KRHO2),1) 735 736 IF (IPRSAV.GT.10) THEN 737 WRITE (LUPRI,*) 'Norm of RHO1(t^0 + delta x t^C): ',RHO1N 738 WRITE (LUPRI,*) 'Norm of RHO2(t^0 + delta x t^C): ',RHO2N 739 END IF 740 741 ! divide by 2*delta and copy to result vector: 742 CALL DSCAL(NT1AMX,HALF*DELINV,WORK(KRHO1),1) 743 IF (.NOT.CCS) CALL DSCAL(NT2AMX,HALF*DELINV,WORK(KRHO2),1) 744 745 CALL DCOPY(NT1AMX,WORK(KRHO1),1,RHODIF1,1) 746 IF (.NOT.CCS) CALL DCOPY(NT2AMX,WORK(KRHO2),1,RHODIF2,1) 747 748*---------------------------------------------------------------------* 749* Substract delta x t^C to cluster amplitudes and recalculate the 750* response intermediates and the F matrix transformed T^B vector: 751* F{t^0 - delta x t^C} t^B 752*---------------------------------------------------------------------* 753 754C ------------------------------------------------------------- 755C add finite displadement to t^0 and recalculate intermediates: 756C ------------------------------------------------------------- 757 CALL DCOPY(NT1AMX,WORK(KT1AMPSAV),1,WORK(KT1AMP0),1) 758 IF (.NOT.CCS) 759 & CALL DCOPY(NT2AMX,WORK(KT2AMPSAV),1,WORK(KT2AMP0),1) 760 CALL DAXPY(NT1AMX,-DELTA,WORK(KT1AMPA),1,WORK(KT1AMP0),1) 761 IF (.NOT.CCS) 762 & CALL DAXPY(NT2AMX,-DELTA,WORK(KT2AMPA),1,WORK(KT2AMP0),1) 763 764 IOPT = 3 765 CALL CC_WRRSP('R0 ',0,1,IOPT,MODEL,WORK(KT0AMP0), 766 & WORK(KT1AMP0),WORK(KT2AMP0),WORK(KEND1A),LWRK1A) 767 768 RSPIM = .TRUE. 769 CALL CCRHSN(WORK(KOMEGA1),WORK(KOMEGA2), 770 & WORK(KT1AMP0),WORK(KT2AMP0),WORK(KEND1A),LWRK1A,'XXX') 771 772C --------------------------------- 773C calculate the transformed vector: 774C --------------------------------- 775 IORDER = 1 776 CALL CC_XIETA(IXETRAN, NXETRAN, IOPTRES, IORDER, LISTL, 777 & FILXI, IDUM, RDUM, 778 & FILETA, IDUM, RDUM, 779 & .FALSE.,0, WORK(KEND0), LWRK0 ) 780 781C IOPT = 3 782C IVEC = IXETRAN(4,1) 783C ISYETA = ISYHOP 784C CALL CC_RDRSP('X1',IVEC,ISYETA,IOPT,MODEL, 785C & WORK(KRHO1),WORK(KRHO2)) 786 787 LEN = NT1AMX + NT2AMX 788 IADRF_ETA = IXETRAN(4,1) 789 CALL WOPEN2(LUETA,FILETA,64,0) 790 CALL GETWA2(LUETA,FILETA,WORK(KRHO1),IADRF_ETA,LEN) 791 CALL WCLOSE2(LUETA,FILETA,'DELETE') 792 793 RHO1N = DDOT(NT1AMX,WORK(KRHO1),1,WORK(KRHO1),1) 794 IF (.NOT.CCS) RHO2N = DDOT(NT2AMX,WORK(KRHO2),1,WORK(KRHO2),1) 795 796 IF (IPRSAV.GT.10) THEN 797 WRITE (LUPRI,*) 'Norm of RHO1(t^0 + delta x t^C): ',RHO1N 798 WRITE (LUPRI,*) 'Norm of RHO2(t^0 + delta x t^C): ',RHO2N 799 END IF 800 801 ! divide by 2*delta and substract from final result: 802 CALL DAXPY(NT1AMX,-HALF*DELINV,WORK(KRHO1),1,RHODIF1,1) 803 IF (.NOT.CCS) 804 & CALL DAXPY(NT2AMX,-HALF*DELINV,WORK(KRHO2),1,RHODIF2,1) 805 806*---------------------------------------------------------------------* 807* fix the scale factor of the diagonal, print some output and 808* restore t^0 amplitudes and response intermediates on file: 809*---------------------------------------------------------------------* 810 811C ----------------------------------------------------------------- 812C scale diagonal with 1/2: (only for right vectors --> A,B,C,D mat) 813C ----------------------------------------------------------------- 814C CALL CCLR_DIASCL(RHODIF2,TWO,1) 815 816 IF (IPRSAV.GT.10) THEN 817 WRITE (LUPRI,*) 'RESULT VECTOR FROM CC_FDFBMAT2:' 818 CALL CC_PRP(RHODIF1,RHODIF2,1,1,1) 819 ENDIF 820 821C -------------------------------------------- 822C Restore the CC reference amplitudes on disk: 823C -------------------------------------------- 824 CALL DCOPY(NT1AMX,WORK(KT1AMPSAV),1,WORK(KT1AMP0),1) 825 IF (.NOT.CCS) 826 & CALL DCOPY(NT2AMX,WORK(KT2AMPSAV),1,WORK(KT2AMP0),1) 827 828 IOPT = 3 829 CALL CC_WRRSP('R0 ',0,1,IOPT,MODEL,WORK(KT0AMP0),WORK(KT1AMP0), 830 & WORK(KT2AMP0),WORK(KEND1A),LWRK1A) 831 832 RSPIM = .TRUE. 833 CALL CCRHSN(WORK(KOMEGA1),WORK(KOMEGA2), 834 & WORK(KT1AMP0),WORK(KT2AMP0),WORK(KEND1A),LWRK1A,'XXX') 835 836 IF (IPRSAV .GT. 5) THEN 837 CALL AROUND(' END OF CC_FDFBMAT2:') 838 ENDIF 839 840 IPRINT = IPRSAV 841 842 RETURN 843 END 844*=====================================================================* 845* END OF SUBROUTINE CC_FDFBMAT3: * 846*=====================================================================* 847 848