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_FTSTNEW(WORK,LWORK,APROXR12) 21 IMPLICIT NONE 22#include "priunit.h" 23#include "ccorb.h" 24#include "ccsdsym.h" 25#include "ccsdinp.h" 26#include "ccfro.h" 27#include "dummy.h" 28 29 CHARACTER*3 APROXR12,LISTL,LISTR 30 CHARACTER*6 FILFMAT 31 CHARACTER*8 CDUMMY 32 CHARACTER*10 MODEL,MODELW 33 34 LOGICAL LOCDBG 35 PARAMETER(LOCDBG = .FALSE.) 36 INTEGER NFTRAN,IFTRAN(3,1),IOPTRES,KRHOANA1,KRHOANA2, 37 & KEND1,KRHOANA12,KRHODIF1,KRHODIF2,KRHODIF12,LWORK, 38 & LWRK1,IDLSTL,IDLSTR,ISYM,IR1TAMP,LENALL,LUFMAT, 39 & iopt,kt0amp0,kt12,lenmod,ISYCTR,ISYAMP,ISYRES, 40 & isymkl,isymij 41 42#if defined (SYS_CRAY) 43 REAL WORK(LWORK), ERRNRM, DDOT 44#else 45 DOUBLE PRECISION WORK(LWORK), ERRNRM, DDOT 46#endif 47 48 !external functions: 49 INTEGER ILSTSYM 50 51 WRITE(LUPRI,*) 'Entered CC_FTSTNEW!' 52 53 ISYM = 1 54 55 LISTL = 'L0 ' 56 IDLSTL = 0 57 58 LISTR = 'R1 ' 59 IDLSTR = IR1TAMP('ZDIPLEN ',.FALSE.,0.0D0,ISYM) 60 NFTRAN = 1 61 62 IFTRAN(1,1) = IDLSTL 63 IFTRAN(2,1) = IDLSTR 64 65 IOPTRES = 0 66 FILFMAT = 'CC_FMAT' 67 68C determine symmetry of result vector: 69 ISYCTR = ILSTSYM(LISTL,IDLSTL) 70 ISYAMP = ILSTSYM(LISTR,IDLSTR) 71 ISYRES = MULD2H(ISYAMP,ISYCTR) 72 if (locdbg) then 73 write(lupri,*) 'ISYCTR, ISYAMP, ISYRES:', isyctr, isyamp, isyres 74 end if 75 76C Models: 77 if (CCS) then 78 modelw = 'CCS ' 79 else if (CC2) then 80 modelw = 'CC2 ' 81 else if (CCSD) then 82 modelw = 'CCSD ' 83 else 84 call quit('unknown CC method in CC_FTSTNEW!') 85 end if 86 CALL CCSD_MODEL(MODELW,LENMOD,24,MODELW,10,APROXR12) 87C write(lupri,*) 'MODEL in CC_FTSTNEW: ',MODELW 88 89C Generate R12-parts for Lagrange-multipliers and trial-vector 90C from ground state amplitudes (which are total-symmetric) 91C (just to have some values...) 92C if (CCR12) then 93C memory allocation 94C kt0amp0 = 1 95C kt12 = kt0amp0 + 2*nallai(1) 96C kend1 = kt12 + ngamma(1) 97C lwrk1 = lwork - kend1 98C if (lwrk1 .lt. 0) then 99C call quit('Insufficient memory in CC_FTSTNEW for R12!') 100C end if 101C 102C check symmetry: 103C if (ISYCTR .NE. 1) then 104C call quit('Symmetry mismatch while generating 105C & Lagrange-multipliers') 106C else if (ISYAMP .NE. 1) then 107C call quit('Symmetry mismatch while generating response 108C & amplitudes') 109C end if 110C 111C iopt=32 112C call cc_rdrsp('R0 ',0,1,iopt,model,dummy,work(kt12)) 113C 114C if (locdbg) then 115C write(lupri,*) 'Norm^2(Ground state R12-ampl.)', 116C & ddot(ngamma(1),work(kt12),1,work(kt12),1) 117C call cc_prpr12(work(kt12),1,1,.TRUE.) 118C end if 119C 120C 121C do i=1, ngamma(1) 122C work(kt12 + i - 1) = 100.0d0 * 123C & work(kt12 + i - 1)*dabs(work(kt12+i-1)) 124C end do 125C 126C if (locdbg) then 127C write(lupri,*) 'Norm^2(R12-Zeta)', 128C & ddot(ngamma(1),work(kt12),1,work(kt12),1) 129C call cc_prpr12(work(kt12),1,1,.TRUE.) 130C end if 131C 132C iopt=32 133C call cc_wrrsp(listl,idlstl,1,iopt,modelw,work(kt0amp0),dummy, 134C & work(kt12),work(kend1),lwrk1) 135C do i=1, ngamma(1) 136C work(kt12 + i - 1) = 100.0d0 * 137C & work(kt12 + i - 1)*dabs(work(kt12+i-1)) 138C end do 139C 140C if (locdbg) then 141C write(lupri,*) 'Norm^2(R12 response ampl.)', 142C & ddot(ngamma(1),work(kt12),1,work(kt12),1) 143C call cc_prpr12(work(kt12),1,1,.TRUE.) 144C end if 145C 146C 147C iopt=32 148C call cc_wrrsp(listr,idlstr,1,iopt,modelw,work(kt0amp0),dummy, 149C & work(kt12),work(kend1),lwrk1) 150C end if 151CCN 152C !TEST: Zero ALL singles amplit., resp. vecs. and lagr. mult.: 153C CALL DZERO(WORK,NT1AM(1)) 154C IOPT = 1 155C CALL CC_WRRSP('R0 ',0,1,IOPT,MODELW,WORK,WORK,DUMMY, 156C & WORK(1+NT1AM(1)),LWORK-NT1AM(1)-1) 157C IOPT = 1 158C CALL CC_WRRSP(LISTR,IDLSTR,1,IOPT,MODELW,WORK,WORK,DUMMY, 159C & WORK(1+NT1AM(1)),LWORK-NT1AM(1)-1) 160C IOPT = 1 161C CALL CC_WRRSP(LISTL,IDLSTL,1,IOPT,MODELW,WORK,WORK,DUMMY, 162C & WORK(1+NT1AM(1)),LWORK-NT1AM(1)-1) 163CCN 164 CALL CC_FMATRIX(IFTRAN,NFTRAN,LISTL,LISTR,IOPTRES,FILFMAT, 165 & IDUMMY,DUMMY,0,WORK,LWORK) 166 167 KRHOANA1 = IFTRAN(3,1) 168 KRHOANA2 = KRHOANA1 + NT1AM(ISYRES) 169 KEND1 = KRHOANA2 + NT2AM(ISYRES) 170 IF (CCR12) THEN 171 KRHOANA12 = KEND1 172 KEND1 = KRHOANA12 + NTR12AM(ISYRES) 173 END IF 174 LWRK1 = LWORK - KEND1 175 IF (LWRK1.LT.0) CALL QUIT('Insufficient memory in CC_FTSTNEW') 176 177 LENALL = NT1AM(ISYRES) + NT2AM(ISYRES) 178 IF (CCS) LENALL = NT1AM(ISYRES) 179 IF (CCR12) LENALL = LENALL + NTR12AM(ISYRES) 180 181 CALL DZERO(WORK(KRHOANA1),LENALL) 182 183 LUFMAT = -1 184 CALL WOPEN2(LUFMAT,FILFMAT,64,0) 185 CALL GETWA2(LUFMAT,FILFMAT,WORK(KRHOANA1),1,LENALL) 186 CALL WCLOSE2(LUFMAT,FILFMAT,'DELETE') 187 188 IF (NSYM .EQ. 1) then 189 KRHODIF1 = KEND1 190 KRHODIF2 = KRHODIF1 + NT1AMX 191 KEND1 = KRHODIF2 + NT2AMX 192 IF (CCR12) THEN 193 KRHODIF12 = KEND1 194 KEND1 = KRHODIF12 + NTR12AM(1) 195 END IF 196 197 LWRK1 = LWORK - KEND1 198 IF (LWRK1.LT.0) CALL QUIT('Insufficient memory in CC_FTSTNEW') 199 CALL DZERO(WORK(KRHODIF1),LENALL) 200 201 CALL CC_FDF(WORK(KRHODIF1),WORK(KRHODIF2),WORK(KRHODIF12), 202 & LISTL,IDLSTL,LISTR,IDLSTR, 203 & WORK(KEND1),LWRK1,APROXR12) 204 ELSE 205 write (lupri,*) 'CC_FDF must run in C1 symmetry!' 206 write (lupri,*) 'Will ignore finite differences....' 207 END IF 208 209 N = 1 210 IF (CCS) N = 0 211 212 IF (NSYM .EQ. 1) THEN 213 WRITE(LUPRI,'(/A)') 'Finite difference Result for FMATRIX:' 214 CALL CC_PRP(WORK(KRHODIF1),WORK(KRHODIF2),1,1,N) 215 IF (CCR12) THEN 216 CALL CC_PRPR12(WORK(KRHODIF12),1,1,.TRUE.) 217 END IF 218 END IF 219 WRITE(LUPRI,'(//A)') 'Analytical Result for FMATRIX:' 220 CALL CC_PRP(WORK(KRHOANA1),WORK(KRHOANA2),ISYRES,1,N) 221 IF (CCR12) THEN 222 CALL CC_PRPR12(WORK(KRHOANA12),ISYRES,1,.TRUE.) 223 WRITE(LUPRI,*) 'Norm^2 of FMATRIX: ', 224 & DDOT(NT1AM(ISYRES),WORK(KRHOANA1),1,WORK(KRHOANA1),1) + 225 & DDOT(NT2AM(ISYRES),WORK(KRHOANA2),1,WORK(KRHOANA2),1) + 226 & DDOT(NTR12AM(ISYRES),WORK(KRHOANA12),1,WORK(KRHOANA12),1) 227 END IF 228 229 IF (NSYM .EQ. 1) THEN 230 CALL DAXPY(NT1AMX,-1.0D0,WORK(KRHODIF1),1,WORK(KRHOANA1),1) 231 CALL DAXPY(NT2AMX,-1.0D0,WORK(KRHODIF2),1,WORK(KRHOANA2),1) 232 ERRNRM = DDOT(NT1AMX,WORK(KRHOANA1),1,WORK(KRHOANA1),1) + 233 & DDOT(NT2AMX,WORK(KRHOANA2),1,WORK(KRHOANA2),1) 234 IF (CCR12) THEN 235 CALL DAXPY(NTR12AM(1),-1.0D0,WORK(KRHODIF12),1, 236 & WORK(KRHOANA12),1) 237 ERRNRM = ERRNRM + 238 & DDOT(NTR12AM(1),WORK(KRHOANA12),1,WORK(KRHOANA12),1) 239 END IF 240 241 WRITE(LUPRI,'(//A)') 'Norm of Difference between Results:' 242 WRITE(LUPRI,*) 'Norm:',ERRNRM 243 WRITE(LUPRI,*) 'Singles part: ', 244 & DDOT(NT1AMX,WORK(KRHOANA1),1,WORK(KRHOANA1),1) 245 WRITE(LUPRI,*) 'Doubles part: ', 246 & DDOT(NT2AMX,WORK(KRHOANA2),1,WORK(KRHOANA2),1) 247 IF (CCR12) WRITE(LUPRI,*) 'R12-Doubles part: ', 248 & DDOT(NTR12AM(1),WORK(KRHOANA12),1,WORK(KRHOANA12),1) 249 CALL CC_PRP(WORK(KRHOANA1),WORK(KRHOANA2),1,1,N) 250 IF (CCR12) THEN 251 CALL CC_PRPR12(WORK(KRHOANA12),1,1,.TRUE.) 252 END IF 253 END IF 254 255 WRITE(LUPRI,*) 'LEAVING CC_FTSTNEW!' 256 RETURN 257 END 258*=====================================================================* 259 SUBROUTINE CC_FDF(RHODIF1,RHODIF2,RHODIF12, 260 & LISTL,IDLSTL,LISTR,IDLSTR, 261 & WORK,LWORK,APROXR12) 262C 263C---------------------------------------------------------------------* 264C Christian Neiss, Christof Haettig, october 2004 265C 266C test routine for F-matrix. calculates contractions of the 267C F-matrix by numerical differentiation of the Jacobian. 268C---------------------------------------------------------------------* 269C 270 IMPLICIT NONE 271#include "priunit.h" 272#include "ccsdsym.h" 273#include "ccorb.h" 274#include "ccsdinp.h" 275#include "ccfro.h" 276#include "dummy.h" 277#include "r12int.h" 278#include "ccr12int.h" 279 280* local parameters: 281 LOGICAL LOCDBG, LTESTV 282 PARAMETER (LOCDBG = .FALSE.) 283 284 INTEGER LWORK 285#if defined (SYS_CRAY) 286 REAL WORK(LWORK) 287 REAL RHODIF1(*), RHODIF2(*), RHODIF12(*) 288 REAL ECURR, FAC, DDOT 289 REAL HALF, ZERO, ONE, TWO, DELTA, DELINV 290#else 291 DOUBLE PRECISION WORK(LWORK) 292 DOUBLE PRECISION RHODIF1(*), RHODIF2(*), RHODIF12(*) 293 DOUBLE PRECISION ECURR, FAC, DDOT 294 DOUBLE PRECISION HALF, ZERO, ONE, TWO, DELTA, DELINV 295#endif 296 297 CHARACTER*(3) LISTR, LISTL, APROXR12 298 CHARACTER*(10) MODEL, MODELW 299 INTEGER IDLSTL, IDLSTR, ISYM, IOPT, ISIDE 300 INTEGER KT1AMPSAV, KT2AMPSAV, KT12AMPSAV, KT1AMP0, KT2AMP0, 301 & KT12AMP0, KT1AMPA, KT2AMPA, KT12AMPA, KOMEGA1, KOMEGA2 302 INTEGER KZETA1, KZETA2, KZETA12, KT0AMP0, ISYM0, ISYMR, ISYML 303 INTEGER KRHO1, KRHO2, KRHO12, KETA1, KETA2, KETA12, LUNIT 304 INTEGER IDUM, LDUM, LENMOD 305 INTEGER KEND1, LWRK1, KEND2, LWRK2 306 INTEGER KLAMDP, KLAMDH, KVABKL, KVDIFF, NVCDKL, KVCDKL 307 308 ! external functions: 309 INTEGER IR1TAMP 310 311 312 PARAMETER (HALF=0.5D0, ZERO=0.0D0, ONE=1.0D0, TWO=2.0D0) 313 314C---------------------------------------------------------------------* 315C displacement parameter for finite difference: 316C---------------------------------------------------------------------* 317 PARAMETER (DELTA = 1.0D-07, DELINV = 1.0D+07) 318C PARAMETER (DELTA = 1.0D-02, DELINV = 1.0D+02) 319 320 321 CALL QENTER('CC_FDF') 322 323C---------------------------------------------------------------------* 324C Initialisation 325C---------------------------------------------------------------------* 326 327 LTESTV = .FALSE..AND.CCR12.AND..NOT.(CCS.OR.CC2) 328 329C calc. dimension: 330 NVCDKL = 0 331 DO I = 1, NSYM 332 NVCDKL = NVCDKL + NVIR(1)*NRHFB(1)*(NVIR(1)*NRHFB(1)+1)/2 333 END DO 334 335C only C1-symmetry: 336 if (nsym.gt.1) then 337 write (lupri,*) 'CC_FDF must run in C1 symmetry!' 338 write (lupri,*) 'Will ignore finite differences....' 339 goto 1000 340 end if 341 isym = 1 342 isym0 = 1 343 isyml = 1 344 isymr = 1 345 346C Models: 347 modelw = 'UNKNOWN ' 348 if (CCS) modelw = 'CCS ' 349 if (CC2) modelw = 'CC2 ' 350 if (CCSD) modelw = 'CCSD ' 351 352 CALL CCSD_MODEL(MODELW,LENMOD,24,MODELW,10,APROXR12) 353 354C memory allocation: 355 kt1ampsav = 1 356 kt2ampsav = kt1ampsav + nt1amx 357 kend1 = kt2ampsav + nt2amx 358 if (CCR12) then 359 kt12ampsav = kend1 360 kend1 = kt12ampsav + ntr12am(isym) 361 end if 362 363 kt1ampa = kend1 364 kt2ampa = kt1ampa + nt1amx 365 kend1 = kt2ampa + nt2amx 366 if (CCR12) then 367 kt12ampa = kend1 368 kend1 = kt12ampa + ntr12am(isym) 369C dirty hack; otherwise problems with g77-compiler! 370 else 371 kt12ampa = lwork-1 372 end if 373 374 keta1 = kend1 375 keta2 = keta1 + nt1amx 376 kend1 = keta2 + nt2amx 377 if (CCR12) then 378 keta12 = kend1 379 kend1 = keta12 + ntr12am(1) 380 end if 381 382 kzeta1 = kend1 383 kzeta2 = kzeta1 + nt1amx 384 kend1 = kzeta2 + nt2amx 385 if (CCR12) then 386 kzeta12 = kend1 387 kend1 = kzeta12 + ntr12am(isym) 388 end if 389 390 kt0amp0 = kend1 391 kt1amp0 = kt0amp0 + 2*nallai(isym) 392 komega1 = kt1amp0 + nt1amx 393 komega2 = komega1 + nt1amx 394 kt2amp0 = komega2 + max(nt2amx,2*nt2ort(1),nt2ao(1)) 395 kend1 = kt2amp0 + max(nt2sq(1),nt2r12(1)) 396 if (CCR12) then 397 kt12amp0 = kend1 398 kend1 = kt12amp0 + ntr12am(isym) 399 end if 400 401 if (ltestv) then 402 kvdiff = kend1 403 kvcdkl = kvdiff + nvcdkl 404 kend1 = kvcdkl + nvcdkl 405 end if 406 407 krho1 = kend1 408 krho2 = krho1 + nt1amx 409 kend1 = krho2 + nt2amx 410 if (CCR12) then 411 krho12 = kend1 412 kend1 = krho12 + ntr12am(isym) 413 end if 414 lwrk1 = lwork - kend1 415 416C test if work space is enough: 417 if (lwrk1 .lt. 0) then 418 write(lupri,*) 'Not enough work memory in cc_FDF:' 419 write(lupri,*) 'Available: LWORK = ',lwork 420 write(lupri,*) 'Needed: ', kend1 421 call quit('Too little work space in CC_FDF!') 422 end if 423 424C zero out result arrays: 425 call dzero(rhodif1,nt1amx) 426 call dzero(rhodif2,nt2amx) 427 if (CCR12) call dzero(rhodif12,ntr12am(isym)) 428 429 if (ltestv) then 430 call dzero(work(kvcdkl),nvcdkl) 431 call dzero(work(kvdiff),nvcdkl) 432 end if 433 434C---------------------------------------------------------------------* 435C Read in Lagrange multipliers 436C---------------------------------------------------------------------* 437 iopt=3 438 call cc_rdrsp(listl,idlstl,isyml,iopt,model, 439 & work(kzeta1),work(kzeta2)) 440 if (locdbg) then 441 write(lupri,*) 'norm of zeta1 at 0:', 442 & ddot(nt1amx,work(kzeta1),1,work(kzeta1),1) 443 end if 444 if (CCR12) then 445 iopt=32 446 call cc_rdrsp(listl,idlstl,isyml,iopt,model, 447 & dummy,work(kzeta12)) 448 end if 449 450C---------------------------------------------------------------------* 451C Read in cluster amplitudes 452C---------------------------------------------------------------------* 453 iopt=3 454 call cc_rdrsp('R0 ',0,isym0,iopt,model,work(kt1ampsav), 455 & work(kt2ampsav)) 456 if (CCR12) then 457 iopt=32 458 call cc_rdrsp('R0 ',0,isym0,iopt,model, 459 & dummy,work(kt12ampsav)) 460 end if 461 462C---------------------------------------------------------------------* 463C Read in trial vector 464C---------------------------------------------------------------------* 465 iopt=3 466 call cc_rdrsp(listr,idlstr,isymr,iopt,model,work(kt1ampa), 467 & work(kt2ampa)) 468C trial vector diagonal elements must be scaled by 2 due to 469C different format than ground state amplitudes: 470C (CCRHSN expects ground state format) 471 call cclr_diascl(work(kt2ampa),two,1) 472 if (CCR12) then 473 iopt=32 474 call cc_rdrsp(listr,idlstr,isymr,iopt,model,dummy, 475 & work(kt12ampa)) 476C trial vector diagonal elements must be scaled by KETSCL due to 477C different format than ground state amplitudes: 478C (CCRHSN expects ground state format) 479 call cclr_diasclr12(work(kt12ampa),ketscl,1) 480 if (locdbg) then 481 write(lupri,*) 'T12AMPA in CC_FDF:' 482 call outpak(work(kt12ampa),nmatki(1),1,lupri) 483 end if 484 end if 485 if (locdbg) then 486 write(lupri,*) 'norm of zeta1 at 1:', 487 & ddot(nt1amx,work(kzeta1),1,work(kzeta1),1) 488 end if 489 490C---------------------------------------------------------------------* 491C Loop for +/- displacement DELTA 492C---------------------------------------------------------------------* 493 do n = 1, 2 494 495C---------------------------------------------------------------------* 496C add/subtract finite displacement Delta*Trialvector to t 497C and recalculate intermediates 498C---------------------------------------------------------------------* 499 call dcopy(nt1amx,work(kt1ampsav),1,work(kt1amp0),1) 500 call dcopy(nt2amx,work(kt2ampsav),1,work(kt2amp0),1) 501 if (CCR12) then 502 call dcopy(ntr12am(isym),work(kt12ampsav),1,work(kt12amp0),1) 503 end if 504 505 if (n.eq.1) fac=one 506 if (n.eq.2) fac=-one 507 508 call daxpy(nt1amx,fac*delta,work(kt1ampa),1,work(kt1amp0),1) 509 call daxpy(nt2amx,fac*delta,work(kt2ampa),1,work(kt2amp0),1) 510 if (CCR12) then 511 call daxpy(ntr12am(isym),fac*delta,work(kt12ampa),1, 512 & work(kt12amp0),1) 513 end if 514 515 iopt=3 516 call cc_wrrsp('R0 ',0,isym0,iopt,modelw,work(kt0amp0), 517 & work(kt1amp0),work(kt2amp0),work(kend1),lwrk1) 518 if (CCR12) then 519 iopt = 32 520 call cc_wrrsp('R0 ',0,isym0,iopt,modelw,work(kt0amp0), 521 & dummy,work(kt12amp0),work(kend1),lwrk1) 522 end if 523 524 if (CCR12.AND..NOT.(CCS.OR.CC2)) then 525 klamdp = kend1 526 klamdh = klamdp + nlamdt 527 kvabkl = klamdh + nlamdt 528 kend2 = kvabkl + nvabkl(1) 529 lwrk2 = lwork - kend2 530 if (lwrk2.lt.0) then 531 call quit('Insufficient work space for Vabkl in CC_FDF') 532 end if 533 534 call lammat(work(klamdp),work(klamdh),work(kt1amp0), 535 & work(kend2),lwrk2) 536 537 lunit = -1 538 call gpopen(lunit,fvabkl,'OLD',' ','UNFORMATTED',idummy, 539 & .FALSE.) 540 read(lunit) (work(kvabkl+i-1), i=1, nvabkl(1)) 541 call gpclose(lunit,'KEEP') 542 543 iopt = 1 544 call cc_r12mkvirt(work(kvabkl),work(klamdp),1,work(klamdp),1, 545 & 'R12VCTDTKL',iopt,work(kend2),lwrk2) 546 end if 547 548 rspim=.TRUE. 549 call ccrhsn(work(komega1),work(komega2), 550 & work(kt1amp0),work(kt2amp0), 551 & work(kend1),lwrk1,aproxr12) 552 553C---------------------------------------------------------------------* 554C Do left hand transformation with Lagrange multipliers 555C---------------------------------------------------------------------* 556 if (locdbg) then 557 write(lupri,*) 'norm^2 of zeta1 at 2:', 558 & ddot(nt1amx,work(kzeta1),1,work(kzeta1),1) 559 end if 560 call dcopy(nt1amx,work(kzeta1),1,work(krho1),1) 561 call dcopy(nt2amx,work(kzeta2),1,work(krho2),1) 562 if (CCR12) then 563 call dcopy(ntr12am(isym),work(kzeta12),1,work(krho12),1) 564 end if 565 566 ecurr=0.0D0 567 iside=-1 568 569 call cc_atrr(ecurr,isym,iside,work(krho1),lwrk1,.FALSE.,dummy, 570 & aproxr12,.FALSE.) 571 572 call cc_eta(work(keta1),work(kend1),lwrk1) 573 574C we do not need to add doubles- and r12-contributions of eta 575C to rho, since these quantities are independent on amplitudes: 576 call daxpy(nt1amx,one,work(keta1),1,work(krho1),1) 577 578 if (locdbg) then 579 write(lupri,*) 'norm of zeta1 at 3:', 580 & ddot(nt1amx,work(kzeta1),1,work(kzeta1),1) 581 end if 582C---------------------------------------------------------------------* 583C divide by 2*delta 584C---------------------------------------------------------------------* 585 call daxpy(nt1amx,fac*half*delinv,work(krho1),1,rhodif1,1) 586 call daxpy(nt2amx,fac*half*delinv,work(krho2),1,rhodif2,1) 587 if (CCR12) then 588 call daxpy(ntr12am(isym),fac*half*delinv,work(krho12),1, 589 & rhodif12,1) 590 end if 591 592 if (ltestv) then 593 lunit = -1 594 call gpopen(lunit,'R12VCTDTKL','OLD',' ','UNFORMATTED',idummy, 595 & .FALSE.) 596 read(lunit) (work(kvcdkl+i-1), i=1, nvcdkl) 597 call gpclose(lunit,'KEEP') 598 599 call daxpy(nvcdkl,fac*half*delinv,work(kvcdkl),1, 600 & work(kvdiff),1) 601 602 if (n.eq.2) then 603 call around('Finite Diff. Result of Vabkl in CC_FDF') 604 call outpkb(work(kvdiff),nvir(1)*nrhfb(1),1,1,lupri) 605 end if 606 end if 607 608C---------------------------------------------------------------------* 609C End of Loop over +/- DELTA 610C---------------------------------------------------------------------* 611 end do 612 613C---------------------------------------------------------------------* 614C restore original amplitudes 615C---------------------------------------------------------------------* 616 call dcopy(nt1amx,work(kt1ampsav),1,work(kt1amp0),1) 617 call dcopy(nt2amx,work(kt2ampsav),1,work(kt2amp0),1) 618 if (CCR12) then 619 call dcopy(ntr12am(isym),work(kt12ampsav),1,work(kt12amp0),1) 620 end if 621 622 iopt=3 623 call cc_wrrsp('R0 ',0,isym0,iopt,modelw,work(kt0amp0), 624 & work(kt1amp0),work(kt2amp0),work(kend1),lwrk1) 625 if (CCR12) then 626 iopt = 32 627 call cc_wrrsp('R0 ',0,isym0,iopt,modelw,work(kt0amp0), 628 & dummy,work(kt12amp0),work(kend1),lwrk1) 629 end if 630 631C---------------------------------------------------------------------* 632C restore original intermediates 633C---------------------------------------------------------------------* 634 if (CCR12.AND..NOT.(CCS.OR.CC2)) then 635 klamdp = kend1 636 klamdh = klamdp + nlamdt 637 kvabkl = klamdh + nlamdt 638 kend2 = kvabkl + nvabkl(1) 639 lwrk2 = lwork - kend2 640 if (lwrk2.lt.0) then 641 call quit('Insufficient work space for Vabkl in CC_FDF') 642 end if 643 644 call lammat(work(klamdp),work(klamdh),work(kt1amp0), 645 & work(kend2),lwrk2) 646 647 lunit = -1 648 call gpopen(lunit,fvabkl,'OLD',' ','UNFORMATTED',idummy, 649 & .FALSE.) 650 read(lunit) (work(kvabkl+i-1), i=1, nvabkl(1)) 651 call gpclose(lunit,'KEEP') 652 653 iopt = 1 654 call cc_r12mkvirt(work(kvabkl),work(klamdp),1,work(klamdp),1, 655 & 'R12VCTDTKL',iopt,work(kend2),lwrk2) 656 end if 657 658 rspim=.TRUE. 659 call ccrhsn(work(komega1),work(komega2), 660 & work(kt1amp0),work(kt2amp0), 661 & work(kend1),lwrk1,aproxr12) 662 6631000 continue 664 CALL QEXIT('CC_FDF') 665 666 RETURN 667 END 668*=====================================================================* 669