1!-------------------------------------------------------------------- 2 SUBROUTINE CC_LANCZOS_LR(LABEL,J,ER,EI,R,L,F,U,V) 3!-------------------------------------------------------------------- 4! Given a set of Lanzcos T eigenvalues and eigenvector this 5! calculates the linear response function for 6! given frequencies and gammas 7! J: Chain Length 8! ER: Eigenvalues, real part 9! EI: Eigenvalues, imag part 10! R: Eigenvectors right (stored as cols) 11! L: Eigenvectors Left (stored as cols, e.i. eigvectors of T^T) 12! F: F-matrix "trans" (ikke double q-transformed) 13! U, V are normalization factors 14! 15! FIXME: todo list: 16! 1) rethink handling of complex roots 17! 2) make loops more efficient (vectors maybe?) 18! 3) "weight" thing 19! 20! Sonia Coriani, 2010-2012 21!-------------------------------------------------------------------- 22 IMPLICIT NONE 23#include "priunit.h" 24#include "cclrlancz.h" 25#include "codata.h" 26 27 INTEGER J,IFREQ,NFREQ 28#if defined (SYS_CRAY) 29 REAL U,V 30 REAL FSTART,FSTOP,FSTEP 31 REAL ER(J),EI(J),R(J,J),L(J,J),F(J,J) 32 REAL ZERO, ONE, TWO, PDIFFI,PSUMI,PSUMK,XGAMMA 33 REAL LR_OM_RE,LR_OM_RE_T,LR_OM_RE_F 34 REAL LR_OM_IM,LR_OM_IM_T,LR_OM_IM_F 35 REAL DPI, DPK,FACTOR,FACTOR2, FREQ, DMI 36#else 37 DOUBLE PRECISION U,V 38 DOUBLE PRECISION FSTART,FSTOP,FSTEP 39 DOUBLE PRECISION ER(J),EI(J),R(J,J),L(J,J),F(J,J) 40 DOUBLE PRECISION ZERO, ONE, TWO, PDIFFI,PSUMI,PSUMK,XGAMMA 41 DOUBLE PRECISION LR_OM_RE,LR_OM_RE_T,LR_OM_RE_F 42 DOUBLE PRECISION LR_OM_IM,LR_OM_IM_T,LR_OM_IM_F 43 DOUBLE PRECISION DPI, DPK,FACTOR,FACTOR2, FREQ, DMI 44 DOUBLE PRECISION XGAM_PLUS,XGAM_MINUS,FACTORR,FACTORI 45 Double precision add_to_real, add_to_im 46 Double precision tot_add_real, tot_add_im 47 double precision weight(j,2) 48#endif 49 LOGICAL NOIMAG 50 INTEGER I, K, IDAMP, IDUMMY 51! 52 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 53 !Specifics of the file containing the Lanczos abs spectrum 54 CHARACTER*7 FNSPECTRUM 55 PARAMETER (FNSPECTRUM='CCSPECX') 56 INTEGER LUSPEC 57 ! 58 CHARACTER*(8) LABEL 59 LOGICAL LOCDBG 60 PARAMETER (LOCDBG=.false.) 61! 62 LOGICAL SKIPNEXT 63 SKIPNEXT=.FALSE. 64! 65! Set these using input data 66! 67 FSTART = FREQ_RANGE(1) 68 FSTOP = FREQ_RANGE(2) 69 FSTEP = FREQ_RANGE(3) 70 XGAMMA = DAMPING(1) 71 LR_OM_RE_T=zero 72 LR_OM_IM_T=zero 73 !CALL DZERO(weight,j*2) 74 75 NFREQ=1+INT((FSTOP-FSTART)/FSTEP) 76 WRITE (LUPRI,*) "1. FREQ", FSTART," NFREQ:",NFREQ, " STEP:", FSTEP 77! 78! Open unit with results 79! 80 LUSPEC=-1 81 CALL GPOPEN(LUSPEC,FNSPECTRUM,'NEW',' ','FORMATTED', 82 & IDUMMY,.FALSE.) 83 84 !initialize to zero 85 add_to_real = zero 86 add_to_im = zero 87 tot_add_real = zero 88 tot_add_im = zero 89C 90C Loop over chain 91C 92 DO IDAMP=1,NDAMP 93 XGAMMA = DAMPING(IDAMP) 94 95 WRITE (LUPRI,*) " FOR OPERATOR LABEL = ", LABEL 96 WRITE (LUPRI,'(1X,A27,F12.8,A3,I10)') 97 * " Response fun with gamma = ",XGAMMA, ' J=',J 98 WRITE(LUPRI,*) " Frequency (eV) & Real (au)& Imaginary(au)", 99 & "& sigma^I (au)" 100 WRITE(LUPRI,'(X,"---------------------------")') 101! 102!Also on file LUSPEC 103! 104 WRITE (LUSPEC,*) " FOR OPERATOR LABEL = ", LABEL 105 WRITE (LUSPEC,'(1X,A27,F12.8,A3,I10)') 106 * " Response fun with gamma = ",XGAMMA, ' J=',J 107 WRITE(LUSPEC,*) " Frequency (eV) & Real (au)& Imaginary(au)", 108 & "& sigma^I (au)" 109 WRITE(LUSPEC,'(X,"---------------------------")') 110 111 DO IFREQ=1,NFREQ 112 113 FREQ = FSTART+FSTEP*(IFREQ-1) 114 LR_OM_RE=ZERO 115 LR_OM_IM=ZERO 116 LR_OM_RE_F=ZERO 117 LR_OM_IM_F=ZERO 118! 119!reinitialize inside freq loop 120! 121 add_to_real = zero 122 add_to_im = zero 123 tot_add_real = zero 124 tot_add_im = zero 125 126 DO I=1,J 127 IF (SKIPNEXT) THEN 128 SKIPNEXT = .FALSE. 129 ELSE 130 !check on imaginary eigenvalue = 0 131 !IF (ABS(EI(I)).EQ.ZERO) THEN 132 IF (EI(I).EQ.ZERO) THEN 133 FACTOR=U*V*R(1,I)*L(1,I) 134C 135C No imaginary eigenvector, calc eta*t conts 136C 137 PDIFFI = FREQ-ER(I) 138 PSUMI = FREQ+ER(I) 139 !denominators 140 DMI = PDIFFI*PDIFFI + XGAMMA*XGAMMA 141 DPI = PSUMI*PSUMI + XGAMMA*XGAMMA 142 LR_OM_RE = LR_OM_RE + FACTOR*(PDIFFI/DMI-PSUMI/DPI) 143 LR_OM_IM = LR_OM_IM - FACTOR*XGAMMA*(ONE/DMI-ONE/DPI) 144 !weight(i,1) = ER(I) 145 !weight(i,2) = R(1,I)*L(1,I) 146C 147C Include F-part cont 148C 149 FACTOR = V*V*L(1,I) 150 DO K=1,J 151 FACTOR2 = FACTOR*L(1,K)*F(I,K) 152 if (abs(ei(k)).eq.zero) then 153 PSUMK = FREQ+ER(K) 154 !denominator 155 DPK = PSUMK*PSUMK + XGAMMA*XGAMMA 156 LR_OM_RE_F = LR_OM_RE_F + FACTOR2* 157 & (PSUMK*PDIFFI-XGAMMA*XGAMMA)/(DMI*DPK) 158 LR_OM_IM_F = LR_OM_IM_F - FACTOR2* 159 & XGAMMA*(PSUMK+PDIFFI)/(DMI*DPK) 160 else 161 !write(lupri,*)'LANCZOS, do not know, index=',k 162 !SHOULD DECIDE HOW TO HANDLE THIS 163! FIXME 1) 164 end if 165 ENDDO 166 167 ELSE 168 169 FACTORR=U*V*(R(1,I)*L(1,I)-R(1,I+1)*L(1,I+1)) 170 FACTORI=U*V*(R(1,I+1)*L(1,I)+R(1,I)*L(1,I+1)) 171! 172! Include also states with imaginary excitation energy 173! 174 PDIFFI = FREQ-ER(I) 175 PSUMI = FREQ+ER(I) 176 xgam_minus = xgamma - EI(I) 177 xgam_plus = xgamma + EI(I) 178 !denominators 179 DMI = PDIFFI*PDIFFI + XGAM_minus*XGAM_minus 180 DPI = PSUMI*PSUMI+ XGAM_plus*XGAM_plus 181 ! 182 !real response function (eta part) 183 ! 184 add_to_real = FACTORR*(PDIFFI/DMI-PSUMI/DPI) 185 & + FACTORI*(xgam_minus/DMI-xgam_plus/DPI) 186 tot_add_real = tot_add_real + add_to_real 187 LR_OM_RE = LR_OM_RE + FACTORR*(PDIFFI/DMI-PSUMI/DPI) 188 & + FACTORI*(xgam_minus/DMI-xgam_plus/DPI) 189 ! 190 !imaginary response function (eta part) 191 ! 192 add_to_im = FACTORR*(XGAM_minus/DMI- 193 & xgam_plus/DPI) + FACTORI*(PDIFFI/DMI 194 & - PSUMI/DPI) 195 tot_add_im = tot_add_im + add_to_im 196 LR_OM_IM = LR_OM_IM - FACTORR*(XGAM_minus/DMI- 197 & xgam_plus/DPI) + FACTORI*(PDIFFI/DMI 198 & - PSUMI/DPI) 199C 200C Include F-part cont 201C 202! FIXME 1) 203! write(lupri,*)'WARNING: CMPLX ROOT, SKIP F PART' 204! FACTOR = V*V*L(1,I) 205! DO K=1,J 206! FACTOR2 = FACTOR*L(1,K)*F(I,K) 207! PSUMK = FREQ+ER(K) 208! !denominator 209! DPK = PSUMK*PSUMK + XGAMMA*XGAMMA 210! LR_OM_RE_F = LR_OM_RE_F + FACTOR2* 211! & (PSUMK*PDIFFI-XGAMMA*XGAMMA)/(DMI*DPK) 212! LR_OM_IM_F = LR_OM_IM_F - FACTOR2* 213! & XGAMMA*(PSUMK+PDIFFI)/(DMI*DPK) 214! ENDDO 215 SKIPNEXT=.true. 216 ENDIF 217 end if !skipnext 218C 219 ENDDO 220 221 if (locdbg) then 222 WRITE (LUPRI,*)'tot_add_real', tot_add_real 223 WRITE (LUPRI,*)'tot_add_im', tot_add_im 224 WRITE (LUPRI,*)'Freq; F contrib Re', FREQ, LR_OM_RE_F 225 WRITE (LUPRI,*)'Freq; F contrib Im', FREQ, LR_OM_IM_F 226 end if 227 228 LR_OM_RE_T = LR_OM_RE - LR_OM_RE_F 229 LR_OM_IM_T = LR_OM_IM - LR_OM_IM_F 230 231 WRITE (LUSPEC,'(6F16.8)') 232 & FREQ*XTEV, LR_OM_RE_T, LR_OM_IM_T, -LR_OM_IM_T*FREQ, 233 & LR_OM_RE,LR_OM_IM 234 235 IF (IFREQ.LE.10) then 236 !dump first 10 values on output file as well 237 WRITE (LUPRI,'(6F16.8)') 238 & FREQ*XTEV, LR_OM_RE_T, LR_OM_IM_T, -LR_OM_IM_T*FREQ, 239 & LR_OM_RE,LR_OM_IM 240 end if 241 242 ENDDO 243 WRITE(LUPRI,'(X,"--- remaining values on LUSPEC file ---")') 244 WRITE(LUSPEC,'(X,"---------------------------")') 245 ENDDO 246 247 CALL GPCLOSE(LUSPEC,'KEEP') 248 249 RETURN 250 END 251!------------------------------------------------------------ 252 SUBROUTINE CC_BIORTOCOMP(J,ER,EI,EVR,EVL,WORK,LWORK) 253!------------------------------------------------------------ 254! From a DGEEV output, take eigenvectors with unit norm and 255! make them biorthogonal so that <R|R>=1,<L|R>=1, by 256! scaling the left eigenvectors with 1/<L|R>. 257! NB: Take special care with complex that are assumed 258! to come with eigenvectors order as 259! Right eigenvector for w_k_r + iw_k_i is R(k) + i R(k+1) 260! Right eigenvector for w_k_r - iw_k-I is R(k) - i R(k+1) 261! and similar for left eigenvectors 262! Ove Christiansen & Sonia Coriani, November 2010 263! 264! J = chain length, equal to length of each eigenvector 265! and to nr of eigenvalues 266! ER = real parts of eigenvalues 267! EI = imaginary parts of eigenvalues 268! EVR = Right eigenvectors (DGEEV solutions) 269! EVL = left eigenvectors (DGEEV solutions) 270C-------------------------------------------------------------------- 271 IMPLICIT NONE 272#include "priunit.h" 273 INTEGER J,LWORK 274#if defined (SYS_CRAY) 275 REAL ER(J),EI(J),EVR(J,J),EVL(J,J),WORK(LWORK) 276 REAL XNORM,THRESH 277 REAL ONE,ZERO,DDOT 278#else 279 DOUBLE PRECISION ER(J),EI(J),EVR(J,J),EVL(J,J),WORK(LWORK) 280 DOUBLE PRECISION XNORM,THRESH 281 DOUBLE PRECISION ONE,ZERO,DDOT 282#endif 283 PARAMETER (ONE=1.0d0,ZERO=0.0d0,THRESH=1.0d-16) 284 INTEGER K, IDAMP, KLRO, KLIO 285 LOGICAL SKIPNEXT 286 287 SKIPNEXT=.FALSE. 288 !write(lupri,*)'I AM INSIDE BIORTO COMP' 289 290 KLRO = 1 !L real Output 291 KLIO = KLRO + J !L imaginary Output 292 IF (KLIO+J .GT. LWORK) 293 & CALL QUIT("Too little space in CC_BIORTOCOMP") 294 295 DO K=1,J 296 IF (SKIPNEXT) THEN 297 SKIPNEXT = .FALSE. 298 ELSE 299 IF (EI(K).EQ.ZERO) THEN 300!=================================== 301! Real eigenvalue 302!=================================== 303 XNORM = DDOT(J,EVL(1,K),1,EVR(1,K),1) 304 !write(lupri,*)'Biortho of RL', xnorm 305 IF (ABS(XNORM).LE. THRESH) THEN 306 CALL QUIT("Error in CC_BIORTOCOMP") 307 ENDIF 308 XNORM = ONE/XNORM 309 !write(lupri,*)'Biortho of RL', xnorm 310 CALL DSCAL(J,XNORM,EVL(1,K),1) 311 !XNORM = DDOT(J,EVL(1,K),1,EVR(1,K),1) 312 !write(lupri,*)'Biortho of RL dopo',xnorm,' K=', K, 'J=',J 313 ELSE 314!============================= 315! complex pair 316!============================= 317 write(lupri,*)'WARNING: COMPLEX PAIR alternative norm' 318 call flshfo(lupri) 319 CALL CC_BIORTOCOMP2(J,EVR(1,K),EVR(1,K+1), 320 * EVL(1,K),EVL(1,K+1), 321 * WORK(KLRO),WORK(KLIO)) 322 !SUBROUTINE CC_BIORTOCOMP2(N,EVRR,EVRI,EVLR,EVLI,EVLRO,EVLIO) 323 call flshfo(lupri) 324 CALL DCOPY(J,WORK(KLRO),1,EVL(1,K),1) 325 CALL DCOPY(J,WORK(KLIO),1,EVL(1,K+1),1) 326 327 SKIPNEXT =.TRUE. 328 ENDIF 329 ENDIF 330 ENDDO 331 332 RETURN 333 END 334!-------- 335 SUBROUTINE CC_BIORTOCOMP1(N,EVRR,EVRI,EVLR,EVLI,EVLRO,EVLIO) 336C-------------------------------------------------------------------- 337C OC,03/11/2010 338C Take eigenvectors in with unit norm. 339C Create new left vectors eigenvectors such that they are biorthog 340C to right eigenvectors. 341C 342C EVRR: Right eigenvector, Real part. 343C EVRI: Right eigenvector, Imaginary part. 344C EVLR: Eigenvector, left, real part. 345C EVLI: Eigenvector, left, imag. part. 346C EVLRO: Eigenvectors, left, real part. OUTPUT 347C EVLIO: Eigenvectors, left, imag. part. OUTPUT 348C N : length of vector 349C (we pass one eigenvector at a time) 350C-------------------------------------------------------------------- 351C 352 IMPLICIT NONE 353#include "priunit.h" 354 INTEGER N 355#if defined (SYS_CRAY) 356 REAL EVRR(N),EVRI(N),EVLR(N),EVLI(N),EVLRO(N),EVLIO(N) 357 REAL CR,CI,CNSQ,DDOT 358#else 359 !input 360 DOUBLE PRECISION EVRR(N),EVRI(N),EVLR(N),EVLI(N) 361 !output 362 Double precision EVLRO(N),EVLIO(N) 363 !local 364 DOUBLE PRECISION CR,CI,CNSQ,DDOT 365#endif 366 367 write(lupri,*)'I AM INSIDE BIORTCOMP1' 368 call flshfo(lupri) 369 write(lupri,*)'<rr|rr>=',ddot(N,EVRR,1,EVRR,1) 370 write(lupri,*)'<ri|ri>=',ddot(N,EVRI,1,EVRI,1) 371 write(lupri,*)'<lr|lr>=',ddot(N,EVLR,1,EVLR,1) 372 write(lupri,*)'<li|li>=',ddot(N,EVLI,1,EVLI,1) 373C 374 write(lupri,*)'<lr|Rr>=',ddot(N,EVLR,1,EVRR,1) 375 write(lupri,*)'<li|Ri>=',ddot(N,EVLI,1,EVRI,1) 376 write(lupri,*)'<lr|Ri>=',ddot(N,EVLR,1,EVRI,1) 377 write(lupri,*)'<li|Rr>=',ddot(N,EVLI,1,EVRR,1) 378 !<ReR|ReL>-<ImR|ImL> 379 CR = DDOT(N,EVLR,1,EVRR,1)-DDOT(N,EVLI,1,EVRI,1) 380 write(lupri,*)'CR=', CR 381 !<ReR|ImL>+<ImR|ReL> 382 CI = DDOT(N,EVLI,1,EVRR,1)+DDOT(N,EVLR,1,EVRI,1) 383 write(lupri,*)'CI=', CI 384 CNSQ=CR*CR+CI*CI 385 write(lupri,*)'CNSQ', CNSQ 386 write(lupri,*)'CI/CNSQ', CI/CNSQ, 'CR/CNSQ', Cr/Cnsq 387!--- 388! !<ReL|ReI>+<ImL|ImR> 389! CRp = DDOT(N,EVLR,1,EVRR,1)+DDOT(N,EVLI,1,EVRI,1) 390! !<ReL|ImR>-<ImL|ReR> 391! CIm = DDOT(N,EVLR,1,EVRI,1)-DDOT(N,EVLI,1,EVRR,1) 392! CNmSQ=CRp*CRp+CIm*CIm 393C 394 CALL DZERO(EVLRO,N) 395 CALL DZERO(EVLIO,N) 396C 397 !write(lupri,*)'CR/CNSQ=',CNSQ 398 !call flshfo(lupri) 399 CALL DCOPY(N,EVLR,1,EVLRO,1) 400 CALL DCOPY(N,EVLI,1,EVLIO,1) 401 CALL DSCAL(N,(CR/CNSQ),EVLRO,1) 402 CALL DSCAL(N,(CR/CNSQ),EVLIO,1) 403C 404 CALL DAXPY(N,(CI/CNSQ),EVLI,1,EVLRO,1) 405 CALL DAXPY(N,(-CI/CNSQ),EVLR,1,EVLIO,1) 406! 407 WRITE(LUPRI,*)"This should be one ", 408 * DDOT(N,EVLRO,1,EVRR,1)-DDOT(N,EVLIO,1,EVRI,1) 409 call flshfo(lupri) 410 WRITE(LUPRI,*)"This should be zero ", 411 * DDOT(N,EVLRO,1,EVRI,1)+DDOT(N,EVLIO,1,EVRR,1) 412 call flshfo(lupri) 413 414 WRITE(LUPRI,*)"<new LrO|Rr>", DDOT(N,EVLRO,1,EVRR,1) 415 WRITE(LUPRI,*)"<new LiO|Ri>", DDOT(N,EVLiO,1,EVRi,1) 416 WRITE(LUPRI,*)"<new LrO|Ri>", DDOT(N,EVLRO,1,EVRi,1) 417 WRITE(LUPRI,*)"<new LiO|Rr>", DDOT(N,EVLiO,1,EVRR,1) 418 call flshfo(lupri) 419 420 RETURN 421 END 422!-------- 423 SUBROUTINE CC_BIORTOCOMP2(N,EVRR,EVRI,EVLR,EVLI,EVLRO,EVLIO) 424C-------------------------------------------------------------------- 425C SC,09/12/2010 426C Biorthonormalization of complex left eigenvectors 427C Create new left vectors eigenvectors such that they are biorthog 428C to right eigenvectors. Based on BIORT1. 429C Complex input left eigenv to dot on R is taken as <Lr|-i<Li| 430C Complex left eigenv in output is taken as <Lr_o|+i<Li_o| 431C 432C EVRR: Right eigenvector, Real part. 433C EVRI: Right eigenvector, Imaginary part. 434C EVLR: Eigenvector, left, real part. 435C EVLI: Eigenvector, left, imag. part. 436C EVLRO: Eigenvectors, left, real part. OUTPUT 437C EVLIO: Eigenvectors, left, imag. part. OUTPUT 438C N : length of vector 439C (we pass one eigenvector at a time) 440C-------------------------------------------------------------------- 441C 442 IMPLICIT NONE 443#include "priunit.h" 444 INTEGER N 445#if defined (SYS_CRAY) 446 REAL EVRR(N),EVRI(N),EVLR(N),EVLI(N),EVLRO(N),EVLIO(N) 447 REAL CR,CI,CNSQ,DDOT 448#else 449 !input 450 DOUBLE PRECISION EVRR(N),EVRI(N),EVLR(N),EVLI(N) 451 !output 452 Double precision EVLRO(N),EVLIO(N) 453 !local 454 DOUBLE PRECISION CR,CI,CNSQ,DDOT 455#endif 456 457 write(lupri,*)'INSIDE BIORTCOMP1' 458 call flshfo(lupri) 459 write(lupri,*)'<rr|rr>=',ddot(N,EVRR,1,EVRR,1) 460 write(lupri,*)'<ri|ri>=',ddot(N,EVRI,1,EVRI,1) 461 write(lupri,*)'<lr|lr>=',ddot(N,EVLR,1,EVLR,1) 462 write(lupri,*)'<li|li>=',ddot(N,EVLI,1,EVLI,1) 463C 464 write(lupri,*)'<lr|Rr>=',ddot(N,EVLR,1,EVRR,1) 465 write(lupri,*)'<li|Ri>=',ddot(N,EVLI,1,EVRI,1) 466 write(lupri,*)'<lr|Ri>=',ddot(N,EVLR,1,EVRI,1) 467 write(lupri,*)'<li|Rr>=',ddot(N,EVLI,1,EVRR,1) 468 !<ReL|ReR>+<ImL|ImR> 469 CR = DDOT(N,EVLR,1,EVRR,1)+DDOT(N,EVLI,1,EVRI,1) 470 write(lupri,*)'CR=', CR 471 !<ReL|ImR>-<ImL|ReR> 472 CI = DDOT(N,EVLR,1,EVRI,1)-DDOT(N,EVLI,1,EVRR,1) 473 write(lupri,*)'CI=', CI 474 CNSQ=CR*CR+CI*CI 475 write(lupri,*)'CNSQ', CNSQ 476 write(lupri,*)'CI/CNSQ', CI/CNSQ, 'CR/CNSQ', Cr/Cnsq 477C 478 CALL DZERO(EVLRO,N) 479 CALL DZERO(EVLIO,N) 480C 481 !call flshfo(lupri) 482 CALL DCOPY(N,EVLR,1,EVLRO,1) 483 CALL DCOPY(N,EVLI,1,EVLIO,1) 484 CALL DSCAL(N,(CR/CNSQ),EVLRO,1) 485 CALL DSCAL(N,(-CR/CNSQ),EVLIO,1) 486C 487 CALL DAXPY(N,(-CI/CNSQ),EVLI,1,EVLRO,1) 488 CALL DAXPY(N,(-CI/CNSQ),EVLR,1,EVLIO,1) 489! 490 WRITE(LUPRI,*)"This should be one ", 491 * DDOT(N,EVLRO,1,EVRR,1)-DDOT(N,EVLIO,1,EVRI,1) 492 call flshfo(lupri) 493 WRITE(LUPRI,*)"This should be zero ", 494 * DDOT(N,EVLRO,1,EVRI,1)+DDOT(N,EVLIO,1,EVRR,1) 495 call flshfo(lupri) 496 497 WRITE(LUPRI,*)"<new LrO|Rr>", DDOT(N,EVLRO,1,EVRR,1) 498 WRITE(LUPRI,*)"<new LiO|Ri>", DDOT(N,EVLiO,1,EVRi,1) 499 WRITE(LUPRI,*)"<new LrO|Ri>", DDOT(N,EVLRO,1,EVRi,1) 500 WRITE(LUPRI,*)"<new LiO|Rr>", DDOT(N,EVLiO,1,EVRR,1) 501 call flshfo(lupri) 502 503 RETURN 504 END 505!-------- 506 SUBROUTINE CC_LANCZOS_Ftrans(ioptrw,Lchain,Lspace,IFQTRAN,XVEC, 507 & EVR,Ftrans,WORK,LWORK,ER,EI,ISYHOP) 508C------------------------------------------------------------------- 509C SC,03/11/2010 510C Read in the FQ vectors and generate the trans^F matrix 511! Lchain is length of Chain 512! Lspace is length of exc. space 513! Xvec is QR 514! EVR is eigenvector 515! ER is (real part of) eigenvalue 516! EI is (imaginary part of) eigenvalue 517C-------------------------------------------------------------------- 518C 519 IMPLICIT NONE 520#include "priunit.h" 521#include "ccsdsym.h" 522 INTEGER Lchain,Lspace,IFQTRAN(3,*),LWORK,ioptrw 523#if defined (SYS_CRAY) 524 REAL XVEC(Lspace,Lchain),EVR(Lchain,Lchain),ER(Lchain), EI(Lchain) 525 REAL Ftrans(Lchain,Lchain), FQX(Lchain,Lchain) 526 REAL DDOT, WORK(LWORK) 527#else 528 DOUBLE PRECISION XVEC(Lspace,Lchain),EVR(Lchain,Lchain),ER(Lchain) 529 DOUBLE PRECISION EI(Lchain) 530 DOUBLE PRECISION Ftrans(Lchain,Lchain), FQX(Lchain,Lchain) 531 DOUBLE PRECISION DDOT, WORK(LWORK) 532 DOUBLE PRECISION ONE,ZERO 533#endif 534 INTEGER kFq1,kFq2,kend,lend,isyhop,ivec, n2vec 535 CHARACTER*(10) MODEL 536 PARAMETER (ONE=1.0d0,ZERO=0.0d0) 537 538 539! CALL AROUND('trans^F print input stuff') 540! DO J = 1, Lchain 541! WRITE (LUPRI,'(/A,I4,A,1P,D15.6/A)') 542! * ' - R EIGENVALUE NO.',J, ' - EIGENVALUE',ER(J), 543! * ' - R EIGENVECTOR :' 544! WRITE (LUPRI,'(5F15.8)') (EVR(I,J),I=1,Lchain) 545! WRITE (LUPRI,'(/A,I4,A,1P,D15.6/A)') 546! * ' - EIGENVALUE NO.',J, ' - EIGENVALUE',ER(J), 547! * ' - QR EIGENVECTOR :' 548! WRITE (LUPRI,'(5F15.8)') (XVEC(I,J),I=1,Lspace) 549! END DO 550 n2vec = 1 551 if (ioptrw.eq.1) n2vec=0 552 553 CALL DZERO(FQX,Lchain*LCHAIN) 554 kFq1 = 1 555 !kfq2 = kFq1 + NT1AMX 556 !kend = kFq2 + NT2AMX 557 kfq2 = kFq1 + NT1AM(isyhop) 558 kend = kFq2 + NT2AM(isyhop) 559 lend = lwork - kend 560 if (lend.lt.0) call quit('Insufficient memory for trans^F') 561 !isyhop=1 562 do i = 1, Lchain 563 ivec = ifqtran(3,i) 564 CALL CC_RDRSP('FQ ',IVEC,ISYHOP,IOPTRW,MODEL, 565 & WORK(KFQ1),WORK(KFQ2)) 566 567! CALL AROUND('FQ transformed vector ') 568! CALL CC_PRP(WORK(KFQ1),WORK(KFQ2),ISYHOP,1,N2VEC) 569 570 do j = 1, Lchain 571 if (EI(j).ne.zero) then 572 !write(lupri,*) 'WARNING: FQX of cmplx, EI(j)=', j, EI(j) 573 end if 574 FQX(j,i) = ddot(Lspace,WORK(KFQ1),1,XVEC(1,j),1) 575 end do 576 end do 577! CALL AROUND('--- The FQQR matrix ---') 578! CALL OUTPUT(FQX,1,Lchain,1,Lchain,Lchain,Lchain,1,LUPRI) 579 580 581 CALL DGEMM('N','N',Lchain,Lchain,Lchain, 582 * ONE,FQX,Lchain, 583 * EVR,Lchain,ZERO, 584 * Ftrans,Lchain) 585! CALL AROUND('--- The trans^F matrix ---') 586! CALL OUTPUT(Ftrans,1,Lchain,1,Lchain,Lchain,Lchain,1,LUPRI) 587 588 RETURN 589 END 590 591C /* Deck cc_pram_lanczos */ 592 SUBROUTINE CC_PRAM_lanczos(CAM,PT1,ISYMTR,LGRS,work,lwork) 593C 594C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 595C Based on Ove's cc_pram 596C 597C Purpose: Writes out vector: 598C %T1 and %T2 and ||T1||/||T2|| 599C makes an analysis of contributions from given 600C occupied orbital 601C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 602C 603#include "implicit.h" 604C 605 dimension work(lwork) 606 PARAMETER (TWO = 2.0D0, THPRT = 1.0D-9, HUNDRED = 100.0D0) 607 PARAMETER (PT1THR = 75.0D0) 608C 609#include "ccorb.h" 610#include "ccsdsym.h" 611#include "ccsdinp.h" 612#include "priunit.h" 613Cholesky 614#include "maxorb.h" 615#include "ccdeco.h" 616C 617 LOGICAL lxray 618 LOGICAL CCSEFF 619Cholesky 620C 621C 622 LOGICAL LGRS 623 DIMENSION CAM(*) 624C 625 INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 626C 627Cholesky 628 CCSEFF = CCS .OR. (CHOINT.AND.CC2) 629Cholesky 630C 631C------------------------ 632C Add up the vectors. 633C------------------------ 634C 635 C1NOSQ = 0.0D0 636 C2NOSQ = 0.0D0 637 KC1 = 1 638 KC2 = KC1 + NT1AM(ISYMTR) 639 !<t1|t1> 640 C1NOSQ = DDOT(NT1AM(ISYMTR),CAM(KC1),1,CAM(KC1),1) 641Chol IF (.NOT. CCS) C2NOSQ = DDOT(NT2AM(ISYMTR),CAM(KC2),1,CAM(KC2),1) 642 IF (.NOT. CCSEFF) 643 !<t2|t2> 644 & C2NOSQ = DDOT(NT2AM(ISYMTR),CAM(KC2),1,CAM(KC2),1) 645 CNOSQ = C1NOSQ + C2NOSQ 646C 647 IF (.NOT. (CCSEFF.OR.CCP2)) THEN 648C 649 WRITE(LUPRI,'(//10X,A)') 650 * 'CC_PRAM:Overall Contribution of the Different Components' 651 WRITE(LUPRI,'(10X,A//)') 652 * '--------------------------------------------------------' 653 WRITE(LUPRI,'(/10X,A,10X,F10.4,A)') 654 * 'Single Excitation Contribution <t1|t1>/<t|t>*100: ', 655 * (C1NOSQ/CNOSQ)*HUNDRED,' %' 656 WRITE(LUPRI,'(/10X,A,10X,F10.4,A)') 657 * 'Double Excitation Contribution <t2|t2>/<t|t>*100: ', 658 * (C2NOSQ/CNOSQ)*HUNDRED,' %' 659 WRITE(LUPRI,'(/10X,A,10X,F10.4)') 660 * '||T1||/||T2||=sqrt(<t1|t1>/<t2|t2>) : ', 661 * SQRT(C1NOSQ/C2NOSQ) 662 IF (LGRS) WRITE(LUPRI,'(/10X,A,10X,F10.4)') 663 * 'Tau1 diagnostic : ', 664 * SQRT(C1NOSQ/(TWO*DFLOAT(NRHFT))) 665 PT1 = (C1NOSQ/CNOSQ)*HUNDRED 666 ELSE 667 PT1 = HUNDRED 668 ENDIF 669 WRITE(LUPRI,'(/10X,A,10X,F10.4)') 670 * 'Norm of Total Amplitude Vector : ',SQRT(CNOSQ) 671C 672 CALL FLSHFO(LUPRI) 673C 674C---------------------------------------------- 675C Initialize threshold etc from Printlevel. 676C---------------------------------------------- 677C 678 NL = MAX(1,2*IPRINT) 679C 680 CNOSQ = MAX(CNOSQ,THPRT) 681C 682 THR1 = SQRT(C1NOSQ/CNOSQ)/NL 683 THR2 = SQRT(C2NOSQ/CNOSQ)/NL 684 THR1 = MAX(THR1,1.0D-02) 685 THR2 = MAX(THR2,1.0D-02) 686 SUMOFP = 0.0D00 687C 688 IF (DEBUG) THR1 = 0.0D0 689C 690C------------------------------------- 691C Loop until a few is Printed out. 692C------------------------------------- 693C 694C 695C--------------------------------------- 696C Loop through One excitation part. 697C--------------------------------------- 698C 699 WRITE(LUPRI,'(//A)') 700 * ' +==============================================' 701 * //'===============================+' 702 WRITE(LUPRI,'(1X,A)') 703 * '| symmetry| orbital index | Excitation Numbers' 704 * //' | Amplitude |' 705 WRITE(LUPRI,'(1X,A)') 706 * '| Index | a b i j | NAI NBJ |' 707 * //' NAIBJ | |' 708 WRITE(LUPRI,'(A)') 709 * ' +==============================================' 710 * //'===============================+' 711C 712 ISYMAI = MULD2H(ISYMTR,ISYMOP) 713C 714 1 CONTINUE 715 N1 = 0 716C 717 DO 100 ISYMA = 1,NSYM 718C 719 ISYMI = MULD2H(ISYMAI,ISYMA) 720C 721 DO 110 I = 1,NRHF(ISYMI) 722C 723 MI = IORB(ISYMI) + I 724C 725 DO 120 A=1,NVIR(ISYMA) 726C 727 NAI = IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I-1) + A 728C 729 MA = IORB(ISYMA) + NRHF(ISYMA) + A 730C 731 IF (ABS(CAM(NAI)) .GE. THR1 ) THEN 732C 733 WRITE(LUPRI,9990) ISYMA,ISYMI,A,I,NAI,CAM(NAI) 734C 735 N1 = N1 + 1 736 SUMOFP = SUMOFP + CAM(NAI)*CAM(NAI) 737C 738 ENDIF 739C 740 120 CONTINUE 741 110 CONTINUE 742 100 CONTINUE 743C 744 IF ((N1.LT.1).AND.(SQRT(C1NOSQ/CNOSQ).GT.1.0D-3)) THEN 745 THR1 = THR1/5.0D0 746 GOTO 1 747 ENDIF 748 749 !skod test with no symmetry 750! maxlength=0 751! do isym=1,nsym 752! maxnocc=max(maxnocc,nrhf(isym)) 753! end do 754 755 IF (PT1.ge.pt1thr) then 756 757! maxocc = 0 758! do isym=1,nsym 759! maxocc = max(maxocc,nrhf(isym)) 760! end do 761 762 kstart = 1 763 kend = nrhft 764 lend = lwork - kend 765C 766 call dzero(work(kstart),nrhft) 767! DO I = 1,NRHFt 768! MI = IORB(1) + I 769! NA = NVIRT*(I-1) + 1 770! work(mi) = ddot(nvirt,CAM(NA),1,CAM(NA),1) 771! write(lupri,*) 'iocc=', mi, ' sum_a=', work(mi) 772! end do 773 WRITE (LUPRI,*) 774 & 'Important occupied orbital contributions (normalized)' 775 WRITE (LUPRI,*) 776 & 'i, mi, contribution contrib/t1norm**2' 777 778 DO ISYMI = 1,NSYM 779 ISYMA = MULD2H(ISYMAI,ISYMI) 780 DO I = 1,NRHF(ISYMI) 781 MI = IORB(ISYMI) + I 782 NA = IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I-1) + 1 783 work(mi) = 784 & ddot(nvir(isyma),CAM(NA),1,CAM(NA),1) 785 write(lupri,*) i, mi, work(mi), work(mi)/C1NOSQ 786 end do 787 end do 788 !call FINDMAXN(X,IXLEN,IMAX,NMAX) 789C 790 CALL FLSHFO(LUPRI) 791 end if 792C 793C-------------------------------------------- 794C Loop through Double excitation vector. 795C If not ccs or ccp2 796C-------------------------------------------- 797C 798 IF (.NOT. ( CCSEFF .OR. CCP2 )) THEN 799C 800 WRITE(LUPRI,'(A)') 801 * ' +----------------------------------------------' 802 * //'-------------------------------+' 803C 804 805 2 CONTINUE 806 N2 = 0 807C 808 DO 200 ISYMAI = 1,NSYM 809C 810 ISYMBJ = MULD2H(ISYMAI,ISYMTR) 811C 812 DO 210 ISYMJ = 1,NSYM 813C 814 ISYMB = MULD2H(ISYMJ,ISYMBJ) 815C 816 DO 220 ISYMI = 1,NSYM 817C 818 ISYMA = MULD2H(ISYMI,ISYMAI) 819C 820 DO 230 J = 1,NRHF(ISYMJ) 821C 822 MJ = IORB(ISYMJ) + J 823C 824 DO 240 B = 1,NVIR(ISYMB) 825C 826 NBJ = IT1AM(ISYMB,ISYMJ) 827 * + NVIR(ISYMB)*(J - 1) + B 828C 829 MB = IORB(ISYMB) + NRHF(ISYMB) + B 830C 831 DO 250 I = 1,NRHF(ISYMI) 832C 833 MI = IORB(ISYMI) + I 834C 835 DO 260 A = 1,NVIR(ISYMA) 836C 837 NAI = IT1AM(ISYMA,ISYMI) 838 * + NVIR(ISYMA)*(I - 1) + A 839C 840 MA = IORB(ISYMA) + NRHF(ISYMA) + A 841C 842 IF (((ISYMAI.EQ.ISYMBJ).AND. 843 * (NAI .LT. NBJ)).OR.(ISYMAI.LT.ISYMBJ)) 844 * GOTO 260 845C 846 IF (ISYMAI.EQ.ISYMBJ) THEN 847 NAIBJ = IT2AM(ISYMAI,ISYMBJ) 848 * + INDEX(NAI,NBJ) 849 ELSE 850 NAIBJ = IT2AM(ISYMAI,ISYMBJ) 851 * + NT1AM(ISYMAI)*(NBJ-1) + NAI 852 ENDIF 853C 854 KAIBJ = NAIBJ + NT1AM(ISYMTR) 855 IF (ABS(CAM(KAIBJ)) .GT. THR2 ) THEN 856C 857 WRITE(LUPRI,9991) ISYMA,ISYMB,ISYMI,ISYMJ, 858 * A,B,I,J,NAI,NBJ,NAIBJ, 859 * CAM(KAIBJ) 860 N2 = N2 + 1 861C 862 SUMOFP = SUMOFP + CAM(KAIBJ)*CAM(KAIBJ) 863C 864 ENDIF 865C 866 260 CONTINUE 867 250 CONTINUE 868 240 CONTINUE 869 230 CONTINUE 870 220 CONTINUE 871 210 CONTINUE 872 200 CONTINUE 873C 874 IF ((N2.LT.1).AND.(SQRT(C2NOSQ/CNOSQ).GT.1.0D-3)) THEN 875 THR2 = THR2/5D00 876 GOTO 2 877 ENDIF 878C 879 ENDIF 880C 881 WRITE(LUPRI,'(A)') 882 * ' +==============================================' 883 * //'===============================+' 884C 885 WRITE(LUPRI,'(//10X,A,8X,F10.4)') 886 * 'Norm of Printed Amplitude Vector : ',SQRT(SUMOFP) 887 WRITE(LUPRI,'(//10X,A43,1X,F9.6)') 888 * 'Printed all single excitations greater than',THR1 889 IF (.NOT. (CCSEFF.OR.CCP2)) THEN 890 WRITE(LUPRI,'(//10X,A43,1X,F9.6)') 891 * 'Printed all double excitations greater than',THR2 892 ENDIF 893C 894 9990 FORMAT(1X,'| ',I1,3X,I1,2X,' | ',I3,5X,I3,4X,' | ',I8,9x, 895 * ' | ',12x,' | ',1x, F10.6,' |') 896 9991 FORMAT(1X,'| ',I1,1X,I1,1X,I1,1X,I1,' | ', 897 * I3,1X,I3,1X,I3,1X,I3,' | ', 898 * I8,1x,I8,' | ',I12,' | ',1x,F10.6,' |') 899C 900 RETURN 901 END 902C================================================================= 903 SUBROUTINE CC_build_Rfull(LUMAT,FMAT, 904 & Lchain,Lspace,XVEC,EVR, 905 & ER,EI, 906 & WORK,LWORK) 907C-------------------------------------------------------------------- 908C SC,10/03/2011: generate pseudo eigenvectors in full space by: 909C - reading in the Q (or PT) matrix 910! - ddotting with Lanczos eigenvector 911! Lchain is length of Chain 912! Lspace is length of full excitation space 913! Xvec is R (or L) eigenvector in full space - OUTPUT 914! EVR is Lanczos eigenvector - INPUT 915! ER is lanczos eigenvalue (real part), EI is lanczos eigenvalue (im part) 916! LUMAT/FMAT are the unit number and file name of the Q (PT) matrix 917C-------------------------------------------------------------------- 918C 919 IMPLICIT NONE 920#include "priunit.h" 921#include "ccsdsym.h" 922 INTEGER Lchain,Lspace,LWORK,LUMAT 923 CHARACTER FMAT*(*) 924#if defined (SYS_CRAY) 925 REAL XVEC(Lspace*Lchain),EVR(Lchain*Lchain) 926 !REAL XVEC(Lspace,Lchain),EVR(Lchain,Lchain) 927 REAL ER(Lchain), EI(Lchain) 928 REAL DDOT, WORK(LWORK) 929#else 930 !DOUBLE PRECISION XVEC(Lspace,Lchain),EVR(Lchain,Lchain) 931 DOUBLE PRECISION XVEC(Lspace*Lchain),EVR(Lchain*Lchain) 932 DOUBLE PRECISION ER(Lchain), EI(Lchain) 933 DOUBLE PRECISION DDOT, WORK(LWORK) 934 DOUBLE PRECISION ONE,ZERO 935#endif 936 INTEGER kend,lend,isyhop,ivec, n2vec 937 INTEGER KQMAT, IOFF 938 CHARACTER*(10) MODEL 939 PARAMETER (ONE=1.0d0,ZERO=0.0d0) 940 941 942 KQMAT = 1 943 KEND = kQMAT + Lspace*Lchain 944 LEND = LWORK - KEND 945 IF (LEND .LT. 0) THEN 946 WRITE(LUPRI,*) 'Needed:', KEND, 'Available:', LWORK 947 CALL QUIT('Insufficient work space in CC_build_Rfull') 948 END IF 949 950 CALL DZERO(WORK(KQMAT),Lspace*Lchain) 951 952 !Read in the entire Q matrix 953 ioff = 0 954 do k=1,Lchain 955 CALL CC_RVEC(LUMAT,FMAT,Lspace,Lspace,K,WORK(KQMAT+ioff)) 956 ioff = ioff + Lspace 957 enddo 958! write(lupri,*)'---- ' 959! write(lupri,*)' The Q (or P) matrix ' 960! call OUTPUT (WORK(KQMAT),1,Lspace,1,Lchain,Lspace,Lchain, 961! * 1,LUPRI) 962! write(lupri,*)'---- ' 963 964 !Compute R_full = Qmat*R_lanczos_mat (Qmat dim N*J, R_l_mat dim J*J) 965 CALL DGEMM('N','N',Lspace,Lchain,Lchain, 966 * ONE,WORK(KQMAT),Lspace, 967 * EVR,Lchain,ZERO, 968 * XVEC,Lspace) 969 970 write(lupri,*)'Un-normalized R_full matrix computed' 971 call flshfo(lupri) 972 973 RETURN 974 END 975!-------------------------------------------------------------------- 976 SUBROUTINE CC_LANCZOS_SUMRULES(LABEL,J,ER,EI,R,L,F,U,V) 977C-------------------------------------------------------------------- 978C Given a set of Lanzcos T eigenvalues and eigenvector this 979C calculates the sum rules 980C J: Chain Length 981C ER: Eigenvalues, real part 982C EI: Eigenvalues, imag part 983C R: Eigenvectors right (stored as cols) 984C L: Eigenvectors Left (stored as cols, e.i. eigvectors of T^T) 985C F: F-matrix "trans" (not double q-transformed) 986C U, V are normalization factors 987C-------------------------------------------------------------------- 988 IMPLICIT NONE 989#include "priunit.h" 990#include "cclrlancz.h" 991#include "codata.h" 992 993 INTEGER J,IFREQ,NFREQ 994 DOUBLE PRECISION U,V 995 DOUBLE PRECISION ER(J),EI(J),R(J,J),L(J,J),F(J,J) 996 DOUBLE PRECISION ZERO, ONE, TWO, PSUMK 997 DOUBLE PRECISION sum_s0,sum_s0_f,sum_s0_t(-6:2) 998 DOUBLE PRECISION sum_l0,sum_l0_f,sum_l0_t(-6:2) 999 DOUBLE PRECISION mean_exc_energy(-6:2) 1000 DOUBLE PRECISION FACTOR,FACTOR2 1001 1002 LOGICAL NOIMAG 1003 INTEGER I, K, IDAMP, N 1004! 1005 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 1006 CHARACTER*(8) LABEL 1007! 1008 LOGICAL SKIPNEXT 1009 SKIPNEXT=.FALSE. 1010C 1011C Loop over chain 1012C 1013 WRITE (LUPRI,*) " FOR OPERATOR LABEL = ", LABEL 1014 1015 !WRITE (LUPRI,*) 'check FREQ=', FREQ 1016 CALL DZERO(SUM_S0_T,9) 1017 CALL DZERO(SUM_L0_T,9) 1018 CALL DZERO(mean_exc_energy,9) 1019 1020 DO N=-6,2 1021 SUM_S0=ZERO 1022 SUM_L0=ZERO 1023 SUM_S0_F=ZERO 1024 SUM_L0_F=ZERO 1025 DO I=1,J 1026 IF (SKIPNEXT) THEN 1027 SKIPNEXT = .FALSE. 1028 ELSE 1029 !check on imaginary eigenvalue = 0 1030 IF (EI(I).EQ.ZERO) THEN 1031 FACTOR=U*V*R(1,I)*L(1,I) 1032C No imaginary eigenvector, calc eta*t conts 1033 SUM_S0 = SUM_S0 + FACTOR*ER(I)**(N+1) 1034 SUM_L0 = SUM_L0 + FACTOR*(ER(I)**(N+1))*DLOG(ER(I)) 1035C 1036C Include F-part cont 1037C 1038 FACTOR = V*V*L(1,I)*ER(I)**(N+1) 1039 DO K=1,J 1040 FACTOR2 = FACTOR*L(1,K)*F(I,K) 1041 if (abs(ei(k)).eq.zero) then 1042 !denominator 1043 PSUMK = ER(I)+ER(K) 1044 SUM_S0_F = SUM_S0_F - FACTOR2/PSUMK 1045 SUM_L0_F = SUM_L0_F - FACTOR2*DLOG(ER(I))/PSUMK 1046 else 1047 write(lupri,*)'LANCZOS, do not know, index=',k 1048 end if 1049 ENDDO 1050 1051 ELSE 1052 SKIPNEXT=.true. 1053 ENDIF 1054 end if !skipnext 1055C 1056 ENDDO 1057 1058 SUM_S0_T(N) = (SUM_S0 + SUM_s0_F)*TWO 1059 SUM_L0_T(N) = (SUM_L0 + SUM_l0_F)*TWO 1060 END DO 1061 DO N=-6,2 1062 if (SUM_S0_T(N).eq.zero) then 1063 mean_exc_energy(N) = zero 1064 else 1065 mean_exc_energy(N) = xtev*dexp(SUM_L0_T(N)/SUM_S0_T(N)) 1066 end if 1067 END DO 1068 1069 CALL HEADER('Oscillator strength sum rules',30) 1070 WRITE (LUPRI,'(//,14X,A,/,6X,A,5X,A,A,A,/)') 1071 & 'S(N) Sum Rules : Dipole Length Approximation in a.u.', 1072 & 'N',LABEL(1:1),LABEL(1:1),' - component' 1073 WRITE (LUPRI,'(9(5X,I3,(4X,G13.6)/))') 1074 & (N,SUM_S0_T(N),N=-6,2) 1075 WRITE (LUPRI,'(//,14X,A,/,6X,A,5X,A,A,A,/)') 1076 & 'L(N) Sum Rules : Dipole Length Approximation in a.u.', 1077 & 'N',LABEL(1:1),LABEL(1:1),' - component' 1078 WRITE (LUPRI,'(9(5X,I3,(4X,G13.6)/))') 1079 & (N,SUM_L0_T(N),N=-6,2) 1080 WRITE (LUPRI,'(//,14X,A,/,6X,A,5X,A,A,A,/)') 1081 & 'I(N) Sum Rules : Dipole Length Approximation in eV ', 1082 & 'N',LABEL(1:1),LABEL(1:1),' - component' 1083 WRITE (LUPRI,'(9(5X,I3,(4X,G13.6)/))') 1084 & (N,mean_exc_energy(N),N=-6,2) 1085 1086 RETURN 1087 END 1088!-------------------------------------------------------------------- 1089 SUBROUTINE CC_LANCZOS_OSCSTR(LABEL,J,ER,EI,R,L,F,U,V) 1090C-------------------------------------------------------------------- 1091C Given a set of Lanzcos T eigenvalues and eigenvector this 1092C calculates the oscillator strengths 1093C J: Chain Length 1094C ER: Eigenvalues, real part 1095C EI: Eigenvalues, imag part 1096C R: Eigenvectors right (stored as cols) 1097C L: Eigenvectors Left (stored as cols, e.i. eigvectors of T^T) 1098C F: F-matrix "trans" (not double q-transformed) 1099C U, V are normalization factors 1100C-------------------------------------------------------------------- 1101 IMPLICIT NONE 1102#include "priunit.h" 1103#include "cclrlancz.h" 1104#include "codata.h" 1105 1106 INTEGER J,IFREQ,NFREQ 1107 DOUBLE PRECISION U,V 1108 DOUBLE PRECISION ER(J),EI(J),R(J,J),L(J,J),F(J,J) 1109 DOUBLE PRECISION F0I(J) 1110 DOUBLE PRECISION ZERO, ONE, TWO, PSUMK, THREE 1111 DOUBLE PRECISION sum_s0,sum_s0_f 1112 DOUBLE PRECISION FACTOR,FACTOR2 1113 1114 LOGICAL NOIMAG 1115 INTEGER I, K, IDAMP, N 1116! 1117 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,THREE=3.0d0) 1118 CHARACTER*7 FNEXCIL 1119 PARAMETER (FNEXCIL='CCEXCIL') 1120 INTEGER IDUMMY, LUEXCIL 1121 1122 CHARACTER*(8) LABEL 1123! 1124 LOGICAL SKIPNEXT 1125 SKIPNEXT=.FALSE. 1126C 1127C Loop over chain 1128C 1129 WRITE (LUPRI,*) " FOR OPERATOR LABEL = ", LABEL 1130 LUEXCIL=-1 1131 CALL GPOPEN(LUEXCIL,FNEXCIL,'NEW',' ','FORMATTED', 1132 & IDUMMY,.FALSE.) 1133! 1134 CALL DZERO(F0I,J) 1135 1136 DO I=1,J 1137 IF (SKIPNEXT) THEN 1138 SKIPNEXT = .FALSE. 1139 ELSE 1140 !check on imaginary eigenvalue = 0 1141 IF (EI(I).EQ.ZERO) THEN 1142 FACTOR=U*V*R(1,I)*L(1,I) 1143C No imaginary eigenvector, calc eta*t conts 1144 F0I(I) = FACTOR*ER(I) 1145C Include F-part cont 1146 FACTOR = V*V*L(1,I)*ER(I) 1147 SUM_S0_F = ZERO 1148 DO K=1,J 1149 FACTOR2 = FACTOR*L(1,K)*F(I,K) 1150 if (ei(k).eq.zero) then 1151 !denominator 1152 PSUMK = ER(I)+ER(K) 1153 SUM_S0_F = SUM_S0_F - FACTOR2/PSUMK 1154 else 1155 write(lupri,*)'LANCZOS, do not know, index=',k 1156 end if 1157 ENDDO 1158 F0I(I) = (F0I(I) + SUM_S0_F)*TWO/THREE 1159 1160 ELSE 1161 SKIPNEXT=.true. 1162 ENDIF 1163 end if !skipnext 1164C 1165 ENDDO 1166 1167 !CALL HEADER('Oscillator strengths',30) 1168 WRITE (LUEXCIL,'(14X,A,/,6X,A,8X,A,A,A)') 1169 & 'Oscillator strengths: Dipole Length Approximation in a.u.', 1170 & 'I',LABEL(1:1),LABEL(1:1),' - component' 1171! WRITE (LUPRI,'(J(5X,I3,(4X,G13.6)/))') 1172! & (I,ER(I),F0I(I),I=1,J) 1173 WRITE (LUEXCIL,*) 'T MATRIX EIGENVALUES (eV) FOR JCHAIN = ',J 1174 WRITE (LUEXCIL,*) 1175 & 'INDEX & REAL (eV) & IMAG (eV) & REAL (au) & IMAG (au) & OSCST' 1176 do i=1,j 1177 write(luexcil,'(X,I4,5F16.8)') i, ER(i)*XTEV, 1178 & EI(I)*XTEV, ER(I), EI(I), F0I(i) 1179 end do 1180 1181 CALL GPCLOSE(LUEXCIL,'KEEP') 1182 RETURN 1183 END 1184!---------------- 1185C================================================================= 1186 SUBROUTINE CC_build_Rfull2(LUMAT,FMAT, 1187 & nexci,Lchain,Lspace, 1188 & XVEC,EVR,ER,EI, 1189 & WORK,LWORK) 1190C-------------------------------------------------------------------- 1191C SC,05/2011: generate pseudo eigenvectors in full space by: 1192C - reading in the Q (or PT) matrix 1193! - multiplying with Lanczos eigenvector 1194! Nexci is the nr of vectors in full space that will be generated 1195! Lchain is length of Chain 1196! Lspace is length of full excitation space 1197! Xvec is R (or L) eigenvector matrix in full space - OUTPUT 1198! EVR is Lanczos eigenvectors - INPUT 1199! ER is lanczos eigenvalue (real part), EI is lanczos eigenvalue (im part) 1200! LUMAT/FMAT are the unit number and file name of the Q (PT) matrix 1201C-------------------------------------------------------------------- 1202 IMPLICIT NONE 1203#include "priunit.h" 1204#include "ccsdsym.h" 1205 INTEGER Lchain,Lspace,LWORK,LUMAT,nexci 1206 CHARACTER FMAT*(*) 1207#if defined (SYS_CRAY) 1208 REAL XVEC(Lspace*nexci),EVR(Lchain*nexci) 1209 REAL ER(nexci), EI(nexci) 1210 REAL DDOT, WORK(LWORK) 1211#else 1212 DOUBLE PRECISION XVEC(Lspace*nexci),EVR(Lchain*nexci) 1213 DOUBLE PRECISION ER(nexci), EI(nexci) 1214 DOUBLE PRECISION DDOT, WORK(LWORK) 1215 DOUBLE PRECISION ONE,ZERO 1216#endif 1217 INTEGER kend,lend,isyhop,ivec, n2vec 1218 INTEGER KQMAT, IOFF 1219 CHARACTER*(10) MODEL 1220 PARAMETER (ONE=1.0d0,ZERO=0.0d0) 1221 KQMAT = 1 1222 KEND = kQMAT + Lspace*Lchain 1223 LEND = LWORK - KEND 1224 IF (LEND .LT. 0) THEN 1225 WRITE(LUPRI,*) 'Needed:', KEND, 'Available:', LWORK 1226 CALL QUIT('Insufficient work space in CC_build_Rfull') 1227 END IF 1228 1229 CALL DZERO(WORK(KQMAT),Lspace*Lchain) 1230 1231 !Read in the entire Q matrix 1232 ioff = 0 1233 do k=1,Lchain 1234 CALL CC_RVEC(LUMAT,FMAT,Lspace,Lspace,K,WORK(KQMAT+ioff)) 1235 ioff = ioff + Lspace 1236 enddo 1237! write(lupri,*)'---- ' 1238! write(lupri,*)' The Q (or P) matrix ' 1239! call OUTPUT (WORK(KQMAT),1,Lspace,1,Lchain,Lspace,Lchain, 1240! * 1,LUPRI) 1241! 1242! write(lupri,*)'---- ' 1243 1244 !Compute R_full = Qmat*R_lanczos_mat (Qmat dim N*J, R_l_mat dim J*nexci) 1245 CALL DGEMM('N','N',Lspace,nexci,Lchain, 1246 * ONE,WORK(KQMAT),Lspace, 1247 * EVR,Lchain,ZERO, 1248 * XVEC,Lspace) 1249 1250 write(lupri,*)'Un-normalized R/L_full matrix block computed' 1251 call flshfo(lupri) 1252 1253 RETURN 1254 END 1255 1256C================================================================= 1257 SUBROUTINE CC_check_resid(ECURR,iside,nexc,Lspace,isyhop, 1258 & XVEC,ER,EI, 1259 & WORK,LWORK,APROXR12,N2VEC, 1260 & tstomega) 1261C-------------------------------------------------------------------- 1262C Check the residual on the pseudo eigenvectors in full space by: 1263C - Performing AR_i-w_iR 1264! - computing norm of resulting vector 1265! Nexci is the nr of vectors in full space that will be generated 1266! Lspace is length of full excitation space 1267! Xvec is the matrix of R (or L) eigenvectors 1268! ER is lanczos eigenvalue (real part), EI is lanczos eigenvalue (im part) 1269C-------------------------------------------------------------------- 1270 IMPLICIT NONE 1271#include "priunit.h" 1272#include "ccsdsym.h" 1273#include "dummy.h" 1274 INTEGER Lspace,LWORK,nexc,iside 1275 DOUBLE PRECISION XVEC(Lspace,nexc) 1276 DOUBLE PRECISION ER(nexc), EI(nexc) 1277 DOUBLE PRECISION DDOT, WORK(LWORK) 1278 DOUBLE PRECISION ONE,ZERO,xnorm,ECURR 1279 DOUBLE PRECISION omega1,Rdag_R,Rdag_A_R,xnorm2 1280 1281 INTEGER katrr,kend,lend,isyhop,ivec,n2vec 1282 INTEGER IOFF 1283 logical tstomega 1284 CHARACTER*(10) MODEL 1285 CHARACTER*(3) APROXR12 1286 PARAMETER (ONE=1.0d0,ZERO=0.0d0) 1287 1288 KATRR = 1 1289 KEND = KATRR + Lspace 1290 LEND = LWORK - KEND 1291 IF (LEND .LT. 0) THEN 1292 WRITE(LUPRI,*) 'Needed:', KEND, 'Available:', LWORK 1293 CALL QUIT('Insufficient work space in CC_check_resid') 1294 END IF 1295 1296 !Transform one excitation vector at a time. 1297 !ioff = 0 1298 do k=1,nexc 1299 call dcopy(Lspace,XVEC(1,k),1,work(katrr),1) 1300 CALL cc_atrr(ECURR,ISYHOP,ISIDE,WORK(Katrr),Lend,.false., 1301 & dummy,APROXR12,.false.) 1302! SUBROUTINE CC_ATRR(ECURR,ISYMV,ISIDE,WORK,LWORK,LCONT,CONT, 1303! & APROXR12) 1304 1305 !WRITE(LUPRI,*) 'AR vector' 1306 !CALL CC_PRP(WORK(KAtrr),WORK(KAtrr+nt1amx),1,1,N2VEC) 1307 !call dcopy(Lspace,XVEC(1,k),1,work(kend),1) 1308 !call dscal(Lspace,ER(k),work(kend),1) 1309 !WRITE(LUPRI,*) 'wR vector' 1310 !CALL CC_PRP(WORK(Kend),WORK(Kend+nt1amx),1,1,N2VEC) 1311 !CALL DAXPY(Lspace,-1.0d0,WORK(Kend),1,WORK(Katrr),1) 1312 !WRITE(LUPRI,*) 'AR-wR vector' 1313 !CALL CC_PRP(WORK(KAtrr),WORK(KAtrr+nt1amx),1,1,N2VEC) 1314 CALL DAXPY(Lspace,-ER(k),XVEC(1,k),1,WORK(Katrr),1) 1315 xnorm=dsqrt(ddot(Lspace,WORK(KAtrr),1,WORK(KAtrr),1)) 1316 write(lupri,*)'Norm residual eival=', er(k), ' is =', xnorm 1317 !ioff = ioff + Lspace 1318 enddo 1319 ! 1320 ! Now try and compute omega'=<R|A|R>/<R|R> and Res=AR-w'R 1321 ! 1322 if (tstomega) then 1323 do k=1,nexc 1324 call dcopy(Lspace,XVEC(1,k),1,work(katrr),1) 1325 CALL cc_atrr(ECURR,ISYHOP,ISIDE,WORK(Katrr),Lend,.false., 1326 & dummy,APROXR12,.false.) 1327 Rdag_A_R=ddot(Lspace,XVEC(1,k),1,WORK(KAtrr),1) 1328 Rdag_R=ddot(Lspace,XVEC(1,k),1,XVEC(1,k),1) 1329 omega1=Rdag_A_R/Rdag_R 1330 !WRITE(LUPRI,*) 'AR vector' 1331 !CALL CC_PRP(WORK(KAtrr),WORK(KAtrr+nt1amx),1,1,N2VEC) 1332 !call dcopy(Lspace,XVEC(1,k),1,work(kend),1) 1333 !call dscal(Lspace,ER(k),work(kend),1) 1334 !WRITE(LUPRI,*) 'wR vector' 1335 !CALL CC_PRP(WORK(Kend),WORK(Kend+nt1amx),1,1,N2VEC) 1336 !CALL DAXPY(Lspace,-1.0d0,WORK(Kend),1,WORK(Katrr),1) 1337 !WRITE(LUPRI,*) 'AR-wR vector' 1338 !CALL CC_PRP(WORK(KAtrr),WORK(KAtrr+nt1amx),1,1,N2VEC) 1339 CALL DAXPY(Lspace,-omega1,XVEC(1,k),1,WORK(Katrr),1) 1340 xnorm2=dsqrt(ddot(Lspace,WORK(KAtrr),1,WORK(KAtrr),1)) 1341 write(lupri,*)'Norm residual omega1=', omega1,' is =',xnorm2 1342 write(lupri,*)'Diff |omega1-ER(k)|=', abs(omega1-ER(k)) 1343 !ioff = ioff + Lspace 1344 enddo 1345 end if 1346 RETURN 1347 END 1348 1349C /* Deck cc_normalize */ 1350 SUBROUTINE cc_normalize(LENGTH,N,ZR,ZL,ER,EI,ISYHOP,IOPT) 1351C 1352C normalize pseudo RE 1353C 1354#include "implicit.h" 1355#include "priunit.h" 1356 DIMENSION ZR(LENGTH,N), ZL(LENGTH,N) 1357 DIMENSION ER(N), EI(N) 1358 PARAMETER ( D0 = 0.0D0, D1 = 1.0D0 ) 1359C 1360C Normalize RE eigenvectors 1361C 1362 if (iopt.eq.1) then 1363 DO I = 1,N 1364 ZNRM = DDOT(LENGTH,ZR(1,I),1,ZR(1,I),1) 1365 !write(lupri,*)'Norm RrRr', ZNRM 1366 ZNRM = D1 / SQRT(ZNRM) 1367 !IMAX = IDAMAX(N,ZR(1,I),1) 1368 !IF (ZR(IMAX,I) .LT. D0) ZNRM = -ZNRM 1369 CALL DSCAL(LENGTH,ZNRM,ZR(1,I),1) 1370 ! 1371! WRITE (LUPRI,'(/A,I4,A,1P,D15.6/A)') 1372! * ' - XR EIGENVALUE NO.',I, ' - EIGENVALUE',ER(I), 1373! * ' - XR EIGENVECTOR :' 1374! CALL CC_PRAM(ZR(1,I),PT1,ISYHOP,.FALSE.) 1375 ! 1376 END DO 1377 end if 1378C 1379C binormalize LE eigenvectors 1380C 1381 if (iopt.eq.2) then 1382 DO I = 1,N 1383 ZNRM = DDOT(LENGTH,ZL(1,I),1,ZR(1,I),1) 1384 write(lupri,*)'Norm LrRr of vector nr:',I,' is=', ZNRM 1385 ZNRM = D1 / (ZNRM) 1386 !IMAX = IDAMAX(N,ZR(1,I),1) 1387 !IF (ZR(IMAX,I) .LT. D0) ZNRM = -ZNRM 1388 CALL DSCAL(LENGTH,ZNRM,ZL(1,I),1) 1389 ! 1390 ZNRM = DDOT(LENGTH,ZL(1,I),1,ZR(1,I),1) 1391 !write(lupri,*)'Check Norm LrRr of vector nr:',I,' is=', ZNRM 1392! WRITE (LUPRI,'(/A,I4,A,1P,D15.6/A)') 1393! & '- LPT EIGENVALUE NO.',I,' - EIGENVALUE', 1394! & ER(I), 1395! & '- LPT EIGENVECTOR :' 1396! CALL CC_PRAM(ZL(1,I),PT1,ISYHOP,.FALSE.) 1397 1398 END DO 1399 end if 1400 RETURN 1401 END 1402C================================================================= 1403 SUBROUTINE CC_dump_onfile(length,nex,XVEC,isyhop,lista, 1404 & work,lwork,model,ioptrw) 1405C-------------------------------------------------------------------- 1406! Dump on file, according to lista type 1407C-------------------------------------------------------------------- 1408 IMPLICIT NONE 1409#include "priunit.h" 1410#include "ccsdsym.h" 1411#include "dummy.h" 1412 INTEGER length,nex,isyhop,lwork,ioptrw 1413#if defined (SYS_CRAY) 1414 REAL XVEC(length,nex),work(lwork) 1415 REAL DDOT 1416#else 1417 DOUBLE PRECISION XVEC(length,nex), work(lwork) 1418 DOUBLE PRECISION DDOT 1419 DOUBLE PRECISION ONE,ZERO 1420#endif 1421 INTEGER kend,lend,ivec,n2vec 1422 INTEGER IOFF 1423 CHARACTER*(10) MODEL 1424 CHARACTER*(3) LISTA 1425 PARAMETER (ONE=1.0d0,ZERO=0.0d0) 1426 1427 !Dump 1428 kend=1 1429 lend=lwork-kend 1430 IF (LEND .LT. 0) THEN 1431 WRITE(LUPRI,*) 'Needed:', KEND, 'Available:', LWORK 1432 CALL QUIT('Insufficient work space in CC_dump_onfile') 1433 END IF 1434 do I=1,nex 1435 CALL CC_WRRSP(lista,I,ISYHOP,IOPTRW,MODEL,WORK(KEND), 1436 * XVEC(1,I),XVEC(1+nt1am(isyhop),I), 1437 * WORK(kend),Lend) 1438 enddo 1439 RETURN 1440 END 1441C /* Deck cc_resid_lanczos */ 1442 SUBROUTINE cc_resid_lanczos(beta,LENGTHQ,N_RE,Qi,ZR,ER,EI) 1443C 1444C Check residual according to Lanczos recipe 1445C N_RE is the number of Lanczos eigenvectors, each of them long N_RE 1446C (squared matrix). We use the last element of each Lanczos eigenvector. 1447C 1448#include "implicit.h" 1449#include "priunit.h" 1450 DIMENSION ZR(N_RE,N_RE), Qi(lengthq) 1451 DIMENSION ER(N_RE), EI(N_RE) 1452 PARAMETER ( D0 = 0.0D0, D1 = 1.0D0 ) 1453C 1454 xqnorm = ddot(lengthq,Qi,1,Qi,1) 1455 xqnorm = sqrt(xqnorm) 1456 write(lupri,*)'norm of Q_i+1 is ', xqnorm 1457 write(lupri,*)'beta_jchain=', beta 1458 1459 1460 DO I = 1,N_RE 1461 residREi = beta*xqnorm*abs(ZR(N_RE,I)) 1462 write(lupri,*)'Residual norm Lanczos wi:',ER(i),' =', 1463 & residREi 1464 END DO 1465C 1466 RETURN 1467 END 1468 1469!---------------------- 1470C /* Deck ccbiort */ 1471 SUBROUTINE CCBIORT(LUQMAT,FQMAT,LUPMAT,FPMAT,ichain, 1472 & NPREV,NCCVAR, 1473 & QNEW,PNEW,ISYMTR,THRLDP,WRK,LWRK) 1474C 1475C Modified version of CCORT to biorthogonalize p^t and q 1476C vectors of Lanczos chain 1477C two-sided modified Gram-Schmidt (TSMGS) process 1478C Notice we assume to enter the routine with just 1 Q_new and 1 P_new 1479C 1480C Purpose: 1481C Bi-orthogonalize new p and q vectors against all previous p 1482C and q vectors and among themselves, and renormalize. 1483C (Orthogonalization is performed twice if round-off is large, 1484C if larger than THRRND). 1485C 1486C Input: 1487C LUQMAT - file with previous Q vectors with name FQMAT 1488C LUPMAT - file with previous P vectors with name FPMAT 1489C QNEW, PNEW - non-orthogonal q and p new vectors 1490C NBPREV - number of previous q and p vectors on FQMAT/FPMAT 1491C THRLDP - threshold for linear dependence 1492C 1493C Output: 1494C QNEW and PNEW on output are biorthogonal to previous ones 1495C 1496C Scratch: 1497C require WRK of length NCCVAR 1498C 1499#include "implicit.h" 1500#include "priunit.h" 1501#include "ccsdinp.h" 1502#include "ccsdsym.h" 1503C 1504 CHARACTER FQMAT*(*),FPMAT*(*) 1505 DIMENSION QNEW(NCCVAR),PNEW(NCCVAR), WRK(NCCVAR) 1506 LOGICAL LOCDBG 1507 PARAMETER (THRRND=1.D-4, THRTT=1.D-6, THRZER=1.0D-35) 1508 PARAMETER (THRTST=1.0D-12) 1509 PARAMETER (D1 = 1.0D0) 1510 PARAMETER (LOCDBG=.false.) 1511 1512 INTEGER SVEC, RVEC, LVEC 1513C 1514C 1515 IF (LOCDBG) THEN 1516 WRITE (LUPRI,'(//A/)') ' --- Debug output from CCORT ---' 1517 END IF 1518 KEND1 = 1 1519 LWRK1 = LWRK - KEND1 1520 IF (LWRK1.LT.0) THEN 1521 WRITE(LUPRI,*) 'LWRK, KEND1: ',LWRK, KEND1 1522 CALL QUIT('Insufficient memory in CCORT') 1523 END IF 1524C 1525C Check bi-norm of input Q and P vectors 1526C Issue warning if norm .le. THRZER 1527C 1528 TTPQ = DDOT(NCCVAR,PNEW(1),1,QNEW(1),1) 1529! WRITE (LUPRI,*) ' Initial overlap <Pnew|Qnew>:' 1530! & ,TTPQ, 'for ichain=', ichain 1531 !THRZER=1.0D-35 1532 IF (TTPQ.LT.THRZER) THEN 1533 WRITE(LUPRI,*) 'CCBIORTHO: WARNING!! <pi|qi> <= 0' 1534 END IF 1535C 1536C work space allocation 1537C 1538 KQPRE = KEND1 1539 KPPRE = KQPRE + NCCVAR 1540 KEND1 = KPPRE + NCCVAR 1541 IROUND = 0 1542 ITURN = 0 1543 1544 1500 CONTINUE 1545 ITURN = ITURN + 1 1546C 1547C Orthogonalize new q and p vectors against previous ones 1548C 1549 DO K = 1, NPREV 1550 CALL CC_RVEC(LUQMAT,FQMAT,nccvar, 1551 * nccvar,K,WRK(KQPRE)) 1552 CALL CC_RVEC(LUPMAT,FPMAT,nccvar, 1553 * nccvar,K,WRK(KPPRE)) 1554 1555 !<Pold|Qnew> 1556 TTQ = DDOT(NCCVAR,WRK(KPPRE),1,QNEW(1),1) 1557 !<Pnew|Qold> 1558 TTP = DDOT(NCCVAR,PNEW(1),1,WRK(KQPRE),1) 1559 !<Pold|Qold> 1560 TTDQP = DDOT(NCCVAR,WRK(KPPRE),1,WRK(KQPRE),1) 1561 IF (LOCDBG) THEN 1562 WRITE (LUPRI,*)'start <P(',K,')|Qnew(',ichain,')> =', TTQ 1563 WRITE (LUPRI,*)'start <Pnew(',ichain,')| Q(',K,')> =', TTP 1564 WRITE (LUPRI,*)'start <P(',K,')|Q(',K,')>=', TTDQP 1565 END IF 1566 !Qnew = |Qnew> - |Qold>*<Pold|Qnew> 1567 CALL DAXPY(NCCVAR,-TTQ,WRK(KQPRE),1,QNEW(1),1) 1568 !<Pnew| = <Pnew| - <Pnew|Qold><Pold| 1569 CALL DAXPY(NCCVAR,-TTP,WRK(KPPRE),1,PNEW(1),1) 1570C ortho test 1571 TTQ = DDOT(NCCVAR,WRK(KPPRE),1,QNEW(1),1) 1572 TTP = DDOT(NCCVAR,PNEW(1),1,WRK(KQPRE),1) 1573 TTDQP = DDOT(NCCVAR,WRK(KPPRE),1,WRK(KQPRE),1) 1574 if (LOCDBG) then 1575 WRITE(LUPRI,*)'BIORT ortho: <P(',k,')|Qnew(',ichain,')>=', TTQ 1576 WRITE(LUPRI,*)'BIORT ortho: <Pnew(',ichain,')|Q(',k,')>=', TTP 1577 WRITE(LUPRI,*)'Biort ortho: <P(',K,')|Q(',K,')>=', TTDQP 1578 end if 1579 END DO 1580C 1581C Bi-normalize new vectors 1582C (and maybe check if vectors are linear dependent?) 1583C 1584 TTN = DDOT(NCCVAR,PNEW(1),1,QNEW(1),1) 1585 WRITE(LUPRI,*) 'BIORT check ortho: <Pnew|Qnew>(',ichain,')=', TTN 1586 TTqnn = DDOT(NCCVAR,QNEW(1),1,QNEW(1),1) 1587 WRITE(LUPRI,*) 'BIORT check ortho: <Qnew|Qnew>(',ichain,')=',TTqnn 1588! IF (abs(TTN) .LT. THRRND) IROUND = IROUND+1 1589! SQRTTN = SQRT(TTN) 1590! IF (SQRTTN.LT.THRTT) THEN 1591! CALL DSCAL(NCCVAR,(D1/THRTT),PNEW(1),1) 1592! TT = DDOT(NCCVAR,PNEW(1),1,QNEW(1),1) 1593! SQRTT = SQRT(TT) 1594! END IF 1595 CALL DSCAL(NCCVAR,(D1/TTN),PNEW(1),1) 1596 !factor=sign(ttn)*(D1/sqrt((abs(ttn))) 1597 !CALL DSCAL(NCCVAR,factor,PNEW(1),1) 1598 !CALL DSCAL(NCCVAR,factor,QNEW(1),1) 1599C 1600 !IF (IROUND.GT.0 .AND. ITURN.EQ.1 ) GO TO 1500 1601 IF (ITURN.EQ.1 ) GO TO 1500 1602 1603C 1604C *** End of subroutine CCBiORT 1605C 1606 IF (LOCDBG) THEN 1607 WRITE (LUPRI,'(/A//)')' --- End debug output from CCBiORT ---' 1608 END IF 1609 1610 RETURN 1611 END 1612*============================================= 1613C /* Deck ccbiortPQ */ 1614 SUBROUTINE CCBIORTPQ(LUQMAT,FQMAT,LUPMAT,FPMAT,ichain, 1615 & NPREV,NCCVAR, 1616 & QNEW,PNEW,ISYMTR,THRLDP,WRK,LWRK) 1617C 1618C Modified version of CCORT to biorthogonalize p^t and q 1619C vectors of Lanczos chain 1620C two-sided modified Gram-Schmidt (TSMGS) process 1621C Notice we assume to enter the routine with just 1 Q_new and 1 P_new 1622C NORMALIZE BY SHARING PQ NORM ON BOTH P AND Q 1623C 1624C Purpose: 1625C Bi-orthogonalize new p and q vectors against all previous p 1626C and q vectors and among themselves, and renormalize. 1627C (Orthogonalization is performed twice if round-off is large, 1628C if larger than THRRND). 1629C 1630C Input: 1631C LUQMAT - file with previous Q vectors with name FQMAT 1632C LUPMAT - file with previous P vectors with name FPMAT 1633C QNEW, PNEW - non-orthogonal q and p new vectors 1634C NBPREV - number of previous q and p vectors on FQMAT/FPMAT 1635C THRLDP - threshold for linear dependence 1636C 1637C Output: 1638C QNEW and PNEW on output are biorthogonal to previous ones 1639C 1640C Scratch: 1641C require WRK of length NCCVAR 1642C 1643#include "implicit.h" 1644#include "priunit.h" 1645#include "ccsdinp.h" 1646#include "ccsdsym.h" 1647C 1648 CHARACTER FQMAT*(*),FPMAT*(*) 1649 DIMENSION QNEW(NCCVAR),PNEW(NCCVAR), WRK(NCCVAR) 1650 PARAMETER (THRRND=1.D-4, THRTT=1.D-6, THRZER=1.0D-35) 1651 PARAMETER (THRTST=1.0D-12) 1652 PARAMETER (D1 = 1.0D0) 1653 LOGICAL LOCDBG 1654 PARAMETER (LOCDBG = .false.) 1655 1656 INTEGER SVEC, RVEC, LVEC 1657C 1658C 1659 IF (LOCDBG) THEN 1660 WRITE (LUPRI,'(//A/)') ' --- Debug output from CCORT ---' 1661 END IF 1662 KEND1 = 1 1663 LWRK1 = LWRK - KEND1 1664 IF (LWRK1.LT.0) THEN 1665 WRITE(LUPRI,*) 'LWRK, KEND1: ',LWRK, KEND1 1666 CALL QUIT('Insufficient memory in CCORT') 1667 END IF 1668C 1669C Check bi-norm of input Q and P vectors 1670C Issue warning if norm .le. THRZER 1671C 1672 TTPQ = DDOT(NCCVAR,PNEW(1),1,QNEW(1),1) 1673! WRITE (LUPRI,*) ' Initial overlap <Pnew|Qnew>:' 1674! & ,TTPQ, 'for ichain=', ichain 1675 !THRZER=1.0D-35 1676 IF (TTPQ.LT.THRZER) THEN 1677 WRITE(LUPRI,*) 'CCBIORTHO: WARNING!! <pi|qi> <= 0' 1678 END IF 1679C 1680C work space allocation 1681C 1682 KQPRE = KEND1 1683 KPPRE = KQPRE + NCCVAR 1684 KEND1 = KPPRE + NCCVAR 1685 IROUND = 0 1686 ITURN = 0 1687 1688 1500 CONTINUE 1689 ITURN = ITURN + 1 1690C 1691C Orthogonalize new q and p vectors against previous ones 1692C 1693 DO K = 1, NPREV 1694 CALL CC_RVEC(LUQMAT,FQMAT,nccvar, 1695 * nccvar,K,WRK(KQPRE)) 1696 CALL CC_RVEC(LUPMAT,FPMAT,nccvar, 1697 * nccvar,K,WRK(KPPRE)) 1698 1699 !<Pold|Qnew> 1700 TTQ = DDOT(NCCVAR,WRK(KPPRE),1,QNEW(1),1) 1701 !<Pnew|Qold> 1702 TTP = DDOT(NCCVAR,PNEW(1),1,WRK(KQPRE),1) 1703 !<Pold|Qold> 1704 TTDQP = DDOT(NCCVAR,WRK(KPPRE),1,WRK(KQPRE),1) 1705 IF (LOCDBG) THEN 1706 WRITE (LUPRI,*)'start <P(',K,')|Qnew(',ichain,')> =', TTQ 1707 WRITE (LUPRI,*)'start <Pnew(',ichain,')| Q(',K,')> =', TTP 1708 WRITE (LUPRI,*)'start <P(',K,')|Q(',K,')>=', TTDQP 1709 END IF 1710 !Qnew = |Qnew> - |Qold>*<Pold|Qnew> 1711 CALL DAXPY(NCCVAR,-TTQ,WRK(KQPRE),1,QNEW(1),1) 1712 !<Pnew| = <Pnew| - <Pnew|Qold><Pold| 1713 CALL DAXPY(NCCVAR,-TTP,WRK(KPPRE),1,PNEW(1),1) 1714C ortho test 1715 TTQ = DDOT(NCCVAR,WRK(KPPRE),1,QNEW(1),1) 1716 TTP = DDOT(NCCVAR,PNEW(1),1,WRK(KQPRE),1) 1717 TTDQP = DDOT(NCCVAR,WRK(KPPRE),1,WRK(KQPRE),1) 1718 if (LOCDBG) then 1719 WRITE(LUPRI,*)'BIORT ortho: <P(',k,')|Qnew(',ichain,')>=', TTQ 1720 WRITE(LUPRI,*)'BIORT ortho: <Pnew(',ichain,')|Q(',k,')>=', TTP 1721 WRITE(LUPRI,*)'Biort ortho: <P(',K,')|Q(',K,')>=', TTDQP 1722 end if 1723 END DO 1724C 1725C Bi-normalize new vectors 1726C (and maybe check if vectors are linear dependent?) 1727C 1728 TTN = DDOT(NCCVAR,PNEW(1),1,QNEW(1),1) 1729 WRITE(LUPRI,*) 'BiORTPQ check <Pnew|Qnew>(',ichain,')=', TTN 1730 TTqnn = DDOT(NCCVAR,QNEW(1),1,QNEW(1),1) 1731 WRITE(LUPRI,*) 'BiORTPQ check <Qnew|Qnew>(',ichain,')=',TTqnn 1732! IF (abs(TTN) .LT. THRRND) IROUND = IROUND+1 1733! SQRTTN = SQRT(ABS(TTN)) 1734! IF (SQRTTN.LT.THRTT) THEN 1735! CALL DSCAL(NCCVAR,(D1/THRTT),PNEW(1),1) 1736! TT = DDOT(NCCVAR,PNEW(1),1,QNEW(1),1) 1737! SQRTT = SQRT(TT) 1738! END IF 1739! CALL DSCAL(NCCVAR,(D1/TTN),PNEW(1),1) 1740! 1741 !distribute norm between the two vectors 1742 xn=D1/sqrt((abs(ttn))) 1743 factor=sign(xn,ttn) 1744 write(lupri,*)'BiOrthPQ: factor 1/sqrt(|<p|q>|)=',factor 1745 CALL DSCAL(NCCVAR,factor,PNEW(1),1) 1746 CALL DSCAL(NCCVAR,factor,QNEW(1),1) 1747 write(lupri,*)'BiOrthPQ 2: <p|q>=', 1748 & DDOT(NCCVAR,PNEW(1),1,QNEW(1),1) 1749 write(lupri,*)'BiOrthPQ 2: <q|q>=', 1750 & DDOT(NCCVAR,QNEW(1),1,QNEW(1),1) 1751C 1752 !IF (IROUND.GT.0 .AND. ITURN.EQ.1 ) GO TO 1500 1753 IF (ITURN.EQ.1 ) GO TO 1500 1754 1755C 1756C *** End of subroutine CCBiORTPQ 1757C 1758 IF (LOCDBG) THEN 1759 WRITE (LUPRI,'(/A//)')' --- End debug output from CCBiORTPQ ---' 1760 END IF 1761 1762 RETURN 1763 END 1764*============================================= 1765C /* Deck rgord2 */ 1766 SUBROUTINE RGORD2(NM,N,WR,WI,ZR,ZL,LSORT) 1767C 1768C 15-Aug-1989 Hans Joergen Aa. Jensen 1769C modified by Ove Christiansen 20-dec 1994 1770C to normalize correct! 1771C modified by Sonia Coriani 28-oct 2010 to 1772C to sort both left and right 1773C 1774C Reorder (no normalize!) eigenpairs from DGEEV 1775C (they are already normalized in DGEEV) 1776C LSORT = .FALSE. --> ascending order 1777C LSORT = .TRUE. --> descending order 1778#include "implicit.h" 1779 1780 LOGICAL LSORT 1781 DIMENSION WR(N), WI(N), ZR(NM,N), ZL(NM,N) 1782 PARAMETER ( D0 = 0.0D0, D1 = 1.0D0 ) 1783 DO 100 I = 1,N-1 1784 ILOW = I 1785 XLOW = WR(I) 1786 DO 50 J = I+1,N 1787 IF (.NOT.LSORT) THEN 1788 !descending order 1789 IF (WR(J) .LT. XLOW) THEN 1790 ILOW = J 1791 XLOW = WR(J) 1792 END IF 1793 ELSE 1794 IF (WR(J) .GT. XLOW) THEN 1795 ILOW = J 1796 XLOW = WR(J) 1797 END IF 1798 END IF 1799 50 CONTINUE 1800 IF (ILOW .NE. I) THEN 1801 WR(ILOW) = WR(I) 1802 WR(I) = XLOW 1803 XLOW = WI(ILOW) 1804 WI(ILOW) = WI(I) 1805 WI(I) = XLOW 1806 CALL DSWAP(N,ZR(1,I),1,ZR(1,ILOW),1) 1807 CALL DSWAP(N,ZL(1,I),1,ZL(1,ILOW),1) 1808 END IF 1809 100 CONTINUE 1810C 1811C Normalize eigenvectors (should already be normalized actually) 1812C 1813! DO 200 I = 1,N 1814! ZNRM = DDOT(N,ZR(1,I),1,ZR(1,I),1) 1815! write(lupri,*)'Norm RrRr', ZNRM 1816! ZNRM = D1 / SQRT(ZNRM) 1817! IMAX = IDAMAX(N,ZR(1,I),1) 1818! IF (ZR(IMAX,I) .LT. D0) ZNRM = -ZNRM 1819! CALL DSCAL(N,ZNRM,ZR(1,I),1) 1820! ! 1821! ZNRM = DDOT(N,ZL(1,I),1,ZL(1,I),1) 1822! write(lupri,*)'Norm LrLr', ZNRM 1823! ZNRM = D1 / SQRT(ZNRM) 1824! IMAX = IDAMAX(N,ZL(1,I),1) 1825! IF (ZL(IMAX,I) .LT. D0) ZNRM = -ZNRM 1826! CALL DSCAL(N,ZNRM,ZL(1,I),1) 1827! 200 CONTINUE 1828 RETURN 1829 END 1830! 1831