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*---------------------------------------------------------------------* 20c/* Deck CC_FBTST */ 21*=====================================================================* 22 SUBROUTINE CC_FBTST(WORK,LWORK) 23*---------------------------------------------------------------------* 24* 25* Purpose: provide some tests for the F{B}T{A} transformation 26* calculated from a (one- and two-electron) derivative operator 27* B refers to the (two-electron) perturbation operator in F{B}, 28* A to the other perturbation in T{A} 29* Sonia Coriani / Christof Haettig, 1998/99 30*---------------------------------------------------------------------* 31#if defined (IMPLICIT_NONE) 32 IMPLICIT NONE 33#else 34# include "implicit.h" 35#endif 36#include "priunit.h" 37#include "ccsdinp.h" 38#include "ccsdsym.h" 39#include "ccorb.h" 40#include "maxorb.h" 41#include "ccroper.h" 42#include "ccr1rsp.h" 43#include "cco1rsp.h" 44#include "ccx1rsp.h" 45#include "ccfro.h" 46 47* local parameters: 48 LOGICAL LOCDBG 49 PARAMETER (LOCDBG = .FALSE.) 50 INTEGER MXFBTRAN 51 PARAMETER (MXFBTRAN = 20) 52 INTEGER LWORK 53 54#if defined (SYS_CRAY) 55 REAL WORK(LWORK) 56 REAL FREQB, FREQA, WRKDLM, DUMMY 57 REAL DDOT 58 REAL ZERO, ONE, TWO, FOUR, FIVE 59 REAL RDUM 60#else 61 DOUBLE PRECISION WORK(LWORK) 62 DOUBLE PRECISION FREQB, FREQA, WRKDLM, DUMMY 63 DOUBLE PRECISION DDOT 64 DOUBLE PRECISION ZERO, ONE, TWO, FOUR, FIVE 65 DOUBLE PRECISION RDUM 66#endif 67 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 68 PARAMETER (FOUR = 4.0D0, FIVE = 5.0D0) 69 70 INTEGER IFBTRAN(5,MXFBTRAN), NFBTRAN 71 INTEGER ISYM0, ISYHOP, ISYOPA, IOPERB 72 INTEGER ITEST, IREAL, IRELAX, IOPT, ILEFT, IRIGHT, IDXR1 73 INTEGER IDETA1, IDRHSR1, IDXR1HF 74 INTEGER IDUM, IDTA1 75 76 CHARACTER*(3) LISTL, LISTR 77 CHARACTER*(8) LABELB, LABELA, FILFBTA 78 CHARACTER*(10) MODEL 79 80 LOGICAL LORXB, LORXA, LRSP, LTWO, FD_ON_FMAT, FD_ON_ETA 81 82 INTEGER KT0, KCMOPQ, KCMOHQ 83 INTEGER KEND0, KEND1A, LWRK1A 84 INTEGER KEND1, KDUM, KEND2, LWRK0, LWRK1, LWRK2, KSCR2 85 INTEGER KKAPPA, KRMAT, KRHO1, KRHO2, KOMEGA1, KOMEGA2 86 INTEGER KRHS1, KRHS2, KT2AMP0, KT1AMP0, IORDER 87 88 INTEGER IOPTRES, N2VEC 89 90* external function: 91 INTEGER ILSTSYM 92 INTEGER IROPER 93 INTEGER IR1TAMP 94 INTEGER IR1KAPPA 95 INTEGER IL1ZETA 96 INTEGER IRHSR1 97 INTEGER IETA1 98 99*---------------------------------------------------------------------* 100* set up information for rhs vectors for the different tests: 101*---------------------------------------------------------------------* 102* number of simultaneous transformations: 103 NFBTRAN = 1 104*---------------------------------------------------------------------* 105* test 1: use zeroth-order Hamiltonian as perturbation. 106* This gives a transformed vector which is 5 * {F T^A} 107* . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108c ITEST = 1 109c LISTL = 'L0' 110c FILFBTA = 'CC_FBMAT' !output file (debug) 111 !The B operator is set uqual to zero-order Hamilton operator 112c LABELB = 'HAM0 ' 113c LORXB = .TRUE. 114c LTWO = .TRUE. 115c FREQB = ZERO 116c ISYHOP = 1 117 !The A operator and T{A} 118c IOPTRES = 1 119c LISTR = 'R1' 120c ISYOPA = 1 121c LABELA = 'ZDIPLEN ' 122c LORXA = .TRUE. 123c FREQA = ZERO 124c 125* . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126* 127* test 2: test the F^B matrix transformation for a non-relaxed 128* one-electron perturbation against the old implementation. 129* (F^A) ---> working 130* (does only require that the integrals for the operator 131* are available on file) 132* . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133c 134c ITEST = 2 135c LISTL = 'L0' 136c IOPTRES = 1 137c FILFBTA = 'CC_FBMAT' 138 !The B operator is a non relaxed one-electron operator 139c LABELB = 'ZDIPLEN ' 140c LORXB = .FALSE. 141c LTWO = .FALSE. 142c FREQB = ZERO 143c ISYHOP = 1 144 !the A operator and T{A} 145c ISYOPA = 1 146c LABELA = 'ZDIPLEN ' 147c LORXA = .FALSE. 148c FREQA = ZERO 149c LISTR = 'R1' 150* . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151* 152* test 3: set differentiated integrals to zero and test only the 153* orbital relaxation & reorthogonalization contributions 154* to the F^A matrix transformed vector. 155* (requires that the CPHF equations for the `reference' 156* operator, specified here, have been solved and the Kappa 157* vector is available on disc) 158* . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 159c 160 ITEST = 3 161 LISTL = 'L0' 162 IOPTRES = 1 163 FILFBTA = 'CC_FBMAT' !output file (debug) 164 !The B operator is DUMMYOP. It will be defined later (vide infra) 165 LABELB = 'ZDIPLEN ' 166 LORXB = .TRUE. 167 LTWO = .FALSE. 168 FREQB = ZERO 169 ISYHOP = 1 170 !the A operator and T{A} 171 ISYOPA = 1 172 LISTR = 'R1' 173 LABELA = 'ZDIPLEN ' 174 LORXA = .TRUE. 175 FREQA = ZERO 176* . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 177* 178* test 4: test the xi and eta vectors for a `orbital-relaxed' 179* one-electron perturbation. 180* (requires that the CPHF equations for the operator specified 181* have been solved and the Kappa vector is available on disc) 182*---------------------------------------------------------------------* 183c ITEST = 4 184c LABEL = 'ZDIPLEN ' 185c LORX = .TRUE. 186c LTWO = .FALSE. 187c FREQ = ZERO 188c ISYHOP = 1 189* . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 190* 191* test 5: test use of London orbitals in FBTA (LAO) 192*---------------------------------------------------------------------* 193c ITEST = 5 194c LISTL = 'L0' 195c IOPTRES = 1 196c FILFBTA = 'CC_FBMAT' !output file (debug) 197c 198c !The B operator is a component of magnetic field (x) 199c LABELB = 'dh/dBX ' 200c LORXB = .TRUE. 201c LTWO = .TRUE. 202c FREQB = ZERO 203c FREQB = 0.09d+00 204c ISYHOP = 1 205 !The A operator and T{A} amplitude 206c ISYOPA = 1 207c LISTR = 'R1' 208c LABELA = 'YDIPLEN ' 209c LORXA = .TRUE. 210c FREQA = ZERO 211c FREQA = 0.09d+00 212*---------------------------------------------------------------------* 213* Special test case ITEST = 5 214* Force calculation of those vectors we need to be able to do FD on eta 215*---------------------------------------------------------------------* 216c ! allow extensions in the vector lists for the next few lines 217c ! to calculate the (X1)^B (eta) vectors ourselves for ITEST 5 218c 219c IF (ITEST.EQ.5) THEN 220c LOPROPN = .TRUE. 221c LO1OPN = .TRUE. 222c LX1OPN = .TRUE. 223c LR1OPN = .TRUE. 224c IRELAX = 0 225c IF (LORXB) THEN 226c IRELAX = IR1KAPPA(LABELB,FREQB,ISYHOP) 227c END IF 228c 229c LPDBSOP(IROPER(LABELB,ISYHOP)) = LTWO 230c 231c ! set the IXETRAN array for one (XI,ETA) pair 232c IXETRAN(1,1) = IROPER(LABELB,ISYHOP) 233c IXETRAN(2,1) = 0 234c IXETRAN(3,1) = IRHSR1(LABELB,LORXB,FREQB,ISYHOP) 235c IXETRAN(4,1) = IETA1(LABELB,LORXB,FREQB,ISYHOP) 236c IXETRAN(5,1) = IRELAX 237c ! disallow again extension in the vector lists 238c LOPROPN = .FALSE. 239c LO1OPN = .FALSE. 240c LX1OPN = .FALSE. 241c 242c 243*---------------------------------------------------------------------* 244* call CC_XIETA to calculate the Xi and Eta vectors: 245*---------------------------------------------------------------------* 246c LISTL = 'L0 ' 247c FILXI = 'O1 ' 248c FILETA = 'X1 ' 249c IOPTRES = 3 250c 251c KEND0 = 1 252c KEND1 = KEND0 + 1 253c LWRK1 = LWORK - KEND1 254c IF (LWRK1 .LT. 0) THEN 255c CALL QUIT('Insufficient work space in CC_XETST. (1)') 256c END IF 257c 258c IDUM = 0 259c WORK(KEND1-1) = WORK(IDUM) 260c WRITE (LUPRI,*) 'WORK(0) before entering CC_XIETA:',WORK(KEND1-1) 261c 262c CALL CC_XIETA(IXETRAN, NXETRAN, IOPTRES, IORDER, LISTL, 263c & FILXI, IXDOTS, XCONS, 264c & FILETA, IEDOTS, ECONS, 265c & .FALSE.,0, WORK(KEND1), LWRK1 ) 266c 267c WRITE (LUPRI,*) 'returned from CC_XIETA... WORK(0) :',WORK(KEND1) 268c CALL FLSHFO(LUPRI) 269c 270c END IF 271*--------------------------------------------------------------------* 272* End of special test case 273*--------------------------------------------------------------------* 274 275c get the index of the kappa vector (kappa{B}) which is saved on the 276c same file on which the corresponding T amplitude responses (T{B}) 277c would be saved a value of 0 indicates `unrelaxed' 278 279 IRELAX = 0 280 IF (LORXB) THEN 281 IRELAX = IR1KAPPA(LABELB,FREQB,ISYHOP) !it could be frequency dep 282 END IF 283 284c set logical for perturbation-dependent basis sets (for B); this will 285c overwrite the default already saved on the common block 286 287 LPDBSOP(IROPER(LABELB,ISYHOP)) = LTWO 288 289c set the IFBTRAN array for one F{B} T{A} result vector 290 291 IFBTRAN(1,1) = IROPER(LABELB,ISYHOP) ! index B operator 292 IFBTRAN(2,1) = 0 ! index left vector (L0) 293 IFBTRAN(3,1) = IR1TAMP(LABELA,LORXA,FREQA,ISYOPA) ! index right vector T{A} (R1) 294 IFBTRAN(4,1) = 0 ! index result vector (return in memory) 295 IFBTRAN(5,1) = IRELAX ! index kappa vector 296 297c choose finite difference path: 298 FD_ON_FMAT = .FALSE. ! f.d. on usual F matrix 299 FD_ON_ETA = .TRUE. ! f.d. on Eta vector 300 301 IF (ITEST.EQ.3) THEN 302 303 ! for test 3, replace now the operator labels on the 304 ! common blocks by 'DUMMYOP ', which should be interpreted 305 ! inside of CC_XIETA/CC_FBTA as a zero one-electr. operator, 306 ! but associated with a orb.-relax. (kappa) vector. 307 308 ! index for a R1 type vect for LABELB operator (B) 309 IDXR1 = IR1TAMP(LABELB,LORXB,FREQB,ISYHOP) 310 ! index for correponding CPHF vector 311 IDXR1HF = IR1KAPPA(LABELB,FREQB,ISYHOP) 312 KT0 = 1 313 KEND1 = KT0 + 2*NALLAI(ISYHOP) 314 LWRK1 = LWORK - KEND1 315 CALL CC_RDHFRSP('R1 ',IDXR1HF,ISYHOP,WORK(KT0)) 316 317 IDRHSR1 = IRHSR1(LABELB,LORXB,FREQB,ISYHOP) !index for O1, RHS of R1 equations (TB), dvs csi^B (O1) 318 IDETA1 = IETA1(LABELB,LORXB,FREQB,ISYHOP) !index for X1, the Eta^B vector (X1) 319 LABELB = 'DUMMYOP ' 320 LBLOPR(IFBTRAN(1,1)) = LABELB ! operator common block 321 LRTLBL(IDXR1) = LABELB ! common block for R1{B}/kappa{B} 322 LBLO1(IDRHSR1) = LABELB !operator label of IDRHSR1 (O1) 323 LBLX1(IDETA1) = LABELB !operator label of IDETA1 (X1) 324* 325* substitute DUMMY to ZDIPLEN on file for Kappa (use IOPT = 4) 326* 327 IOPT = 4 328 CALL CC_WRRSP('R1 ',IDXR1HF,ISYHOP,IOPT,MODEL,WORK(KT0), 329 & DUMMY,DUMMY,WORK(KEND1),LWRK1) 330 331 WRITE (LUPRI,*) 'CC_FBTST: orbital relaxation vector:' 332 CALL OUTPUT(WORK(KT0),1,2*NALLAI(ISYHOP),1,1, 333 & 2*NALLAI(ISYHOP),1,1,LUPRI) 334 335 END IF 336*---------------------------------------------------------------------* 337* Print general infos 338*---------------------------------------------------------------------* 339 WRITE (LUPRI,*) 'Test case nr ITEST =',ITEST 340 WRITE (LUPRI,*) 'B operator LABELB =',LABELB 341 WRITE (LUPRI,*) '2-el operator? LTWO =',LTWO 342 WRITE (LUPRI,*) 'relaxed oper? LORXB =',LORXB 343 WRITE (LUPRI,*) 'of frequency? FREQB =',FREQB 344 WRITE (LUPRI,*) 'IROPER =',IFBTRAN(1,1) 345 WRITE (LUPRI,*) 'ILEFT =',IFBTRAN(2,1) 346 WRITE (LUPRI,*) 'IRIGHT =',IFBTRAN(3,1) 347 WRITE (LUPRI,*) 348 & 'IRELAX =',IFBTRAN(5,1),' LRTLBL(RELAX):', LRTLBL(IRELAX) 349 WRITE (LUPRI,*) 'A op. in T{A} LABELA =',LABELA, 350 & ' LRTLBL(R1):', LRTLBL(IFBTRAN(3,1)) 351 WRITE (LUPRI,*) 'relaxed? LORXA =',LORXA 352 WRITE (LUPRI,*) 'of frequency? FREQA =',FREQA 353 354 N2VEC = 1 355 IF (CCS) N2VEC = 0 356 357*---------------------------------------------------------------------* 358* call CC_FBTA to calculate F{B}T{A} transformed vector: 359*---------------------------------------------------------------------* 360 361 KEND0 = 1 362 LWRK0 = LWORK 363 CALL CC_FBTA(IFBTRAN, NFBTRAN, IOPTRES, LISTL, LISTR, 364 & FILFBTA, IDUM, RDUM,0, 365 & WORK(KEND0), LWRK0) 366 367c CALL CC_FBTA(IFBTRAN, NFBTRAN, IOPTRES, LISTL, LISTR, 368c & FILFBTA, IFBDOTS, FBCON, MXVEC, 369c & WORK(KEND0), LWRK0) 370 371*---------------------------------------------------------------------* 372* calculate the reference vector and compare 373* with the results from the CC_FBTA routine: 374*---------------------------------------------------------------------* 375 376* ------------------------------------------------------------------- * 377* ITEST = 1 Test against F matrix 378* ------------------------------------------------------------------- * 379 IF (ITEST.EQ.1) THEN 380 ILEFT = IFBTRAN(2,1) 381 IRIGHT = IFBTRAN(3,1) 382 CALL AROUND('ITEST =1, 5 times the Fmatrix transformation:') 383 CALL CCFTRANSC(LISTL, ILEFT, LISTR, IRIGHT, WORK(KEND0), 384c CALL CC_FTRAN(LISTL, ILEFT, LISTR, IRIGHT, WORK(KEND0), 385 & LWRK0) 386 END IF 387 388 IF (ITEST.EQ.1) RETURN 389 390* -------------------------------------------------------------------- * 391* ITEST > 1 Test against others 392* -------------------------------------------------------------------- * 393 KKAPPA = 1 394 KRMAT = KKAPPA + 2*NALLAI(ISYHOP) 395 KCMOPQ = KRMAT + N2BST(ISYHOP) 396 KCMOHQ = KCMOPQ + NGLMDT(ISYHOP) 397 KT1AMP0 = KCMOHQ + NGLMDT(ISYHOP) 398 KOMEGA1 = KT1AMP0 + NT1AMX 399 KOMEGA2 = KOMEGA1 + NT1AM(ISYHOP) 400 KRHS1 = KOMEGA2 + NT2AM(ISYHOP) 401 KRHS2 = KRHS1 + NT1AM(ISYHOP) 402 KEND1 = KRHS2 + NT2AM(ISYHOP) 403 LWRK1 = LWORK - KEND1 404 405 KT2AMP0 = KOMEGA2 + MAX(NT2AMX,2*NT2ORT(1),NT2AO(1)) 406 KSCR2 = KT2AMP0 + NT2AMX 407 KEND1A = KSCR2 + NT2AMX + NT1AMX 408 LWRK1A = LWORK - KEND1A 409 410 IF (LWRK1.LT.0 .OR. LWRK1A.LT.0) THEN 411 CALL QUIT('Insufficient memory in CC_FBTST.') 412 END IF 413 414c 415 IF (LORXB) THEN 416 IF (LABELB.EQ.'HAM0 ') THEN 417 CALL DZERO(WORK(KKAPPA),2*NALLAI(ISYHOP)) 418 ELSE 419 CALL FLSHFO(LUPRI) 420 CALL CC_RDHFRSP('R1 ',IRELAX,ISYHOP,WORK(KKAPPA)) 421 CALL FLSHFO(LUPRI) 422 END IF 423 ELSE 424 CALL DZERO(WORK(KKAPPA),2*NALLAI(ISYHOP)) 425 END IF 426 427* ------------------------------ 428* get orbital connection matrix: 429* ------------------------------ 430 IF (LTWO) THEN 431 IOPERB = IROPER(LABELB,ISYHOP) 432 IORDER= 1 433 CALL CC_GET_RMAT(WORK(KRMAT),IOPERB,IORDER, 434 & ISYHOP,WORK(KEND1),LWRK1) 435 END IF 436 437* ------------------------------------ 438* construct the derivative CMO vector: 439* ------------------------------------ 440 IF (LORXB.OR.LTWO) THEN 441 IREAL = ISYMAT(IROPER(LABELB,ISYHOP)) 442 IOPT = 0 443 CALL CC_LAMBDAQ(DUMMY,DUMMY,WORK(KCMOPQ),WORK(KCMOHQ),ISYHOP, 444 & DUMMY,WORK(KKAPPA),WORK(KRMAT),IREAL,IOPT, 445 & WORK(KEND1),LWRK1) 446 END IF 447 448 IF (ITEST.EQ.2) THEN 449* ------------------------------------------------------------ 450* call the usual FA matrix routine to calculate the reference 451* result for the FA transformed vector: 452* ------------------------------------------------------------ 453c ILEFT = IFBTRAN(2,1) 454c IRIGHT = IFBTRAN(3,1) 455c CALL CCLR_FA(LABELB, ISYHOP, ! inp: label/symmetry A 456c & LISTR, IRIGHT, ! inp: B resp. amplit. 457c & LISTL, ILEFT, ! inp: C resp. zeta vec. 458c & WORK(KEND1), LWRK1 )! work space 459c 460 END IF 461 462 IF (LORXB.OR.LTWO) THEN 463* 464* ------------------------------------------------------------ 465* evaluate orbital relaxation and reorthogonalization contrib. 466* to the F matrix transformed vector by 467* a) finite difference on the usual F matrix (mht C) 468* b) finite difference on the Eta{A} vector (CC_XIETA) (mht T) 469* 470* (does only work for totally symmetric perturbations... 471* actually only tested without symmetry...) 472* ------------------------------------------------------------ 473 IF (ISYHOP.NE.1) THEN 474 WRITE (LUPRI,*) 'finite difference test of CC_FBTA not ' 475 WRITE (LUPRI,*) 476 & 'available for non-totally sym. perturbations.' 477 CALL QUIT( 478 & 'CC_FDFBMAT only implemented for tot. sym. perturb.') 479 END IF 480 481 IF (FD_ON_FMAT) THEN 482 ILEFT = IFBTRAN(2,1) 483 IRIGHT = IFBTRAN(3,1) 484 485 IF (IREAL.EQ.-1) THEN 486 WRITE (LUPRI,*) 487 & 'CC_FDFBMAT not applicable for IREAL = -1.' 488 END IF 489 WRITE(LUPRI,*) 'NOW CALL CC_FDFBMAT' 490 CALL CC_FDFBMAT(WORK(KCMOHQ),LISTL,ILEFT,LISTR,IRIGHT, 491 & WORK(KOMEGA1),WORK(KEND1),LWRK1) 492 END IF 493 494 IF (FD_ON_ETA) THEN 495 IRIGHT = IFBTRAN(3,1) 496 CALL CC_FDFBMAT2(LISTR,IRIGHT,WORK(KOMEGA1),WORK(KOMEGA2), 497 & LABELB,IRELAX,WORK(KEND1),LWRK1) 498 CALL AROUND('fin. diff. orbital relaxation F^B T^A vector:') 499 CALL CC_PRP(WORK(KOMEGA1),WORK(KOMEGA2),ISYHOP,1,N2VEC) 500 END IF 501 502 CALL CC_PRP(WORK(KOMEGA1),WORK(KOMEGA2),ISYHOP,1,N2VEC) 503 504 ELSE 505 WRITE (LUPRI,*) 'no orb. relax. or reorth. contribution '// 506 & 'to F^B T^A.' 507 CALL DZERO(WORK(KOMEGA1),NT1AM(ISYHOP)) 508 CALL DZERO(WORK(KOMEGA2),NT2AM(ISYHOP)) 509 END IF 510 511C ------------------------------------------------------ 512C calculate the contribution from one-electron h^1: 513C (not needed, if we do finite difference on Eta vector) 514C ------------------------------------------------------ 515 IF (FD_ON_FMAT .AND. 516 & (LABELB(1:8).NE.'DUMMYOP '.AND. LABELB(1:8).NE.'HAM0 ')) THEN 517 518 CALL QUIT('one-electron contr. not yet available in CC_FBTA.') 519 520 CALL AROUND('one-electron h^1 contribution to Xi vector:') 521 CALL CC_PRP(WORK(KRHS1),WORK(KRHS2),ISYHOP,1,N2VEC) 522 CALL DAXPY(NT1AM(ISYHOP),ONE,WORK(KRHS1),1,WORK(KOMEGA1),1) 523 CALL DAXPY(NT2AM(ISYHOP),ONE,WORK(KRHS2),1,WORK(KOMEGA2),1) 524 525 CALL AROUND('num. orb.-relax. + analyt. one-electr. F^A T^B:') 526 CALL CC_PRP(WORK(KOMEGA1),WORK(KOMEGA2),ISYHOP,1,N2VEC) 527 528 END IF 529 530C ------------------------------------------------- 531C read the result vector from CC_FBTA and compare: 532C ------------------------------------------------- 533c DO ITRAN = 1, NFBTRAN 534c IOPT = 3 535c IVEC = IFBTRAN(3,ITRAN) 536c PRINT *,'CC_FBTST: IVEC = ',IVEC 537c CALL CC_RDRSP('O1 ',IVEC,ISYHOP,IOPT,MODEL, 538c & WORK(KRHS1),WORK(KRHS2)) 539c 540c CALL AROUND('analytical F^B T^A vector:') 541c CALL CC_PRP(WORK(KRHS1),WORK(KRHS2),ISYHOP,1,N2VEC) 542c 543c CALL DAXPY(NT1AM(ISYHOP),-1.0d0,WORK(KOMEGA1),1,WORK(KRHS1),1) 544c CALL DAXPY(NT2AM(ISYHOP),-1.0d0,WORK(KOMEGA2),1,WORK(KRHS2),1) 545c 546c WRITE(LUPRI,*) 'Norm of difference between analytical F^B T^A ' 547c > // 'vector and the numerical result:' 548c WRITE (LUPRI,*) 'singles excitation part:', 549c > DSQRT(DDOT(NT1AM(ISYHOP),WORK(KRHS1),1,WORK(KRHS1),1)) 550c IF (.NOT.CCS) WRITE (LUPRI,*) 'double excitation part: ', 551c > DSQRT(DDOT(NT2AM(ISYHOP),WORK(KRHS2),1,WORK(KRHS2),1)) 552c 553c WRITE (LUPRI,*) 'difference vector:' 554c CALL CC_PRP(WORK(KRHS1),WORK(KRHS2),ISYHOP,1,N2VEC) 555c END DO 556 557 RETURN 558 END 559**=====================================================================* 560