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_FDD(NC1VEC,NC2VEC,LISTA,ITAMPA,LISTB,ITAMPB, 21 & LISTC,ITAMPC,TYAM,RESULT,WORK,LWORK) 22C 23C--------------------------------------------------------------------- 24C Test routine for calculating the CC D matrix by finite difference on 25C the C matrix transformation. 26C Ch. Haettig, maj 1997, based on Oves CCLR_FDF routine 27C--------------------------------------------------------------------- 28C 29#include "implicit.h" 30#include "priunit.h" 31#include "dummy.h" 32#include "maxorb.h" 33#include "iratdef.h" 34#include "ccorb.h" 35#include "aovec.h" 36#include "ccsdinp.h" 37#include "cclr.h" 38#include "ccsdsym.h" 39#include "ccsdio.h" 40#include "leinf.h" 41C 42 DIMENSION WORK(LWORK),ITADR(2),RESULT(*) 43 PARAMETER (XHALF = 0.5D00,XMTWO = -2.0D00, DELTA = 1.0D-07) 44 PARAMETER (ONE = 1.0d0, ZERO = 0.0d0, TWO = 2.0d0) 45 CHARACTER MODEL*10 46 CHARACTER FILCMA*8 47 LOGICAL L1TST,L2TST, LETST 48 INTEGER ICTRAN(4) 49C 50 INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 51C 52 MODEL = 'CCSD ' 53 IF (CCS) MODEL = 'CCS ' 54 IF (CC2) MODEL = 'CC2 ' 55C 56 IF (CCR12) CALL QUIT('Finite-difference D-matrix for CCR12 '// 57 & 'not adapted') 58C 59 IF (IPRINT.GT.5) THEN 60 CALL AROUND( 'IN CC_FDD : MAKING FINITE DIFF. CC D Matrix') 61 ENDIF 62C 63C---------------------------- 64C Work space allocations. 65C---------------------------- 66C 67 ISYMTR = 1 68 ISYMOP = 1 69C 70 NTAMP = NT1AM(ISYMTR) + NT2AM(ISYMTR) 71 NTAMP2 = NTAMP*(NC1VEC + NC2VEC ) 72 KF = 1 73 KRHO1 = KF + NTAMP2 74 KRHO2 = KRHO1 + NT1AMX 75 KC1AM = KRHO2 + MAX(NT2AMX,NT2AM(ISYMTR)) 76 KC2AM = KC1AM + NT1AM(ISYMTR) 77 KEND1 = KC2AM 78 * + MAX(NT2AMX,NT2AM(ISYMTR),NT2AO(ISYMTR), 79 * 2*NT2ORT(ISYMTR)) 80 LWRK1 = LWORK - KEND1 81C 82 KRHO1D = KEND1 83 KRHO2D = KRHO1D + NT1AMX 84 KEND2 = KRHO2D 85 * + MAX(NT2AMX,NT2AM(ISYMTR),NT2AO(ISYMTR), 86 * 2*NT2ORT(ISYMTR)) 87 LWRK2 = LWORK - KEND1 88C 89 IF (IPRINT .GT. 100 ) THEN 90 WRITE(LUPRI,*) ' IN CC_FDD: KF = ',KF 91 WRITE(LUPRI,*) ' IN CC_FDD: KRHO1 = ',KRHO1 92 WRITE(LUPRI,*) ' IN CC_FDD: KRHO2 = ',KRHO2 93 WRITE(LUPRI,*) ' IN CC_FDD: KC1AM = ',KC1AM 94 WRITE(LUPRI,*) ' IN CC_FDD: KC2AM = ',KC2AM 95 WRITE(LUPRI,*) ' IN CC_FDD: KRHO1D = ',KRHO1D 96 WRITE(LUPRI,*) ' IN CC_FDD: KRHO2D = ',KRHO2D 97 WRITE(LUPRI,*) ' IN CC_FDD: KEND2 = ',KEND2 98 WRITE(LUPRI,*) ' IN CC_FDD: LWRK2 = ',LWRK2 99 ENDIF 100 IF (LWRK2.LT.0 ) THEN 101 WRITE(LUPRI,*) 'Too little work space in CC_FDD ' 102 WRITE(LUPRI,*) 'AVAILABLE: LWORK = ',LWORK 103 WRITE(LUPRI,*) 'NEEDED (AT LEAST) = ',KEND2 104 CALL QUIT('TOO LITTLE WORKSPACE IN CC_FDD ') 105 ENDIF 106 KF2 = KF + NC1VEC*NTAMP 107C 108C--------------------- 109C Initializations. 110C--------------------- 111C 112 CALL DZERO(WORK(KC1AM),NT1AMX) 113 CALL DZERO(WORK(KC2AM),NT2AMX) 114 CALL DZERO(WORK(KF),NTAMP2) 115 IF (ABS(DELTA) .GT. 1.0D-15 ) THEN 116 DELTAI = 1.0D00/DELTA 117 ELSE 118 DELTAI = 1 119 ENDIF 120 X11 = 0.0D00 121 X12 = 0.0D00 122 X21 = 0.0D00 123 X22 = 0.0D00 124 XNJ = 0.0D00 125C 126C------------------------------------------------ 127C Read the CC reference amplitudes From disk. 128C------------------------------------------------ 129C 130 IOPT = 3 131 CALL CC_RDRSP('R0',0,1,IOPT,MODEL,WORK(KC1AM),WORK(KC2AM)) 132C 133C---------------------------------------------- 134C Save the CC reference amplitudes on disk. 135C---------------------------------------------- 136C 137 LUTAM = -1 138 CALL GPOPEN(LUTAM,'TAM_SAV','UNKNOWN',' ','UNFORMATTED',IDUMMY, 139 * .FALSE.) 140 REWIND(LUTAM) 141 WRITE(LUTAM) (WORK(KC1AM + I -1 ), I = 1, NT1AMX) 142 WRITE(LUTAM) (WORK(KC2AM + I -1 ), I = 1, NT2AMX) 143 CALL GPCLOSE(LUTAM,'KEEP') 144C 145 IF (IPRINT.GT.125) THEN 146 RHO1N = DDOT(NT1AMX,WORK(KC1AM),1,WORK(KC1AM),1) 147 RHO2N = DDOT(NT2AMX,WORK(KC2AM),1,WORK(KC2AM),1) 148 WRITE(LUPRI,*) 'Norm of T1AM: ',RHO1N 149 WRITE(LUPRI,*) 'Norm of T2AM: ',RHO2N 150 CALL CC_PRP(WORK(KC1AM),WORK(KC2AM),1,1,1) 151 ENDIF 152 RSPIM = .TRUE. 153C 154C------------------------------------------ 155C Calculate reference A*T vector. 156C------------------------------------------ 157C 158 ICTRAN(1) = ITAMPA 159 ICTRAN(2) = ITAMPB 160 ICTRAN(3) = ITAMPC 161 ICTRAN(4) = 0 162 163 NCTRAN = 1 164 IOPT = 1 165 FILCMA = 'CCCMAT' 166 167 CALL CC_CMAT(ICTRAN,NCTRAN,LISTA,LISTB,LISTC,IOPT, 168 & FILCMA,IDUM,RDUM,0,WORK(KRHO1D),LWORK-KRHO1D) 169 170C 171C------------------------- 172C Zero out components. 173C------------------------- 174C 175 IF (LCOR .OR. LSEC) THEN 176C 177 CALL CC_CORE(WORK(KRHO1D),WORK(KRHO2D),ISYMTR) 178C 179 ENDIF 180C 181 IF (IPRINT.GT.2) THEN 182 RHO1N = DDOT(NT1AMX,WORK(KRHO1D),1,WORK(KRHO1D),1) 183 RHO2N = DDOT(NT2AMX,WORK(KRHO2D),1,WORK(KRHO2D),1) 184 WRITE(LUPRI,*) 'Norm of RHO1: ',RHO1N,'ref' 185 WRITE(LUPRI,*) 'Norm of RHO2: ',RHO2N,'ref' 186 ENDIF 187 IF (IPRINT.GT.125) THEN 188 CALL CC_PRP(WORK(KRHO1D),WORK(KRHO2D),1,1,1) 189 ENDIF 190 191 CALL DCOPY(NT1AMX,WORK(KRHO1D),1,WORK(KRHO1),1) 192 CALL DCOPY(NT2AMX,WORK(KRHO2D),1,WORK(KRHO2),1) 193C 194C============================================= 195C Calculate C-matrix by finite difference. 196C============================================= 197C 198 DO 100 I = 1, NC1VEC 199 WRITE (LUPRI,*) 'singles index:',I 200C 201C---------------------------------------- 202C Add finite displadement to t and 203C calculate new intermediates. 204C---------------------------------------- 205C 206 LUTAM = -1 207 CALL GPOPEN(LUTAM,'TAM_SAV','UNKNOWN',' ','UNFORMATTED',IDUMMY, 208 * .FALSE.) 209 READ(LUTAM) (WORK(KC1AM + J -1 ) , J = 1, NT1AMX) 210 READ(LUTAM) (WORK(KC2AM + J -1 ) , J = 1, NT2AMX) 211 CALL GPCLOSE(LUTAM,'KEEP') 212C 213 TI = SECOND() 214 WORK(KC1AM +I -1) = WORK(KC1AM +I -1 ) + DELTA 215 IF (LCOR .OR. LSEC) THEN 216 CALL CC_CORE(WORK(KC1AM),WORK(KC2AM),ISYMTR) 217 ENDIF 218C 219 IOPT = 3 220 CALL CC_WRRSP('R0',0,1,IOPT,MODEL,DUMMY,WORK(KC1AM), 221 * WORK(KC2AM),WORK(KEND2),LWRK2) 222C 223 RSPIM = .TRUE. 224 CALL CCRHSN(WORK(KRHO1D),WORK(KRHO2D),WORK(KC1AM), 225 * WORK(KC2AM),WORK(KEND2),LWRK2,'XXX') 226C 227C------------------ 228C Transform. 229C------------------ 230C 231 ICTRAN(1) = ITAMPA 232 ICTRAN(2) = ITAMPB 233 ICTRAN(3) = ITAMPC 234 ICTRAN(4) = 0 235 236 NCTRAN = 1 237 IOPT = 1 238 FILCMA = 'CCCMAT' 239 240 CALL CC_CMAT(ICTRAN,NCTRAN,LISTA,LISTB,LISTC,IOPT, 241 & FILCMA,IDUM,RDUM,0,WORK(KRHO1D),LWORK-KRHO1D) 242 243 244 IF (LCOR .OR. LSEC) THEN 245 CALL CC_CORE(WORK(KRHO1D),WORK(KRHO2D),ISYMTR) 246 ENDIF 247C 248 IF (IPRINT.GT.2) THEN 249 RHO1N = DDOT(NT1AMX,WORK(KRHO1D),1,WORK(KRHO1D),1) 250 RHO2N = DDOT(NT2AMX,WORK(KRHO2D),1,WORK(KRHO2D),1) 251 WRITE(LUPRI,*) 'Norm of RHO1: ',RHO1N,'ai=',I 252 WRITE(LUPRI,*) 'Norm of RHO2: ',RHO2N,'ai=',I 253 ENDIF 254 IF (IPRINT.GT.125) THEN 255 CALL CC_PRP(WORK(KRHO1D),WORK(KRHO2D),1,1,1) 256 ENDIF 257 CALL DAXPY(NT1AMX,-1.0D00,WORK(KRHO1),1,WORK(KRHO1D),1) 258 CALL DAXPY(NT2AMX,-1.0D00,WORK(KRHO2),1,WORK(KRHO2D),1) 259 CALL DSCAL(NT1AMX,DELTAI,WORK(KRHO1D),1) 260 CALL DSCAL(NT2AMX,DELTAI,WORK(KRHO2D),1) 261 CALL DCOPY(NT1AMX,WORK(KRHO1D),1, 262 * WORK(KF+NTAMP*(I-1)),1) 263 CALL DCOPY(NT2AMX,WORK(KRHO2D),1, 264 * WORK(KF+NTAMP*(I-1)+NT1AMX),1) 265 X11 = X11 + DDOT(NT1AMX,WORK(KRHO1D),1,WORK(KRHO1D),1) 266 X21 = X21 + DDOT(NT2AMX,WORK(KRHO2D),1,WORK(KRHO2D),1) 267C 268 TI = SECOND() - TI 269 IF (IPRINT.GT.5 ) THEN 270 WRITE(LUPRI,*) ' ' 271 WRITE(LUPRI,*) 'FDB ROW NR. ',I,' DONE IN ',TI,' SEC.' 272 ENDIF 273C 274 100 CONTINUE 275C 276C---------------------------------------------------------------- 277C Loop over T2 amplitudes. Take care of diagonal t2 elements 278C is in a different convention in the energy code. 279C Factor 1/2 from right , and factor 2 from left. 280C---------------------------------------------------------------- 281C 282 DO 200 NAI = 1, NT1AMX 283 DO 300 NBJ = 1, NAI 284 I = INDEX(NAI,NBJ) 285C 286 IF (I.LE.NC2VEC) THEN 287 WRITE (LUPRI,*) 'doubles index:',I 288C 289C-------------------------------------------- 290C Add finite displacement to t and 291C calculate new intermediates. 292C------------------------------------------- 293C 294 LUTAM = -1 295 CALL GPOPEN(LUTAM,'TAM_SAV','UNKNOWN',' ','UNFORMATTED', 296 * IDUMMY,.FALSE.) 297 READ(LUTAM) (WORK(KC1AM + J -1 ) , J = 1, NT1AMX) 298 READ(LUTAM) (WORK(KC2AM + J -1 ) , J = 1, NT2AMX) 299 CALL GPCLOSE(LUTAM,'KEEP') 300C 301 TI = SECOND() 302 DELT = DELTA 303 IF (NAI.EQ.NBJ) DELT = 2*DELTA 304 WORK(KC2AM + I -1) = WORK(KC2AM+I -1) + DELT 305 IF (LCOR .OR. LSEC) THEN 306 CALL CC_CORE(WORK(KC1AM),WORK(KC2AM),ISYMTR) 307 ENDIF 308C 309 IOPT = 3 310 CALL CC_WRRSP('R0',0,1,IOPT,MODEL,DUMMY,WORK(KC1AM), 311 * WORK(KC2AM),WORK(KEND2),LWRK2) 312C 313 RSPIM = .TRUE. 314 CALL CCRHSN(WORK(KRHO1D),WORK(KRHO2D),WORK(KC1AM), 315 * WORK(KC2AM),WORK(KEND2),LWRK2,'XXX') 316C 317C-------------------- 318C Transform. 319C-------------------- 320C 321 ICTRAN(1) = ITAMPA 322 ICTRAN(2) = ITAMPB 323 ICTRAN(3) = ITAMPC 324 ICTRAN(4) = 0 325 326 NCTRAN = 1 327 IOPT = 1 328 FILCMA = 'CCCMAT' 329 330 CALL CC_CMAT(ICTRAN,NCTRAN,LISTA,LISTB,LISTC,IOPT, 331 & FILCMA,IDUM,RDUM,0,WORK(KRHO1D),LWORK-KRHO1D) 332 333 334 IF (LCOR .OR. LSEC) THEN 335 CALL CC_CORE(WORK(KRHO1D),WORK(KRHO2D),ISYMTR) 336 ENDIF 337C 338 IF (IPRINT.GT.2) THEN 339 RHO1N = DDOT(NT1AMX,WORK(KRHO1D),1,WORK(KRHO1D),1) 340 RHO2N = DDOT(NT2AMX,WORK(KRHO2D),1,WORK(KRHO2D),1) 341 WRITE(LUPRI,*) 'Norm of RHO1: ',RHO1N,'aibj=',I 342 WRITE(LUPRI,*) 'Norm of RHO2: ',RHO2N,'aibj=',I 343 ENDIF 344 IF (IPRINT.GT.125) THEN 345 CALL CC_PRP(WORK(KRHO1D),WORK(KRHO2D),1,1,1) 346 ENDIF 347C 348 CALL DAXPY(NT1AMX,-1.0D00,WORK(KRHO1),1,WORK(KRHO1D),1) 349 CALL DAXPY(NT2AMX,-1.0D00,WORK(KRHO2),1,WORK(KRHO2D),1) 350 CALL DSCAL(NT1AMX,DELTAI,WORK(KRHO1D),1) 351 CALL DSCAL(NT2AMX,DELTAI,WORK(KRHO2D),1) 352 CALL DCOPY(NT1AMX,WORK(KRHO1D),1, 353 * WORK(KF2+NTAMP*(I-1)),1) 354 CALL DCOPY(NT2AMX,WORK(KRHO2D),1, 355 * WORK(KF2+NTAMP*(I-1)+NT1AMX),1) 356C 357 X12 = X12 + DDOT(NT1AMX,WORK(KRHO1D),1,WORK(KRHO1D),1) 358 X22 = X22 + DDOT(NT2AMX,WORK(KRHO2D),1,WORK(KRHO2D),1) 359 TI = SECOND() - TI 360 IF (IPRINT.GT.5 ) THEN 361 WRITE(LUPRI,*) ' ' 362 WRITE(LUPRI,*) 'FDB ROW NR. ',I+NT1AMX, 363 * ' DONE IN ',TI,' SEC.' 364 ENDIF 365C 366 ENDIF 367C 368 300 CONTINUE 369 200 CONTINUE 370C 371 WRITE(LUPRI,*) ' ' 372 WRITE(LUPRI,*) '** FINITE DIFF WITH DELTA ',DELTA, '**' 373 WRITE(LUPRI,*) ' ' 374 IF ((IPRINT .GT. 4).AND.(.TRUE.)) THEN 375 CALL AROUND('FINITE DIFF. CC D*TA*TB*TC-Matrix - 11 & 21 PART' ) 376 CALL OUTPUT(WORK(KF),1,NTAMP,1,NC1VEC,NTAMP,NC1VEC,1,LUPRI) 377 CALL AROUND('FINITE DIFF. CC D*TA*TB*TC-Matrix - 12 & 22 PART' ) 378 CALL OUTPUT(WORK(KF+NTAMP*NC1VEC),1,NTAMP,1,NC2VEC, 379 * NTAMP,NC2VEC,1,LUPRI) 380 ENDIF 381 IF (.TRUE.) THEN 382 XNJ = X11 + X12 + X21 + X22 383 WRITE(LUPRI,*)' ' 384 WRITE(LUPRI,*)' NORM OF FIN. DIFF. D*TA*TB*TC-Matrix.', SQRT(XNJ) 385 WRITE(LUPRI,*)' ' 386 WRITE(LUPRI,*)' NORM OF 11 PART OF FD. D*TA*TB*TC-mat.: ', 387 & SQRT(X11) 388 WRITE(LUPRI,*)' NORM OF 21 PART OF FD. D*TA*TB*TC-mat.: ', 389 & SQRT(X21) 390 WRITE(LUPRI,*)' NORM OF 12 PART OF FD. D*TA*TB*TC-mat.: ', 391 & SQRT(X12) 392 WRITE(LUPRI,*)' NORM OF 22 PART OF FD. D*TA*TB*TC-mat.: ', 393 & SQRT(X22) 394 ENDIF 395C 396C-------------------------------------- 397C Calculate Matrix times Ty vector. 398C-------------------------------------- 399C 400 CALL DGEMV('N',NTAMP,NTAMP,ONE,WORK(KF),NTAMP,TYAM,1, 401 * ZERO,RESULT,1) 402 403C-------------------------------------- 404C scale diagonal with 1/2: 405C-------------------------------------- 406C CALL CCLR_DIASCL(RESULT(NT1AM(1)+1),TWO,1) 407 408 WRITE (LUPRI,*) 'NTAMP:',NTAMP 409 WRITE (LUPRI,*) 'NORM^2 OF TXAM VECTOR:', 410 * DDOT(NT1AM(1)+NT2AM(1),TXAM,1,TXAM,1) 411 WRITE (LUPRI,*) 'NORM^2 OF TYAM VECTOR:', 412 * DDOT(NT1AM(1)+NT2AM(1),TYAM,1,TYAM,1) 413 WRITE (LUPRI,*) 'NORM^2 OF RESULT VECTOR:', 414 * DDOT(NTAMP,RESULT,1,RESULT,1) 415 416C 417C------------------------------------------------- 418C Restore the CC reference amplitudes on disk. 419C------------------------------------------------- 420C 421 LUTAM = -1 422 CALL GPOPEN(LUTAM,'TAM_SAV','UNKNOWN',' ','UNFORMATTED',IDUMMY, 423 * .FALSE.) 424 REWIND(LUTAM) 425 READ(LUTAM) (WORK(KC1AM + I -1 ) , I = 1, NT1AMX) 426 READ(LUTAM) (WORK(KC2AM + I -1 ) , I = 1, NT2AMX) 427 CALL GPCLOSE(LUTAM,'DELETE') 428C 429 IOPT = 3 430 CALL CC_WRRSP('R0',0,1,IOPT,MODEL,DUMMY,WORK(KC1AM), 431 * WORK(KC2AM),WORK(KEND2),LWRK2) 432C 433 RSPIM = .TRUE. 434 CALL CCRHSN(WORK(KRHO1D),WORK(KRHO2D),WORK(KC1AM), 435 & WORK(KC2AM), 436 * WORK(KEND2),LWRK2,'XXX') 437C 438 IF (IPRINT .GT. 10) THEN 439 CALL AROUND(' END OF CC_FDD ') 440 ENDIF 441C 442 RETURN 443 END 444*=====================================================================* 445