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_XETST */ 21*=====================================================================* 22 SUBROUTINE CC_XETST(WORK,LWORK) 23*---------------------------------------------------------------------* 24* 25* Purpose : provide some tests for the CC_XIETA routine 26* for more detailed information of the tests see below 27* 28* noddy implementation for programmers use only... 29* ... to switch between the different test or for changing 30* the parameters the program must be recompiled... 31* 32* Christof Haettig, 1998/99 33*---------------------------------------------------------------------* 34#if defined (IMPLICIT_NONE) 35 IMPLICIT NONE 36#else 37# include "implicit.h" 38#endif 39#include "priunit.h" 40#include "ccsdinp.h" 41#include "ccsdsym.h" 42#include "ccorb.h" 43#include "maxorb.h" 44#include "ccroper.h" 45#include "ccr1rsp.h" 46#include "cco1rsp.h" 47#include "ccx1rsp.h" 48#include "ccfro.h" 49#include "cclists.h" 50#include "dummy.h" 51 52* local parameters: 53 CHARACTER MSGDBG*(18) 54 PARAMETER (MSGDBG='[debug] CC_XETST> ') 55 56 LOGICAL LOCDBG 57 PARAMETER (LOCDBG = .FALSE.) 58 INTEGER MXXETRAN, MXVEC 59 PARAMETER (MXXETRAN = 20, MXVEC = 10) 60 INTEGER LWORK 61#if defined (SYS_CRAY) 62 REAL WORK(LWORK) 63 REAL FREQ, WRKDLM 64 REAL XCONS(MXVEC,MXXETRAN), ECONS(MXVEC,MXXETRAN) 65 REAL DDOT, PROPRSP, PROPAVE, FF 66 REAL ZERO, ONE, TWO, FOUR, FIVE, FACUP 67#else 68 DOUBLE PRECISION WORK(LWORK) 69 DOUBLE PRECISION FREQ, FREQ1, FREQ2, WRKDLM 70 DOUBLE PRECISION XCONS(MXVEC,MXXETRAN), ECONS(MXVEC,MXXETRAN) 71 DOUBLE PRECISION DDOT, PROPRSP, PROPAVE, FF 72 DOUBLE PRECISION ZERO, ONE, TWO, FOUR, FIVE, FACUP 73#endif 74 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 75 PARAMETER (FOUR = 4.0D0, FIVE = 5.0D0) 76 77 CHARACTER*(3) FILXI, FILETA, LISTL, APROXR12 78 CHARACTER*(8) LABEL, LABEL1, LABEL2 79 CHARACTER*(10) MODEL, CROUT 80 LOGICAL LORX, LORX1, LORX2, LTWO, LTWO1, LTWO2 81 INTEGER IOPTRES, N2VEC 82 INTEGER IXETRAN(MXDIM_XEVEC,MXXETRAN), NXETRAN, IDXR1HF 83 INTEGER IXDOTS(MXVEC,MXXETRAN), IEDOTS(MXVEC,MXXETRAN) 84 INTEGER IADRLQ(MXXETRAN), IADRDQ(MXXETRAN), ISYAMP, IOPTCC2 85 INTEGER ISYM0, ISYHOP, IOPT, IOPER, ITRAN, IVEC, IDUM, ISYETA 86 INTEGER KEND0, ITEST, KEND1A, LWRK1A, IORDER, ILSTETA, ISIGN 87 INTEGER KIT2DEL0, KIT2DELQ, KEND1, KDUM, KEND2, LWRK1, LWRK2 88 INTEGER KXKSI1, KXKSI2, KRHS1, KRHS2, KT1AMP0, KT2AMP0,KT0AMP0 89 INTEGER NDENSQ, NLAMDQ, IPRS, KOMEGA1, KOMEGA2, KSCR2, KT0 90 INTEGER KCMO,KCMO0,KCMOPQ,KCMOHQ,IDLSTL,ISYCTR,IDXR1,INUM,IREAL 91 INTEGER KKAPPA, KRMAT, KLEFT1, KLEFT2, KETA1, KETA2, KRHO1, KRHO2 92 INTEGER IRELAX, IRELAX1, IRELAX2, ISYM1, ISYM2, KCTR1, KCTR2 93 INTEGER KETA12 94 95 96* external function: 97 INTEGER ILSTSYM 98 INTEGER IROPER 99 INTEGER IR1TAMP 100 INTEGER IR1KAPPA 101 INTEGER IL1ZETA 102 INTEGER IRHSR1 103 INTEGER IETA1 104 105*---------------------------------------------------------------------* 106* set up information for rhs vectors for the different tests: 107*---------------------------------------------------------------------* 108* number of simultaneous transformations: 109 NXETRAN = 1 110*---------------------------------------------------------------------* 111* some initialization 112 LORX1 = .FALSE. 113 LORX2 = .FALSE. 114*---------------------------------------------------------------------* 115* test 1: use zeroth-order Hamiltonian as perturbation. this gives 116* a xi vector which is 5 x vector function (omega) and a 117* a eta vector which is 5 x (eta(0) + zeta(0) A). if the 118* cluster amplitude and lagrangian multiplier equations are 119* converged, both vectors should be zero. 120* (does not require that anything is done before CC_XETST) 121* . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 ITEST = 1 123 IORDER = 1 124 LABEL = 'HAM0 ' 125 LORX = .TRUE. 126 LTWO = .TRUE. 127 FREQ = ZERO 128 ISYHOP = 1 129* . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130* 131* test 2: test the xi and eta vector for a non-relaxed one-elctron 132* perturbation against the old implemenation. 133* (does only require that the integrals for the operator 134* are available on file) 135* . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136c ITEST = 2 137c IORDER = 1 138c LABEL = 'ZDIPLEN ' 139c LORX = .FALSE. 140c LTWO = .FALSE. 141c FREQ = ZERO 142c ISYHOP = 1 143* . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144* 145* test 3: set differentiated integrals to zero and test only the 146* orbital relaxation & reorthogonalization part of the xi 147* and eta vectors. 148* (requires that the CPHF equations for the `reference' 149* operator, specified here, have been solved and the Kappa 150* vector is available on disc) 151* . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152c ITEST = 3 153c IORDER = 1 154c LABEL = 'ZDIPLEN ' 155c LORX = .TRUE. 156c LTWO = .FALSE. 157c FREQ = ZERO 158c ISYHOP = 1 159* . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160* 161* test 4: test the xi and eta vectors for a `orbital-relaxed' 162* one-electron perturbation. 163* (requires that the CPHF equations for the operator specified 164* have been solved and the Kappa vector is available on disc) 165* . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 166c ITEST = 4 167c IORDER = 1 168c LABEL = 'ZDIPLEN ' 169c LORX = .TRUE. 170c LTWO = .FALSE. 171c FREQ = ZERO 172c ISYHOP = 1 173* . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 174* 175* test 5: test the xi and eta vectors for a second-order perturbation 176* operator made up from a relaxed field with pert.-dep. basis 177* and an unrelaxed one-electron perturbation 178* (requires that the CPHF equations for the operator specified 179* have been solved and the Kappa vector is available on disc) 180* . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 181c ITEST = 5 182c IORDER = 2 183c LABEL1 = '1DHAM003' 184c LORX1 = .TRUE. 185c LTWO1 = .TRUE. 186c FREQ1 = ZERO 187c LABEL2 = 'ZDIPLEN ' 188c LORX2 = .FALSE. 189c LTWO2 = .FALSE. 190c FREQ2 = ZERO 191c CALL CC_FIND_SO_OP(LABEL1,LABEL2,LABEL,ISYHOP,ISIGN, 192c & INUM,WORK,LWORK) 193c LORX = ( LORX1 .AND. LORX2 ) 194c LTWO = ( LTWO1 .AND. LTWO2 ) 195c FREQ = FREQ1 + FREQ2 196*---------------------------------------------------------------------* 197 198 ! allow extensions in the vector lists for the next few lines 199 LOPROPN = .TRUE. 200 LO1OPN = .TRUE. 201 LX1OPN = .TRUE. 202 203 IRELAX = 0 204 IF (LORX) THEN 205 IRELAX = IR1TAMP(LABEL,LORX,FREQ,ISYHOP) 206 END IF 207 208 209 LPDBSOP(IROPER(LABEL,ISYHOP)) = LTWO 210 211 ! set the IXETRAN array for one (XI,ETA) pair 212 IXETRAN(1,1) = IROPER(LABEL,ISYHOP) 213 IXETRAN(2,1) = 0 214 IXETRAN(3,1) = IRHSR1(LABEL,LORX,FREQ,ISYHOP) 215 IXETRAN(4,1) = IETA1(LABEL,LORX,FREQ,ISYHOP) 216 IXETRAN(5,1) = IRELAX 217 218 IF (IORDER.EQ.2) THEN 219 IRELAX1 = 0 220 IRELAX2 = 0 221 IF (LORX1) IRELAX1 = IR1TAMP(LABEL1,LORX1,FREQ1,ISYM1) 222 IF (LORX2) IRELAX2 = IR1TAMP(LABEL2,LORX2,FREQ2,ISYM2) 223 IXETRAN(5,1) = IRELAX1 224 IXETRAN(6,1) = IRELAX2 225 END IF 226 227 ! disallow again extension in the vector lists 228 LOPROPN = .FALSE. 229 LO1OPN = .FALSE. 230 LX1OPN = .FALSE. 231 232 ! for test 3, replace now the operator labels on the common blocks 233 ! 'DUMMYOP ', interpreted inside of CC_XIETA as a zero one-electr. 234 ! operator, but associated with a orb.-relax. (kappa) vector. 235 IF (ITEST.EQ.3) THEN 236C 237 IDXR1 = IR1TAMP(LABEL,LORX,FREQ,ISYHOP) 238 IDXR1HF = IR1KAPPA(LABEL,FREQ,ISYHOP) 239 KT0 = 1 240 KEND1 = KT0 + 2*NALLAI(ISYHOP) 241 LWRK1 = LWORK - KEND1 242 CALL CC_RDHFRSP('R1 ',IDXR1HF,ISYHOP,WORK(KT0)) 243C 244 LABEL = 'DUMMYOP ' 245 LBLOPR(IXETRAN(1,1)) = LABEL 246 LBLO1(IXETRAN(3,1)) = LABEL 247 LBLX1(IXETRAN(4,1)) = LABEL 248 LRTLBL(IDXR1) = LABEL 249C 250 IOPT = 4 251 CALL CC_WRRSP('R1 ',IDXR1,ISYHOP,IOPT,MODEL,WORK(KT0), 252 & DUMMY,DUMMY,WORK(KEND1),LWRK1) 253C 254 END IF 255C 256 CALL AROUND('CC_XETST: test of CC_XIETA subroutine') 257 WRITE (LUPRI,*) 'ITEST =',ITEST 258 WRITE (LUPRI,*) 'LABEL =',LABEL 259 WRITE (LUPRI,*) 'LTWO =',LTWO 260 WRITE (LUPRI,*) 'LORX =',LORX 261 WRITE (LUPRI,*) 'FREQ =',FREQ 262 WRITE (LUPRI,*) 'IROPER =',IXETRAN(1,1) 263 WRITE (LUPRI,*) 'ILEFT =',IXETRAN(2,1) 264 WRITE (LUPRI,*) 'IRHSR1 =',IXETRAN(3,1) 265 WRITE (LUPRI,*) 'IETA1 =',IXETRAN(4,1) 266 WRITE (LUPRI,*) 'IRELAX =',IXETRAN(5,1) 267C 268 DO I = 2, NXETRAN 269 IXETRAN(1,I) = IXETRAN(1,1) 270 IXETRAN(2,I) = IXETRAN(2,1) 271 IXETRAN(3,I) = IXETRAN(3,1) + I-1 272 IXETRAN(4,I) = IXETRAN(4,1) + I-1 273 IXETRAN(5,I) = IXETRAN(5,1) 274 IXETRAN(6,I) = IXETRAN(6,1) 275 LBLO1(IXETRAN(3,1)+I-1) = LABEL 276 ISYO1(IXETRAN(3,1)+I-1) = ISYHOP 277 LBLX1(IXETRAN(4,1)+I-1) = LABEL 278 ISYX1(IXETRAN(4,1)+I-1) = ISYHOP 279 END DO 280 NO1LBL = MAX(NO1LBL,IXETRAN(3,1)+NXETRAN-1) 281 NX1LBL = MAX(NX1LBL,IXETRAN(4,1)+NXETRAN-1) 282C 283 N2VEC = 1 284 IF (CCS) N2VEC = 0 285 286*---------------------------------------------------------------------* 287* call CC_XIETA to calculate the Xi and Eta vectors: 288*---------------------------------------------------------------------* 289 LISTL = 'L0 ' 290 FILXI = 'O1 ' 291 FILETA = 'X1 ' 292 IOPTRES = 3 293 294 KEND0 = 1 295 KEND1 = KEND0 + 1 296 LWRK1 = LWORK - KEND1 297 IF (LWRK1 .LT. 0) THEN 298 CALL QUIT('Insufficient work space in CC_XETST. (1)') 299 END IF 300 301 IDUM = 0 302 WORK(KEND1-1) = WORK(IDUM) 303 WRITE (LUPRI,*) 'WORK(0) before entering CC_XIETA:',WORK(KEND1-1) 304 305 CALL CC_XIETA(IXETRAN, NXETRAN, IOPTRES, IORDER, LISTL, 306 & FILXI, IXDOTS, XCONS, 307 & FILETA, IEDOTS, ECONS, 308 & .FALSE.,0, WORK(KEND1), LWRK1 ) 309 310 WRITE (LUPRI,*) 'returned from CC_XIETA... WORK(0) :',WORK(KEND1) 311 CALL FLSHFO(LUPRI) 312 313*---------------------------------------------------------------------* 314* calculate the reference values xi vector and compare 315* with the results from the CC_XIETA routine: 316*---------------------------------------------------------------------* 317 KKAPPA = 1 318 KRMAT = KKAPPA + 2*NALLAI(ISYHOP) 319 KCMOPQ = KRMAT + N2BST(ISYHOP) 320 KCMOHQ = KCMOPQ + NGLMDT(ISYHOP) 321 KT1AMP0 = KCMOHQ + NGLMDT(ISYHOP) 322 KOMEGA1 = KT1AMP0 + NT1AMX 323 KOMEGA2 = KOMEGA1 + NT1AM(ISYHOP) 324 KRHS1 = KOMEGA2 + NT2AM(ISYHOP) 325 KRHS2 = KRHS1 + NT1AM(ISYHOP) 326 KEND1 = KRHS2 + NT2AM(ISYHOP) 327 LWRK1 = LWORK - KEND1 328 329 KT2AMP0 = KOMEGA2 + MAX(NT2AMX,2*NT2ORT(1),NT2AO(1)) 330 KSCR2 = KT2AMP0 + NT2AMX 331 KEND1A = KSCR2 + NT2AMX + NT1AMX 332 LWRK1A = LWORK - KEND1A 333 334 KCTR1 = KEND1 335 KCTR2 = KCTR1 + NT1AMX 336 KEND2 = KCTR2 + NT2AMX 337 LWRK2 = LWORK - KEND2 338 339 IF (LWRK1.LT.0 .OR. LWRK1A.LT.0 .OR. LWRK2.LT.0) THEN 340 CALL QUIT('Insufficient memory in CC_XETST.') 341 END IF 342 343C ------------------------------ 344C get orbital relaxation vector: 345C ------------------------------ 346 IF (LORX) THEN 347 IF (LABEL.EQ.'HAM0 ') THEN 348 CALL DZERO(WORK(KKAPPA),2*NALLAI(ISYHOP)) 349 ELSE 350 IDXR1HF = IR1KAPPA(LABEL,FREQ,ISYHOP) 351 CALL CC_RDHFRSP('R1 ',IDXR1HF,ISYHOP,WORK(KKAPPA)) 352 END IF 353 ELSE 354 CALL DZERO(WORK(KKAPPA),2*NALLAI(ISYHOP)) 355 END IF 356 357C ------------------------------ 358C get orbital connection matrix: 359C ------------------------------ 360 IF (LORX.OR.LTWO) THEN 361 IOPER = IROPER(LABEL,ISYHOP) 362 CALL CC_GET_RMAT(WORK(KRMAT),IOPER,1,ISYHOP,WORK(KEND1),LWRK1) 363 END IF 364 365C ------------------------------------ 366C construct the derivative CMO vector: 367C ------------------------------------ 368 IF (LORX.OR.LTWO) THEN 369 IOPT = 0 370 IREAL = ISYMAT(IOPER) 371 CALL CC_LAMBDAQ(DUMMY,DUMMY,WORK(KCMOPQ),WORK(KCMOHQ),ISYHOP, 372 & DUMMY,WORK(KKAPPA),WORK(KRMAT),IREAL,IOPT, 373 & WORK(KEND1),LWRK1) 374 WRITE (LUPRI,*) 'CC_XIETA> CMOPQ vector:' 375 CALL CC_PRONELAO(WORK(KCMOPQ),ISYHOP) 376 WRITE (LUPRI,*) 'CC_XIETA> CMOHQ vector:' 377 CALL CC_PRONELAO(WORK(KCMOHQ),ISYHOP) 378 END IF 379 380 IF (LORX.OR.LTWO) THEN 381 IF (ITEST.EQ.1) THEN 382C ------------------------------------------------------------ 383C call the vector function routine to calculate the reference 384C result for the xi vector: 385C ------------------------------------------------------------ 386 IOPT = 3 387 CALL CC_RDRSP('R0',0,1,IOPT,MODEL, 388 & WORK(KT1AMP0),WORK(KT2AMP0)) 389 RSPIM = .FALSE. 390 CALL CCRHSN(WORK(KOMEGA1),WORK(KOMEGA2),WORK(KT1AMP0), 391 & WORK(KT2AMP0),WORK(KEND1A),LWRK1A,'XXX') 392 CALL DSCAL(NT1AMX,5.0D0,WORK(KOMEGA1),1) 393 IF (.NOT.CCS) CALL DSCAL(NT2AMX,5.0D0,WORK(KOMEGA2),1) 394 CALL AROUND('result from CCRHSN (x 5):') 395 ELSE 396C ------------------------------------------------------------ 397C evaluate orbital relaxation and reorthogonalization contrib. 398C to the xi vector by finite difference on the vector function 399C (does only work for totally symmetric perturbations... 400C actually only tested without symmetry...) 401C ------------------------------------------------------------ 402 IF (ISYHOP.NE.1) THEN 403 WRITE (LUPRI,*) 'finite difference test of Xi vector not ' 404 WRITE (LUPRI,*) 'available for non-totally sym. '// 405 & 'perturbations.' 406 CALL QUIT( 407 & 'CC_FDXI only implemented for tot. sym. perturb.') 408 END IF 409 IF (IREAL.NE.1) THEN 410 WRITE (LUPRI,*) 'finite difference test of Xi vector not ' 411 WRITE (LUPRI,*) 'available for imaginary perturbations.' 412 CALL QUIT('CC_FDXI only implemented for real perturb.') 413 END IF 414 CALL CC_FDXI(WORK(KCMOHQ),WORK(KOMEGA1),WORK(KEND1),LWRK1) 415 CALL AROUND('fin. diff. orbital relaxation Xi vector:') 416 END IF 417 CALL CC_PRP(WORK(KOMEGA1),WORK(KOMEGA2),ISYHOP,1,N2VEC) 418 ELSE 419 WRITE (LUPRI,*) 'no orb. relax. or reorth. contribution to Xi.' 420 CALL DZERO(WORK(KOMEGA1),NT1AM(ISYHOP)) 421 IF (.NOT.CCS) CALL DZERO(WORK(KOMEGA2),NT2AM(ISYHOP)) 422 END IF 423 424C ---------------------------------------------------- 425C calculate the contribution from one-electron h^1: 426C ---------------------------------------------------- 427 IF (LABEL(1:8).NE.'DUMMYOP '.AND.LABEL(1:8).NE.'HAM0 ') THEN 428 IOPTCC2 = 0 429 IF (LORX) IOPTCC2 = 1 430 CALL CC_XKSI(WORK(KRHS1),LABEL,ISYHOP,IOPTCC2,DUMMY, 431 & WORK(KEND1),LWRK1) 432 CALL AROUND('one-electron h^1 contribution to Xi vector:') 433 CALL CC_PRP(WORK(KRHS1),WORK(KRHS2),ISYHOP,1,N2VEC) 434 CALL DAXPY(NT1AM(ISYHOP),ONE,WORK(KRHS1),1,WORK(KOMEGA1),1) 435 IF (.NOT.CCS) 436 & CALL DAXPY(NT2AM(ISYHOP),ONE,WORK(KRHS2),1,WORK(KOMEGA2),1) 437 END IF 438 439 CALL AROUND('num. orb.-relax. + analyt. one-electr. Xi vector:') 440 CALL CC_PRP(WORK(KOMEGA1),WORK(KOMEGA2),ISYHOP,1,N2VEC) 441 442C ------------------------------------------------- 443C read the result vector from CC_XIETA and compare: 444C ------------------------------------------------- 445 DO ITRAN = 1, NXETRAN 446 IOPT = 3 447 IVEC = IXETRAN(3,ITRAN) 448 WRITE (LUPRI,*) 'CC_XETST: IVEC = ',IVEC 449 CALL CC_RDRSP('O1 ',IVEC,ISYHOP,IOPT,MODEL, 450 & WORK(KRHS1),WORK(KRHS2)) 451 452 CALL AROUND('analytical Xi vector:') 453 CALL CC_PRP(WORK(KRHS1),WORK(KRHS2),ISYHOP,1,N2VEC) 454 455 CALL DAXPY(NT1AM(ISYHOP),-1.0d0,WORK(KOMEGA1),1,WORK(KRHS1),1) 456 IF (.NOT.CCS) 457 & CALL DAXPY(NT2AM(ISYHOP),-1.0d0,WORK(KOMEGA2),1,WORK(KRHS2),1) 458 459 WRITE (LUPRI,*) 'Norm of difference between analytical Xi ' 460 & // 'vector and the numerical result:' 461 WRITE (LUPRI,*) 'singles excitation part:', 462 & DSQRT(DDOT(NT1AM(ISYHOP),WORK(KRHS1),1,WORK(KRHS1),1)) 463 IF (.NOT.CCS) WRITE (LUPRI,*) 'double excitation part: ', 464 & DSQRT(DDOT(NT2AM(ISYHOP),WORK(KRHS2),1,WORK(KRHS2),1)) 465 466 WRITE (LUPRI,*) 'difference vector:' 467 CALL CC_PRP(WORK(KRHS1),WORK(KRHS2),ISYHOP,1,N2VEC) 468 469C ------------------------------------------------- 470C test first-order property result: 471C ------------------------------------------------- 472 IF (ISYHOP.EQ.1) THEN 473 IOPT = 3 474 IVEC = IXETRAN(3,ITRAN) 475 CALL CC_RDRSP('O1 ',IVEC,ISYHOP,IOPT,MODEL, 476 & WORK(KRHS1),WORK(KRHS2)) 477 IOPT = 3 478 IVEC = 0 479 ISYM0 = 1 480 CALL CC_RDRSP('L0',IVEC,ISYM0,IOPT,MODEL, 481 & WORK(KCTR1),WORK(KCTR2)) 482 PROPRSP = DDOT(NT1AMX+NT2AMX,WORK(KRHS1),1,WORK(KCTR1),1) 483 IOPER = IROPER(LABEL,ISYHOP) 484 IF (LORX.OR.LTWO .OR. IORDER.GT.1) THEN 485 ILSTETA = IETA1(LABEL,LORX,FREQ,ISYHOP) 486 PROPAVE = AVEX1(ILSTETA) 487 WRITE (LUPRI,*) '...average contribution taken from '// 488 & 'AVEX1...' 489 ELSE 490 FF = 1.0D0 491 CALL CC_AVE(PROPAVE,LABEL,WORK(KEND1),LWRK1,FF) 492 WRITE (LUPRI,*) '...average contribution calculated '// 493 & 'by CC_AVE...' 494 END IF 495 WRITE (LUPRI,*) ' OPERATOR : ',LABEL 496 WRITE (LUPRI,*) ' AVERAGE CONTRIBUTION :',PROPAVE 497 WRITE (LUPRI,*) ' RESPONSE CONTRIBUTION :',PROPRSP 498 WRITE (LUPRI,*) ' TOTAL ELECTRONIC CONTRIBUTION :', 499 & PROPRSP + PROPAVE 500 END IF 501 502 END DO 503 504 505 IF (ITEST.EQ.1) THEN 506 ! recalculate the response intermediates 507 IOPT = 3 508 CALL CC_RDRSP('R0',0,1,IOPT,MODEL, 509 & WORK(KT1AMP0),WORK(KT2AMP0)) 510 RSPIM = .TRUE. 511 CALL CCRHSN(WORK(KOMEGA1),WORK(KOMEGA2),WORK(KT1AMP0), 512 & WORK(KT2AMP0),WORK(KEND1A),LWRK1A,'XXX') 513 END IF 514 515*---------------------------------------------------------------------* 516* calculate the reference value for the eta vector and compare 517* with the results from the CC_XIETA routine: 518*---------------------------------------------------------------------* 519 ISYCTR = ILSTSYM(LISTL,IDLSTL) 520 ISYETA = MULD2H(ISYCTR,ISYHOP) 521 522 KKAPPA = 1 523 KRMAT = KKAPPA + 2*NALLAI(ISYHOP) 524 KCMOPQ = KRMAT + N2BST(ISYHOP) 525 KCMOHQ = KCMOPQ + NGLMDT(ISYHOP) 526 KLEFT1 = KCMOHQ + NGLMDT(ISYHOP) 527 KLEFT2 = KLEFT1 + NT1AM(ISYCTR) 528 KRHO1 = KLEFT2 + NT2AM(ISYCTR) 529 KRHO2 = KRHO1 + NT1AM(ISYETA) 530 KETA1 = KRHO2 + NT2AM(ISYETA) 531 KETA2 = KETA1 + NT1AM(ISYETA) 532 KEND1 = KETA2 + NT2AM(ISYETA) 533 if (CCR12) then 534 keta12 = kend1 535 kend1 = keta12 + ntr12am(isyeta) 536 end if 537 LWRK1 = LWORK - KEND1 538 539 IF (LWRK1.LT.0 .OR. LWRK1A.LT.0) THEN 540 CALL QUIT('Insufficient memory in CC_XETST.') 541 END IF 542 543 IOPT = 3 544 Call CC_RDRSP(LISTL,IDLSTL,ISYCTR,IOPT,MODEL, 545 & WORK(KLEFT1),WORK(KLEFT2)) 546 547 IF (LORX.OR.LTWO) THEN 548 IF (ITEST.EQ.1) THEN 549C ------------------------------------------------------------ 550C call the left transformation to calculate the reference 551C result for the eta vector: 552C ------------------------------------------------------------ 553 CALL DCOPY(NT1AM(ISYCTR),WORK(KLEFT1),1,WORK(KRHO1),1) 554 CALL DCOPY(NT2AM(ISYCTR),WORK(KLEFT2),1,WORK(KRHO2),1) 555 CALL CC_ATRR(0.0D0,ISYCTR,-1,WORK(KRHO1),LWRK1,.FALSE.,DUMMY, 556 & APROXR12,.FALSE.) 557 IF (LISTL(1:2).EQ.'L0') THEN 558 CALL CC_ETA(WORK(KETA1),WORK(KEND1),LWRK1) 559 CALL DAXPY(NT1AMX,1.0D0,WORK(KETA1),1,WORK(KRHO1),1) 560 IF (.NOT.CCS) 561 & CALL DAXPY(NT2AMX,1.0D0,WORK(KETA2),1,WORK(KRHO2),1) 562 END IF 563 CALL DSCAL(NT1AMX,5.0D0,WORK(KRHO1),1) 564 IF (.NOT.CCS) CALL DSCAL(NT2AMX,5.0D0,WORK(KRHO2),1) 565 CALL AROUND('result of CC_ETA & CC_LHTR (x 5):') 566 ELSE 567C --------------------------------------------------------- 568C evaluate orbital relaxation and reorthogonalization contrib. 569C to the eta vector by finite difference on the jacobian 570C left transformation + zeroth-order eta vector. 571C (does only work for totally symmetric perturbations... 572C actually only tested without symmetry...) 573C --------------------------------------------------------- 574 IF (ISYETA.NE.1) THEN 575 CALL QUIT 576 * ('CC_FDETA only implemented for total symmetric VECL.') 577 END IF 578 IF (IREAL.NE.1) THEN 579 WRITE (LUPRI,*) 'finite difference test of Xi vector '// 580 * 'not available for imaginary' 581 WRITE (LUPRI,*) 'perturbations.' 582 CALL QUIT('CC_FDETA only implemented for real perturb.') 583 END IF 584 IF (ISYHOP.NE.1) THEN 585 WRITE (LUPRI,*) 'finite difference test of Xi '// 586 & 'vector not available' 587 WRITE (LUPRI,*) 'for non-totally symmetric perturbations.' 588 CALL QUIT('CC_FDETA only implemented for tot. '// 589 & 'sym. perturb.') 590 END IF 591 CALL CC_FDETA(WORK(KCMOHQ),LISTL,WORK(KLEFT1),WORK(KLEFT2), 592 & WORK(KRHO1),WORK(KEND1),LWRK1) 593 CALL AROUND('fin. diff. orbital relaxation Eta vector:') 594 END IF 595 CALL CC_PRP(WORK(KRHO1),WORK(KRHO2),ISYETA,1,N2VEC) 596 ELSE 597 WRITE (LUPRI,*) 'no orb. relax. or reorth. '// 598 & 'contribution to Eta.' 599 CALL DZERO(WORK(KRHO1),NT1AM(ISYETA)) 600 IF (.NOT.CCS) CALL DZERO(WORK(KRHO2),NT2AM(ISYETA)) 601 END IF 602 603C ---------------------------------------------------- 604C calculate the contribution from one-electron h^1: 605C ---------------------------------------------------- 606 IF (LABEL(1:8).NE.'DUMMYOP '.AND.LABEL(1:8).NE.'HAM0 ') THEN 607 IOPTCC2 = 0 608 IF (LORX) IOPTCC2 = 1 609 CALL CC_ETAC(ISYHOP,LABEL,WORK(KETA1),LISTL,IDLSTL,IOPTCC2, 610 & DUMMY,WORK(KEND1),LWRK1) 611 CALL AROUND('one-electron h^1 contribution to Eta vector:') 612 CALL CC_PRP(WORK(KETA1),WORK(KETA2),ISYETA,1,N2VEC) 613 CALL DAXPY(NT1AM(ISYETA),ONE,WORK(KETA1),1,WORK(KRHO1),1) 614 IF (.NOT.CCS) 615 & CALL DAXPY(NT2AM(ISYETA),ONE,WORK(KETA2),1,WORK(KRHO2),1) 616 END IF 617 618 CALL AROUND('num. orb.-relax. + analyt. one-electr. Eta vector:') 619 CALL CC_PRP(WORK(KRHO1),WORK(KRHO2),ISYETA,1,N2VEC) 620 621C ------------------------------------------------- 622C read the result vector from CC_XIETA and compare: 623C ------------------------------------------------- 624 DO ITRAN = 1, NXETRAN 625 IOPT = 3 626 IVEC = IXETRAN(4,ITRAN) 627 WRITE (LUPRI,*) 'CC_XETST: IVEC = ',IVEC 628 CALL CC_RDRSP('X1 ',IVEC,ISYETA,IOPT,MODEL, 629 & WORK(KETA1),WORK(KETA2)) 630 631 CALL AROUND('analytical Eta vector:') 632 CALL CC_PRP(WORK(KETA1),WORK(KETA2),ISYETA,1,N2VEC) 633 634 CALL DAXPY(NT1AM(ISYETA),-1.0d0,WORK(KRHO1),1,WORK(KETA1),1) 635 IF (.NOT.CCS) 636 & CALL DAXPY(NT2AM(ISYETA),-1.0d0,WORK(KRHO2),1,WORK(KETA2),1) 637 638 WRITE (LUPRI,*) 'Norm of difference between analytical ' 639 & // 'Eta vector and the numerical result:' 640 WRITE (LUPRI,*) 'singles excitation part:', 641 & DSQRT(DDOT(NT1AM(ISYETA),WORK(KETA1),1,WORK(KETA1),1)) 642 IF (.NOT.CCS) WRITE (LUPRI,*) 'double excitation part: ', 643 & DSQRT(DDOT(NT2AM(ISYETA),WORK(KETA2),1,WORK(KETA2),1)) 644 645 WRITE (LUPRI,*) 'difference vector:' 646 CALL CC_PRP(WORK(KETA1),WORK(KETA2),ISYETA,1,N2VEC) 647 END DO 648 649 RETURN 650 END 651*=====================================================================* 652*=====================================================================* 653 SUBROUTINE CC_FDXI(CMO1,OMEGADIF,WORK,LWORK) 654C 655C--------------------------------------------------------------------- 656C Test routine for calculating the orbital relaxation contribution 657C to the CC Xi vector by finite difference on the vector function. 658C Ch. Haettig, januar 1999 659C--------------------------------------------------------------------- 660C 661#include "implicit.h" 662#include "priunit.h" 663#include "dummy.h" 664#include "maxorb.h" 665#include "iratdef.h" 666#include "ccorb.h" 667#include "aovec.h" 668#include "ccsdinp.h" 669#include "cclr.h" 670#include "ccsdsym.h" 671#include "ccsdio.h" 672#include "inftap.h" 673#include "ccinftap.h" 674#include "leinf.h" 675#include "ccfield.h" 676C 677 DIMENSION WORK(LWORK), OMEGADIF(*), CMO1(*) 678 PARAMETER (DELTA = 1.0D-07, DELINV = 1.0D+7) 679 PARAMETER (ONE = 1.0d0, ZERO = 0.0d0, TWO = 2.0d0, XHALF = 0.5D0) 680 LOGICAL NONHFSAVE,LHTF 681 CHARACTER MODEL*10 682C 683 IF (IPRINT.GT.10) THEN 684 CALL AROUND( 'IN CC_FDXI : MAKING FINITE DIFF. CC Xi Vector') 685 END IF 686C 687 IF (CCR12) CALL QUIT('Finite-difference Xi-vector for CCR12 '// 688 & 'not adapted') 689C 690 IF (CC2) THEN 691 IF (NFIELD.GT.0 .AND. (.NOT.NONHF)) THEN 692 CALL QUIT('CC_FDXI incompatible with relaxed finite fields.') 693 END IF 694 NONHFSAVE = NONHF 695 NONHF = .TRUE. 696 END IF 697C 698C---------------------------- 699C Work space allocations. 700C---------------------------- 701C 702 ISYMTR = 1 703 ISYMOP = 1 704 ISYM0 = 1 705C 706 NTAMP = NT1AMX + NT2AMX 707 IF (CCS) NTAMP = NT1AMX 708C 709 KCMOREF = 1 710 KCMO = KCMOREF + NLAMDS 711 KOMREF1 = KCMO + NLAMDS 712 KOMREF2 = KOMREF1 + NT1AMX 713 KXIMAT = KOMREF2 + NT2AMX 714 KFCKDIF = KXIMAT + NTAMP * NLAMDS 715 KFCKREF = KFCKDIF + N2BST(ISYM0) 716 KFCKMAT = KFCKREF + N2BST(ISYM0) 717 KEND0 = KFCKMAT + N2BST(ISYM0) * NLAMDS 718C 719 KT1AMP0 = KEND0 720 KOMEGA1 = KT1AMP0 + NT1AMX 721 KOMEGA2 = KOMEGA1 + NT1AMX 722 KT2AMP0 = KOMEGA2 + MAX(NT2AMX,2*NT2ORT(1),NT2AO(1)) 723 KSCR2 = KT2AMP0 + NT2AMX 724 KEND1 = KSCR2 + NT2AMX + NT1AMX 725 LWRK1 = LWORK - KEND1 726C 727 IF (LWRK1.LT.0 ) THEN 728 WRITE(LUPRI,*) 'Too little work space in CC_FDXI ' 729 WRITE(LUPRI,*) 'AVAILABLE: LWORK = ',LWORK 730 WRITE(LUPRI,*) 'NEEDED (AT LEAST) = ',KEND1 731 CALL QUIT('TOO LITTLE WORKSPACE IN CC_FDXI ') 732 ENDIF 733C 734C--------------------- 735C Initializations. 736C--------------------- 737C 738 X1 = 0.0D0 739 X2 = 0.0D0 740 CALL DZERO(OMEGADIF,NTAMP) 741C 742C--------------------------------------------- 743C Read the reference CMO vector from disk: 744C--------------------------------------------- 745C 746C 747 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY, 748 & .FALSE.) 749 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 750 READ(LUSIFC) 751 READ(LUSIFC) 752 READ(LUSIFC) (WORK(KCMOREF-1+I),I=1,NLAMDS) 753 CALL GPCLOSE(LUSIFC,'KEEP') 754C 755 CALL DZERO(WORK(KT1AMP0),NT1AMX) 756 LHTF = .FALSE. 757 CALL CCSD_IAJB(WORK(KT2AMP0),WORK(KT1AMP0),LHTF, 758 & .FALSE.,.FALSE.,WORK(KEND1),LWRK1) 759 REWIND(LUIAJB) 760 CALL WRITI(LUIAJB,IRAT*NT2AMX,WORK(KT2AMP0)) 761C 762C---------------------------- 763C Read the CC amplitudes: 764C---------------------------- 765C 766 IOPT = 3 767 CALL CC_RDRSP('R0 ',0,1,IOPT,MODEL,WORK(KT1AMP0),WORK(KT2AMP0)) 768C 769C--------------------------------------------------------- 770C Calculate reference vector function and fock matrix: 771C--------------------------------------------------------- 772C 773 RSPIM = .FALSE. 774 CALL CCRHSN(WORK(KOMEGA1),WORK(KOMEGA2), 775 & WORK(KT1AMP0),WORK(KT2AMP0), 776 & WORK(KEND1),LWRK1,'XXX') 777C 778 IF (CCS) CALL DZERO(WORK(KOMEGA2),NT2AMX) 779C 780 OMEGA1N = DDOT(NT1AMX,WORK(KOMEGA1),1,WORK(KOMEGA1),1) 781 OMEGA2N = DDOT(NT2AMX,WORK(KOMEGA2),1,WORK(KOMEGA2),1) 782C 783 IF (IPRINT.GT.10) THEN 784 WRITE (LUPRI,*) 'CC_FDXI: norm of reference OMEGA1:',OMEGA1N 785 WRITE (LUPRI,*) 'CC_FDXI: norm of reference OMEGA2:',OMEGA2N 786 END IF 787C 788 CALL DCOPY(NT1AMX,WORK(KOMEGA1),1,WORK(KOMREF1),1) 789 CALL DCOPY(NT2AMX,WORK(KOMEGA2),1,WORK(KOMREF2),1) 790C 791 LUFOCK = -1 792 CALL GPOPEN(LUFOCK,'CC_FCKH','OLD',' ','UNFORMATTED',IDUMMY, 793 & .FALSE.) 794 REWIND(LUFOCK) 795 READ(LUFOCK) (WORK(KFCKREF-1+I),I=1,N2BST(ISYM0)) 796 CALL GPCLOSE(LUFOCK,'KEEP') 797 798*----------------------------------------------------------------------* 799* Calculate the derivative of the vector function with respect 800* to the CMO vector by finite difference: 801*----------------------------------------------------------------------* 802 803 DO IDXDIF = 1, NLAMDS 804C 805C ------------------------------- 806C add finite displacement to CMO: 807C ------------------------------- 808C 809 CALL DCOPY(NLAMDS,WORK(KCMOREF),1,WORK(KCMO),1) 810 WORK(KCMO-1 + IDXDIF) = WORK(KCMO-1 + IDXDIF) + DELTA 811C 812 CMON = DDOT(NLAMDS,WORK(KCMO),1,WORK(KCMO),1) 813C 814 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY, 815 & .FALSE.) 816 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 817 READ(LUSIFC) 818 READ(LUSIFC) 819 WRITE(LUSIFC) (WORK(KCMO-1+I),I=1,NLAMDS) 820 CALL GPCLOSE(LUSIFC,'KEEP') 821C 822 CALL DZERO(WORK(KT1AMP0),NT1AMX) 823 LHTF = .FALSE. 824 CALL CCSD_IAJB(WORK(KT2AMP0),WORK(KT1AMP0), 825 & LHTF,.FALSE.,.FALSE.,WORK(KEND1),LWRK1) 826 REWIND(LUIAJB) 827 CALL WRITI(LUIAJB,IRAT*NT2AMX,WORK(KT2AMP0)) 828C 829C ------------------------------ 830C calculate the vector function: 831C ------------------------------ 832C 833 IOPT = 3 834 CALL CC_RDRSP('R0 ',0,1,IOPT,MODEL,WORK(KT1AMP0),WORK(KT2AMP0)) 835 836 RSPIM = .FALSE. 837 CALL CCRHSN(WORK(KOMEGA1),WORK(KOMEGA2), 838 & WORK(KT1AMP0),WORK(KT2AMP0), 839 & WORK(KEND1),LWRK1,'XXX') 840C 841 IF (CCS) CALL DZERO(WORK(KOMEGA2),NT2AMX) 842C 843 CALL GPOPEN(LUFOCK,'CC_FCKH','OLD',' ','UNFORMATTED',IDUMMY, 844 & .FALSE.) 845 REWIND(LUFOCK) 846 READ(LUFOCK) (WORK(KFCKDIF-1+I),I=1,N2BST(ISYM0)) 847 CALL GPCLOSE(LUFOCK,'KEEP') 848C 849C -------------------------------------------------- 850C construct the row nb. IDXDIF of the result matrix: 851C -------------------------------------------------- 852C 853 OMEGA1N = DDOT(NT1AMX,WORK(KOMEGA1),1,WORK(KOMEGA1),1) 854 OMEGA2N = DDOT(NT2AMX,WORK(KOMEGA2),1,WORK(KOMEGA2),1) 855C 856 IF (IPRINT.GT.10) THEN 857 WRITE (LUPRI,*) 'CMO index:',IDXDIF 858 WRITE (LUPRI,*) 'Norm of OMEGA1: ',OMEGA1N 859 WRITE (LUPRI,*) 'Norm of OMEGA2: ',OMEGA2N 860 END IF 861C 862 CALL DAXPY(NT1AMX,-1.0D0,WORK(KOMREF1),1,WORK(KOMEGA1),1) 863 CALL DAXPY(NT2AMX,-1.0D0,WORK(KOMREF2),1,WORK(KOMEGA2),1) 864 CALL DSCAL(NT1AMX,DELINV,WORK(KOMEGA1),1) 865 CALL DSCAL(NT2AMX,DELINV,WORK(KOMEGA2),1) 866C 867 IF (IPRINT.GT.100) THEN 868 CALL CC_PRP(WORK(KOMEGA1),WORK(KOMEGA2),1,1,1) 869 WRITE (LUPRI,*) 'CMO1(index) = ', CMO1(IDXDIF) 870 CALL DAXPY(NT1AMX,CMO1(IDXDIF),WORK(KOMEGA1),1, 871 & OMEGADIF(1),1) 872 CALL DAXPY(NT2AMX,CMO1(IDXDIF),WORK(KOMEGA2),1, 873 & OMEGADIF(1+NT1AMX),1) 874 WRITE (LUPRI,*) 'accumulated result vector:' 875 CALL CC_PRP(OMEGADIF,OMEGADIF(NT1AMX+1),1,1,1) 876 ENDIF 877C 878 KOFF1 = KXIMAT + NTAMP*(IDXDIF-1) 879 KOFF2 = KOFF1 + NT1AMX 880 CALL DCOPY(NT1AMX,WORK(KOMEGA1),1,WORK(KOFF1),1) 881 CALL DCOPY(NT2AMX,WORK(KOMEGA2),1,WORK(KOFF2),1) 882C 883 X1 = X1 + DDOT(NT1AMX,WORK(KOMEGA1),1,WORK(KOMEGA1),1) 884 X2 = X2 + DDOT(NT2AMX,WORK(KOMEGA2),1,WORK(KOMEGA2),1) 885C 886 CALL DAXPY(N2BST(ISYM0),-1.0D0,WORK(KFCKREF),1,WORK(KFCKDIF),1) 887 CALL DSCAL(N2BST(ISYM0),1.0D0/DELTA,WORK(KFCKDIF),1) 888C 889 KOFF1 = KFCKMAT + N2BST(ISYM0)*(IDXDIF-1) 890 CALL DCOPY(N2BST(ISYM0),WORK(KFCKDIF),1,WORK(KOFF1),1) 891C 892 END DO 893C 894 IF (IPRINT .GT. 10) THEN 895 WRITE(LUPRI,*) ' ' 896 WRITE(LUPRI,*) '** FINITE DIFF WITH DELTA ',DELTA, '**' 897 WRITE(LUPRI,*) ' ' 898 IF (IPRINT .GT. 20) THEN 899 CALL AROUND( 'FINITE DIFF. CC Xi-Matrix - 10 & 20 PART' ) 900 CALL OUTPUT(WORK(KXIMAT),1,NTAMP,1,NLAMDS,NTAMP,NLAMDS,1, 901 & LUPRI) 902 END IF 903 XN = X1 + X2 904 WRITE(LUPRI,*) ' ' 905 WRITE(LUPRI,*) ' NORM OF 10 PART OF FD. Xi-mat.: ', SQRT(X1) 906 WRITE(LUPRI,*) ' NORM OF 20 PART OF FD. Xi-mat.: ', SQRT(X2) 907 WRITE(LUPRI,*) ' NORM OF COMPLETE FD. Xi-mat.: ', SQRT(XN) 908 ENDIF 909C 910C------------------------------------------- 911C calculate Xi-matrix times CMO1 vector: 912C------------------------------------------- 913C 914 CALL DGEMV('N',NTAMP,NLAMDS,ONE,WORK(KXIMAT),NTAMP,CMO1,1, 915 * ZERO,OMEGADIF,1) 916C 917C----------------------------- 918C scale diagonal with 1/2: 919C----------------------------- 920 CALL CCLR_DIASCL(OMEGADIF(NT1AM(1)+1),XHALF,1) 921 922C 923C---------------------------------------------- 924C restore the reference CMO vector on disk: 925C---------------------------------------------- 926C 927 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY, 928 & .FALSE.) 929 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 930 READ(LUSIFC) 931 READ(LUSIFC) 932 WRITE(LUSIFC) (WORK(KCMOREF-1+I),I=1,NLAMDS) 933 CALL GPCLOSE(LUSIFC,'KEEP') 934 IF (CC2) THEN 935 NONHF = NONHFSAVE 936 END IF 937C 938C------------------------------------------------------------------ 939C make sure that all intermediates on file are calculated with 940C the reference CMO vector: 941C------------------------------------------------------------------ 942C 943C 944 CALL DZERO(WORK(KT1AMP0),NT1AMX) 945 LHTF = .FALSE. 946 CALL CCSD_IAJB(WORK(KT2AMP0),WORK(KT1AMP0),LHTF, 947 & .FALSE.,.FALSE.,WORK(KEND1),LWRK1) 948 REWIND(LUIAJB) 949 CALL WRITI(LUIAJB,IRAT*NT2AMX,WORK(KT2AMP0)) 950C 951 RSPIM = .TRUE. 952 CALL CCRHSN(WORK(KOMEGA1),WORK(KOMEGA2), 953 & WORK(KT1AMP0),WORK(KT2AMP0), 954 & WORK(KEND1),LWRK1,'XXX') 955C 956 IF (IPRINT.GT.10) THEN 957 CALL AROUND(' END OF CC_FDXI ') 958 END IF 959C 960 RETURN 961 END 962*=====================================================================* 963*======================================================================* 964 SUBROUTINE CC_FDETA(CMO1,LISTL,VECL1,VECL2,RHODIF,WORK,LWORK) 965C 966C---------------------------------------------------------------------- 967C Test routine for calculating the orbital relaxation contribution 968C to the CC ETA{O} vector by finite difference on the jacob. left tran. 969C Ch. Haettig, januar 1999 970C---------------------------------------------------------------------- 971C 972#include "implicit.h" 973#include "priunit.h" 974#include "dummy.h" 975#include "maxorb.h" 976#include "iratdef.h" 977#include "ccorb.h" 978#include "aovec.h" 979#include "ccsdinp.h" 980#include "cclr.h" 981#include "ccsdsym.h" 982#include "ccsdio.h" 983#include "ccinftap.h" 984#include "inftap.h" 985#include "leinf.h" 986#include "ccfield.h" 987C 988C------------------------------------------------------------ 989C the displacement for the finite difference calculation: 990C------------------------------------------------------------ 991 PARAMETER (DELTA = 1.0D-07, DELINV = 1.0D+7) 992C------------------------------------------------------------ 993C 994 DIMENSION WORK(LWORK), RHODIF(*), CMO1(*), VECL1(*), VECL2(*) 995 PARAMETER (ONE = 1.0d0, ZERO = 0.0d0, TWO = 2.0d0, XHALF = 0.5D0) 996 CHARACTER MODEL*10, LISTL*(*) 997 LOGICAL NONHFSAVE,LHTF 998C 999 IF (IPRINT.GT.10) THEN 1000 CALL AROUND('IN CC_FDETA : MAKING FIN. DIFF. CC Eta{O} Vector') 1001 END IF 1002C 1003 IF (CC2) THEN 1004 IF (NFIELD.GT.0 .AND. (.NOT.NONHF)) THEN 1005 CALL QUIT('CC_FDETA incompatible with relaxed finite fields.') 1006 END IF 1007 NONHFSAVE = NONHF 1008 NONHF = .TRUE. 1009 END IF 1010C 1011C---------------------------- 1012C Work space allocations. 1013C---------------------------- 1014C 1015 ISYMTR = 1 1016 ISYMOP = 1 1017 ISYM0 = 1 1018C 1019 NTAMP = NT1AMX + NT2AMX 1020 IF (CCS) NTAMP = NT1AMX 1021C 1022 KCMOREF = 1 1023 KCMO = KCMOREF + NLAMDS 1024 KRHREF1 = KCMO + NLAMDS 1025 KRHREF2 = KRHREF1 + NT1AMX 1026 KXIMAT = KRHREF2 + NT2AMX 1027 KEND0 = KXIMAT + NTAMP * NLAMDS 1028 LWRK0 = LWORK - KEND0 1029C 1030 KT1AMP0 = KEND0 1031 KOMEGA1 = KT1AMP0 + NT1AMX 1032 KOMEGA2 = KOMEGA1 + NT1AMX 1033 KT2AMP0 = KOMEGA2 + MAX(NT2AMX,2*NT2ORT(1),NT2AO(1)) 1034 KSCR2 = KT2AMP0 + NT2AMX 1035 KEND1 = KSCR2 + NT2AMX + NT1AMX 1036 LWRK1 = LWORK - KEND1 1037C 1038 KRHO1 = KEND0 1039 KRHO2 = KRHO1 + NT1AMX 1040 KEND1A = KRHO2 + NT2AMX 1041 LWRK1A = LWORK - KEND1A 1042C 1043 KETA1 = KEND1A 1044 KETA2 = KETA1 + NT1AMX 1045 KEND1B = KETA2 + NT2AMX 1046 if (CCR12) then 1047 keta12 = kend1b 1048 kend1b = keta12 + ntr12am(1) 1049 end if 1050 LWRK1B = LWORK - KEND1B 1051C 1052 IF ( LWRK1.LT.0 .OR. LWRK1B.LT.0 .OR. LWRK1A.LT.0) THEN 1053 WRITE(LUPRI,*) 'Too little work space in CC_FDETA ' 1054 WRITE(LUPRI,*) 'AVAILABLE: LWORK = ',LWORK 1055 WRITE(LUPRI,*) 'NEEDED (AT LEAST) = ',MAX(KEND1,KEND1B) 1056 CALL QUIT('TOO LITTLE WORKSPACE IN CC_FDETA ') 1057 ENDIF 1058C 1059C--------------------- 1060C Initializations. 1061C--------------------- 1062C 1063 X1 = 0.0D0 1064 X2 = 0.0D0 1065 CALL DZERO(RHODIF,NTAMP) 1066C 1067C--------------------------------------------- 1068C Read the reference CMO vector from disk: 1069C--------------------------------------------- 1070C 1071 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY, 1072 & .FALSE.) 1073 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 1074 READ(LUSIFC) 1075 READ(LUSIFC) 1076 READ(LUSIFC) (WORK(KCMOREF-1+I),I=1,NLAMDS) 1077 CALL GPCLOSE(LUSIFC,'KEEP') 1078C 1079C------------------------------------------------------------------- 1080C make sure that the correct response intermediates are on disc: 1081C------------------------------------------------------------------- 1082C 1083 CALL DZERO(WORK(KT1AMP0),NT1AMX) 1084 LHTF = .FALSE. 1085 CALL CCSD_IAJB(WORK(KT2AMP0),WORK(KT1AMP0),LHTF, 1086 & .FALSE.,.FALSE.,WORK(KEND1),LWRK1) 1087 REWIND(LUIAJB) 1088 CALL WRITI(LUIAJB,IRAT*NT2AMX,WORK(KT2AMP0)) 1089C 1090 IOPT = 3 1091 CALL CC_RDRSP('R0 ',0,1,IOPT,MODEL,WORK(KT1AMP0),WORK(KT2AMP0)) 1092C 1093 RSPIM = .TRUE. 1094 CALL CCRHSN(WORK(KRHO1),WORK(KRHO2), 1095 & WORK(KT1AMP0),WORK(KT2AMP0), 1096 & WORK(KEND1),LWRK1,'XXX') 1097C 1098C------------------------------------------------------------------- 1099C calculate the reference vector: 1100C------------------------------------------------------------------- 1101C 1102 CALL DCOPY(NT1AMX,VECL1,1,WORK(KRHO1),1) 1103 IF (.NOT.CCS) CALL DCOPY(NT2AMX,VECL2,1,WORK(KRHO2),1) 1104C 1105 CALL CC_ATRR(0.0D0,ISYM0,-1,WORK(KRHO1),LWRK0,.FALSE.,DUMMY, 1106 & APROXR12,.FALSE.) 1107C 1108 IF (LISTL(1:2).EQ.'L0') THEN 1109 CALL CC_ETA(WORK(KETA1),WORK(KEND1B),LWRK1B) 1110 CALL DAXPY(NT1AMX,ONE,WORK(KETA1),1,WORK(KRHO1),1) 1111 CALL DAXPY(NT2AMX,ONE,WORK(KETA2),1,WORK(KRHO2),1) 1112 END IF 1113C 1114 IF (CCS) CALL DZERO(WORK(KRHO2),NT2AMX) 1115C 1116 RHO1N = DDOT(NT1AMX,WORK(KRHO1),1,WORK(KRHO1),1) 1117 RHO2N = DDOT(NT2AMX,WORK(KRHO2),1,WORK(KRHO2),1) 1118C 1119 IF (IPRINT.GT.10) THEN 1120 WRITE (LUPRI,*) 'CC_FDETA: norm of reference RHO1:',RHO1N 1121 WRITE (LUPRI,*) 'CC_FDETA: norm of reference RHO2:',RHO2N 1122 END IF 1123C 1124 CALL DCOPY(NT1AMX,WORK(KRHO1),1,WORK(KRHREF1),1) 1125 CALL DCOPY(NT2AMX,WORK(KRHO2),1,WORK(KRHREF2),1) 1126 1127*----------------------------------------------------------------------* 1128* Calculate the derivative of the vector function with respect 1129* to the CMO vector by finite difference: 1130*----------------------------------------------------------------------* 1131 1132 DO IDXDIF = 1, NLAMDS 1133C 1134C ------------------------------- 1135C add finite displacement to CMO: 1136C ------------------------------- 1137C 1138 CALL DCOPY(NLAMDS,WORK(KCMOREF),1,WORK(KCMO),1) 1139 WORK(KCMO-1 + IDXDIF) = WORK(KCMO-1 + IDXDIF) + DELTA 1140C 1141 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY, 1142 & .FALSE.) 1143 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 1144 READ(LUSIFC) 1145 READ(LUSIFC) 1146 WRITE(LUSIFC) (WORK(KCMO-1+I),I=1,NLAMDS) 1147 CALL GPCLOSE(LUSIFC,'KEEP') 1148C 1149C ------------------------------------- 1150C calculate new response intermediates: 1151C ------------------------------------- 1152C 1153 CALL DZERO(WORK(KT1AMP0),NT1AMX) 1154 LHTF = .FALSE. 1155 CALL CCSD_IAJB(WORK(KT2AMP0),WORK(KT1AMP0), 1156 & LHTF,.FALSE.,.FALSE.,WORK(KEND1),LWRK1) 1157 REWIND(LUIAJB) 1158 CALL WRITI(LUIAJB,IRAT*NT2AMX,WORK(KT2AMP0)) 1159C 1160 IOPT = 3 1161 CALL CC_RDRSP('R0 ',0,1,IOPT,MODEL,WORK(KT1AMP0),WORK(KT2AMP0)) 1162C 1163 RSPIM = .TRUE. 1164 CALL CCRHSN(WORK(KRHO1),WORK(KRHO2), 1165 & WORK(KT1AMP0),WORK(KT2AMP0), 1166 & WORK(KEND1),LWRK1,'XXX') 1167C 1168C --------------------------------- 1169C calculate the transformed vector: 1170C --------------------------------- 1171C 1172 CALL DCOPY(NT1AMX,VECL1,1,WORK(KRHO1),1) 1173 IF (.NOT.CCS) CALL DCOPY(NT2AMX,VECL2,1,WORK(KRHO2),1) 1174C 1175 CALL CC_ATRR(0.0D0,ISYM0,-1,WORK(KRHO1),LWRK0,.FALSE.,DUMMY, 1176 & APROXR12,.FALSE.) 1177C 1178 IF (LISTL(1:2).EQ.'L0') THEN 1179 CALL CC_ETA(WORK(KETA1),WORK(KEND1B),LWRK1B) 1180 CALL DAXPY(NT1AMX,ONE,WORK(KETA1),1,WORK(KRHO1),1) 1181 CALL DAXPY(NT2AMX,ONE,WORK(KETA2),1,WORK(KRHO2),1) 1182 END IF 1183C 1184 IF (CCS) CALL DZERO(WORK(KRHO2),NT2AMX) 1185C 1186C -------------------------------------------------- 1187C construct the row nb. IDXDIF of the result matrix: 1188C -------------------------------------------------- 1189C 1190 RHO1N = DDOT(NT1AMX,WORK(KRHO1),1,WORK(KRHO1),1) 1191 RHO2N = DDOT(NT2AMX,WORK(KRHO2),1,WORK(KRHO2),1) 1192C 1193 IF (IPRINT.GT.10) THEN 1194 WRITE (LUPRI,*) 'CMO index:',IDXDIF 1195 WRITE (LUPRI,*) 'Norm of RHO1: ',RHO1N 1196 WRITE (LUPRI,*) 'Norm of RHO2: ',RHO2N 1197 END IF 1198C 1199 CALL DAXPY(NT1AMX,-1.0D0,WORK(KRHREF1),1,WORK(KRHO1),1) 1200 CALL DAXPY(NT2AMX,-1.0D0,WORK(KRHREF2),1,WORK(KRHO2),1) 1201 CALL DSCAL(NT1AMX,DELINV,WORK(KRHO1),1) 1202 CALL DSCAL(NT2AMX,DELINV,WORK(KRHO2),1) 1203C 1204 IF (IPRINT.GT.100) THEN 1205 CALL CC_PRP(WORK(KRHO1),WORK(KRHO2),1,1,1) 1206 WRITE (LUPRI,*) 'CMO1(index) = ', CMO1(IDXDIF) 1207 CALL DAXPY(NT1AMX,CMO1(IDXDIF),WORK(KRHO1),1, 1208 & RHODIF(1),1) 1209 CALL DAXPY(NT2AMX,CMO1(IDXDIF),WORK(KRHO2),1, 1210 & RHODIF(1+NT1AMX),1) 1211 WRITE (LUPRI,*) 'accumulated result vector:' 1212 CALL CC_PRP(RHODIF,RHODIF(NT1AMX+1),1,1,1) 1213 ENDIF 1214C 1215 KOFF1 = KXIMAT + NTAMP*(IDXDIF-1) 1216 KOFF2 = KOFF1 + NT1AMX 1217 CALL DCOPY(NT1AMX,WORK(KRHO1),1,WORK(KOFF1),1) 1218 CALL DCOPY(NT2AMX,WORK(KRHO2),1,WORK(KOFF2),1) 1219C 1220 X1 = X1 + DDOT(NT1AMX,WORK(KRHO1),1,WORK(KRHO1),1) 1221 X2 = X2 + DDOT(NT2AMX,WORK(KRHO2),1,WORK(KRHO2),1) 1222C 1223 END DO 1224C 1225 IF (IPRINT.GT.10) THEN 1226 WRITE(LUPRI,*) ' ' 1227 WRITE(LUPRI,*) '** FINITE DIFF WITH DELTA ',DELTA, '**' 1228 WRITE(LUPRI,*) ' ' 1229 IF (IPRINT .GT. 20) THEN 1230 CALL AROUND( 'FINITE DIFF. CC Eta-Matrix - 10 & 20 PART' ) 1231 CALL OUTPUT(WORK(KXIMAT),1,NTAMP,1,NLAMDS,NTAMP,NLAMDS,1, 1232 & LUPRI) 1233 ENDIF 1234 XN = X1 + X2 1235 WRITE(LUPRI,*) ' ' 1236 WRITE(LUPRI,*) ' NORM OF 10 PART OF FD. Eta-mat.: ', SQRT(X1) 1237 WRITE(LUPRI,*) ' NORM OF 20 PART OF FD. Eta-mat.: ', SQRT(X2) 1238 WRITE(LUPRI,*) ' NORM OF COMPLETE FD. Eta-mat.: ', SQRT(XN) 1239 END IF 1240C 1241C------------------------------------------- 1242C calculate Eta-matrix times CMO1 vector: 1243C------------------------------------------- 1244C 1245 CALL DGEMV('N',NTAMP,NLAMDS,ONE,WORK(KXIMAT),NTAMP,CMO1,1, 1246 * ZERO,RHODIF,1) 1247C 1248C---------------------------------------------- 1249C restore the reference CMO vector on disk: 1250C---------------------------------------------- 1251C 1252 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY, 1253 & .FALSE.) 1254 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 1255 READ(LUSIFC) 1256 READ(LUSIFC) 1257 WRITE(LUSIFC) (WORK(KCMOREF-1+I),I=1,NLAMDS) 1258 CALL GPCLOSE(LUSIFC,'KEEP') 1259C 1260C------------------------------------------------------------------ 1261C make sure that all intermediates on file are calculated with 1262C the reference CMO vector: 1263C------------------------------------------------------------------ 1264C 1265 CALL DZERO(WORK(KT1AMP0),NT1AMX) 1266 LHTF = .FALSE. 1267 CALL CCSD_IAJB(WORK(KT2AMP0),WORK(KT1AMP0),LHTF, 1268 & .FALSE.,.FALSE.,WORK(KEND1),LWRK1) 1269 REWIND(LUIAJB) 1270 CALL WRITI(LUIAJB,IRAT*NT2AMX,WORK(KT2AMP0)) 1271C 1272 IOPT = 3 1273 CALL CC_RDRSP('R0 ',0,1,IOPT,MODEL,WORK(KT1AMP0),WORK(KT2AMP0)) 1274C 1275 IF (CC2) THEN 1276 NONHF = NONHFSAVE 1277 END IF 1278C 1279 RSPIM = .TRUE. 1280 CALL CCRHSN(WORK(KRHO1),WORK(KRHO2), 1281 & WORK(KT1AMP0),WORK(KT2AMP0), 1282 & WORK(KEND1),LWRK1,'XXX') 1283C 1284 IF (IPRINT.GT.10) THEN 1285 CALL AROUND(' END OF CC_FDETA ') 1286 END IF 1287C 1288 RETURN 1289 END 1290*=====================================================================* 1291