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 19c*DECK CC_QR2RSD 20 SUBROUTINE CC_QR2RSD(WORK,LWORK) 21C 22C----------------------------------------------------------------------------- 23C 24C Purpose: Direct calculation of Coupled Cluster 25C quadratic response second residue calculation. 26C 27C CIS, CCS, CC2, CCSD 28C 29C Ove Christiansen April 1997. 30C 31C----------------------------------------------------------------------------- 32C 33#include "implicit.h" 34#include "priunit.h" 35#include "dummy.h" 36#include "maxorb.h" 37 PARAMETER (TOLFRQ=1.0D-08,ONE=1.0D0,XMONE=-1.0D0,THR=1.0D-08) 38C 39#include "iratdef.h" 40#include "cclr.h" 41#include "ccorb.h" 42#include "ccsdsym.h" 43#include "ccsdio.h" 44#include "ccinftap.h" 45#include "ccsdinp.h" 46#include "cclrinf.h" 47#include "ccexci.h" 48#include "cclres.h" 49#include "ccroper.h" 50#include "ccqr2r.h" 51C 52 LOGICAL LTST,LCALC,LDIP 53 DIMENSION WORK(LWORK) 54 CHARACTER MODEL*10,MODELP*10,MODEL1*10,CHSYM*2,FILEX*10,FILOLD*10 55 CHARACTER LABELA*8,LABELB*8 56C 57#include "leinf.h" 58C 59#include "mxcent.h" 60#include "maxaqn.h" 61#include "symmet.h" 62C 63C------------------------------------ 64C Header of Property calculation. 65 66C 67 WRITE (LUPRI,'(1X,A,/)') ' ' 68 WRITE (LUPRI,'(1X,A)') 69 *'*********************************************************'// 70 *'**********' 71 WRITE (LUPRI,'(1X,A)') 72 *'* '// 73 *' *' 74 WRITE (LUPRI,'(1X,A)') 75 *'*---------- OUTPUT FROM COUPLED CLUSTER QUADRATIC RESPONSE >'// 76 *'--------*' 77 IF ( XOSCST ) THEN 78 WRITE (LUPRI,'(1X,A)') 79 * '* '// 80 * ' *' 81 WRITE (LUPRI,'(1X,A)') 82 * '*----- CALCULATION OF CC EXCITED STATE TRANSITION STRENGT'// 83 * 'HS ------*' 84 ENDIF 85 WRITE (LUPRI,'(1X,A)') 86 *'* '// 87 *' *' 88 WRITE (LUPRI,'(1X,A,/)') 89 *'*********************************************************'// 90 *'**********' 91C 92 MODEL = 'CCSD ' 93 IF (CC2) THEN 94 MODEL = 'CC2' 95 ENDIF 96 IF (CCS) THEN 97 MODEL = 'CCS' 98 ENDIF 99 IF (CC3 ) THEN 100 MODEL = 'CC3' 101 WRITE(LUPRI,'(/,1x,A)') 102 * 'CC3 Oscillator strengths not implemented yet' 103 RETURN 104 ENDIF 105 IF (CC1A) THEN 106 MODEL = 'CCSDT-1a' 107 WRITE(LUPRI,'(/,1x,A)') 108 * 'CC1A Oscillator strengths not implemented yet' 109 RETURN 110 ENDIF 111 IF (CCSD) THEN 112 MODEL = 'CCSD' 113 ENDIF 114C 115 IF (CIS) THEN 116 MODELP = 'CIS' 117 ELSE 118 MODELP = MODEL 119 ENDIF 120C 121 CALL AROUND( 'Calculation of '//MODELP// ' residues ') 122C 123 IF (IPRINT.GT.10) WRITE(LUPRI,*) 'CC_QR2RSD Workspace:',LWORK 124C 125C----------------------- 126C Calculate property 127C----------------------- 128C 129 CALL FLSHFO(LUPRI) 130C 131 NAQR2 = NQR2OP*NXQR2ST 132 NBQR2 = NQR2OP*NXQR2ST 133C 134 KOSCS = 1 135 KOSCSF = KOSCS + NAQR2 136 KEND1 = KOSCSF + NBQR2 137 LEND1 = LWORK - KEND1 138 CALL DZERO(WORK(KOSCS),NAQR2) 139 CALL DZERO(WORK(KOSCSF),NBQR2) 140C 141C---------------------------------------- 142C Calculate linear response residues. 143C---------------------------------------- 144C 145 DO 1000 IQR2 = 1, NXQR2ST 146 ISTI = IQR2STI(IQR2) 147 ISTF = IQR2STF(IQR2) 148 ISYMI = ISYEXC(ISTI) 149 ISYMF = ISYEXC(ISTF) 150 ISTSYI = ISTI - ISYOFE(ISYMI) 151 ISTSYF = ISTF - ISYOFE(ISYMF) 152 EIGVI = EIGVAL(ISTI) 153 EIGVF = EIGVAL(ISTF) 154 IF (IPRINT .GT. 5) THEN 155 WRITE(LUPRI,'(/,1x,A,2(/,1X,A,I3,1X,A,I3,A,F16.8))') 156 * 'Calculating quadratic response 2. residues for: ', 157 * 'State nr. ',ISTSYI,'of symmetry ',ISYMI, 158 * ' and eigenvalue: ',EIGVI, 159 * 'State nr. ',ISTSYF,'of symmetry ',ISYMF, 160 * ' and eigenvalue: ',EIGVF 161 ENDIF 162C 163 DO 2000 IOPER = 1, NQR2OP 164 ISYMA = ISYOPR(IAQR2OP(IOPER)) 165 ISYMB = ISYOPR(IBQR2OP(IOPER)) 166 ISYMAI = MULD2H(ISYMA,ISYMI) 167 ISYMBF = MULD2H(ISYMB,ISYMF) 168C 169C------------------------------ 170C Calculate residues. 171C------------------------------ 172C 173 LABELA = LBLOPR(IAQR2OP(IOPER)) 174 LABELB = LBLOPR(IBQR2OP(IOPER)) 175 IF ((ISYMF.EQ.ISYMAI).AND.(ISYMBF.EQ.ISYMI)) THEN 176 KOSCSOF = NQR2OP*(IQR2-1) + IOPER + KOSCS - 1 177 CALL CC_QR2(LABELA,ISYMA, 178 * ISTI,ISYMI,ISTF,ISYMF, 179 * WORK(KOSCSOF),WORK(KEND1),LEND1) 180 KOSCSOF2 = NQR2OP*(IQR2-1) + IOPER + KOSCSF - 1 181 CALL CC_QR2(LABELB,ISYMB, 182 * ISTF,ISYMF,ISTI,ISYMI, 183 * WORK(KOSCSOF2),WORK(KEND1),LEND1) 184 ENDIF 185 2000 CONTINUE 186 1000 CONTINUE 187C 188C--------------------------------------------------- 189C Output quadratic response residue properties. 190C IF XOSCST put into oscillator strength tensor. 191C--------------------------------------------------- 192C 193 KOSCS2 = KEND1 194 KTRS = KOSCS2 + NEXCI*NEXCI*3*3 195 KVELST = KTRS + NEXCI*NEXCI*3*3 196 KVELST2= KVELST + NEXCI*NEXCI*3*3 197 KEND2 = KVELST2+ NEXCI*NEXCI*3*3 198 LEND2 = LWORK - KEND2 199 CALL DZERO(WORK(KOSCS2),NEXCI*NEXCI*3*3) 200 CALL DZERO(WORK(KTRS),NEXCI*NEXCI*3*3) 201 CALL DZERO(WORK(KVELST),NEXCI*NEXCI*3*3) 202 CALL DZERO(WORK(KVELST2),NEXCI*NEXCI*3*3) 203C 204 WRITE(LUPRI,'(//,1X,A6,A)') MODELP(1:6), 205 * 'Transition moments between excited states in atomic units(L):' 206 WRITE(LUPRI,'(1X,A6,A,/)') '------', 207 * '-------------------------------------------------------------' 208C 209 DO IOPER = 1, NQR2OP 210 ISYMA = ISYOPR(IAQR2OP(IOPER)) 211 LABELA = LBLOPR(IAQR2OP(IOPER)) 212 ISYMB = ISYOPR(IBQR2OP(IOPER)) 213 LABELB = LBLOPR(IBQR2OP(IOPER)) 214 DO IQR2 = 1, NXQR2ST 215 ISTI = IQR2STI(IQR2) 216 ISTF = IQR2STF(IQR2) 217 ISYMI = ISYEXC(ISTI) 218 ISYMF = ISYEXC(ISTF) 219 ISYMAI = MULD2H(ISYMA,ISYMI) 220 ISYMBF = MULD2H(ISYMB,ISYMF) 221 ISTSYI = ISTI - ISYOFE(ISYMI) 222 ISTSYF = ISTF - ISYOFE(ISYMF) 223 EIGVI = EIGVAL(ISTI) 224 EIGVF = EIGVAL(ISTF) 225 IF ((ISYMF.EQ.ISYMAI).AND.(ISYMBF.EQ.ISYMI)) THEN 226 K1 = NQR2OP*(IQR2-1) + IOPER + KOSCS - 1 227 WRITE(LUPRI,'(1X,I2,1X,I2,1X,2(I2,F10.6),2X,A3, 228 & A8,A6,1X,F15.8)') 229 * IOPER,IQR2,ISTI,EIGVI,ISTF,EIGVF, 230 * '<i|',LABELA,'|f> = ',WORK(K1) 231 ENDIF 232 END DO 233 END DO 234C 235 WRITE(LUPRI,'(//,1X,A6,A)') MODELP(1:6), 236 * 'Transition moments between excited states in atomic units(R):' 237 WRITE(LUPRI,'(1X,A6,A,/)') '------', 238 * '-------------------------------------------------------------' 239C 240 DO IOPER = 1, NQR2OP 241 ISYMA = ISYOPR(IAQR2OP(IOPER)) 242 LABELA = LBLOPR(IAQR2OP(IOPER)) 243 ISYMB = ISYOPR(IBQR2OP(IOPER)) 244 LABELB = LBLOPR(IBQR2OP(IOPER)) 245 DO IQR2 = 1, NXQR2ST 246 ISTI = IQR2STI(IQR2) 247 ISTF = IQR2STF(IQR2) 248 ISYMI = ISYEXC(ISTI) 249 ISYMF = ISYEXC(ISTF) 250 ISYMBF = MULD2H(ISYMB,ISYMF) 251 ISYMAI = MULD2H(ISYMA,ISYMI) 252 ISTSYI = ISTI - ISYOFE(ISYMI) 253 ISTSYF = ISTF - ISYOFE(ISYMF) 254 EIGVI = EIGVAL(ISTI) 255 EIGVF = EIGVAL(ISTF) 256 IF ((ISYMF.EQ.ISYMAI).AND.(ISYMBF.EQ.ISYMI)) THEN 257 K1 = NQR2OP*(IQR2-1) + IOPER + KOSCSF - 1 258 WRITE(LUPRI,'(1X,I2,1X,I2,1X,2(I2,F10.6),2X,A3,A8, 259 & A6,1X,F15.8)') 260 * IOPER,IQR2,ISTF,EIGVF,ISTI,EIGVI, 261 * '<f|',LABELB,'|i> = ',WORK(K1) 262 ENDIF 263 END DO 264 END DO 265C 266 WRITE(LUPRI,'(//,1X,A6,A)') MODELP(1:6), 267 * 'transition strength between excited states in atomic units:' 268C 269 WRITE(LUPRI,'(1X,A6,A,/)') '------', 270 * '-----------------------------------------------------------' 271C 272 DO IOPER = 1, NQR2OP 273 ISYMA = ISYOPR(IAQR2OP(IOPER)) 274 ISYMB = ISYOPR(IBQR2OP(IOPER)) 275 LABELA = LBLOPR(IAQR2OP(IOPER)) 276 LABELB = LBLOPR(IBQR2OP(IOPER)) 277 DO IQR2 = 1, NXQR2ST 278 ISTI = IQR2STI(IQR2) 279 ISTF = IQR2STF(IQR2) 280 ISYMI = ISYEXC(ISTI) 281 ISYMF = ISYEXC(ISTF) 282 ISYMBF = MULD2H(ISYMB,ISYMF) 283 ISYMAI = MULD2H(ISYMA,ISYMI) 284 ISTSYI = ISTI - ISYOFE(ISYMI) 285 ISTSYF = ISTF - ISYOFE(ISYMF) 286 EIGVI = EIGVAL(ISTI) 287 EIGVF = EIGVAL(ISTF) 288 IF ((ISYMF.EQ.ISYMAI).AND.(ISYMBF.EQ.ISYMI)) THEN 289 K1 = NQR2OP*(IQR2-1) + IOPER + KOSCS - 1 290 K2 = NQR2OP*(IQR2-1) + IOPER + KOSCSF - 1 291 RESIDUE = WORK(K1)*WORK(K2) 292 IF (RESIDUE.GE.0.0D0) THEN 293 SQRRES=SQRT(RESIDUE) 294 ELSE 295 SQRRES=-SQRT(-RESIDUE) 296 ENDIF 297 WRITE(LUPRI,'(1X,A,A8,A1,A8,A2,F9.6,A,F9.6,A,F14.8, 298 & A,F12.8,A)') 299 * 'S{',LABELA,',',LABELB,'}(',EIGVI,' -> ',EIGVF,') =', 300 * RESIDUE,' ( ',SQRRES,')' 301 IF (XOSCST) THEN 302 IADR1 = 0 303 IADR2 = 0 304 IF (LABELA(1:5).EQ.'XDIPL') IADR1 = 1 305 IF (LABELA(1:5).EQ.'YDIPL') IADR1 = 2 306 IF (LABELA(1:5).EQ.'ZDIPL') IADR1 = 3 307 IF (LABELB(1:5).EQ.'XDIPL') IADR2 = 1 308 IF (LABELB(1:5).EQ.'YDIPL') IADR2 = 2 309 IF (LABELB(1:5).EQ.'ZDIPL') IADR2 = 3 310 IF ((IADR1+IADR2).GE.2) THEN 311 IEXAD = NEXCI*(ISTF - 1) + ISTI 312 IOSCS2 = 3*3*(IEXAD-1)+3*(IADR2-1)+IADR1+KOSCS2-1 313 WORK(IOSCS2) = RESIDUE 314 ENDIF 315 ENDIF 316 IF (XVELST) THEN 317 IADR1 = 0 318 IADR2 = 0 319 IF (LABELA(1:5).EQ.'XDIPV') IADR1 = 1 320 IF (LABELA(1:5).EQ.'YDIPV') IADR1 = 2 321 IF (LABELA(1:5).EQ.'ZDIPV') IADR1 = 3 322 IF (LABELB(1:5).EQ.'XDIPV') IADR2 = 1 323 IF (LABELB(1:5).EQ.'YDIPV') IADR2 = 2 324 IF (LABELB(1:5).EQ.'ZDIPV') IADR2 = 3 325 IF ((IADR1+IADR2).GE.2) THEN 326 IEXAD = NEXCI*(ISTF - 1) + ISTI 327 IOSCS2 = 3*3*(IEXAD-1)+3*(IADR2-1)+IADR1+KVELST-1 328 WORK(IOSCS2) = RESIDUE 329 ENDIF 330 ENDIF 331 ELSE 332 SQRRES = 0.0D0 333 RESIDUE = 0.0D0 334 ENDIF 335 IF (LABELA.EQ.LABELB) THEN 336 CALL WRIPRO(SQRRES,MODELP,-2, 337 * LABELA,LABELB,LABELA,LABELB, 338 * EIGVI,EIGVF,EIGVF-EIGVI,MULD2H(ISYMAI,ISYMF), 339 * 0,0,0) 340 OSCCON = (EIGVF-EIGVI)*SQRRES*SQRRES 341 CALL WRIPRO(OSCCON,MODEL,-22, 342 & LABELA,LABELB,LABELA,LABELB, 343 & EIGVI,EIGVF,EIGVF-EIGVI,MULD2H(ISYMAI,ISYMF), 344 & ISYME,1,ISTATE) 345 ENDIF 346 END DO 347c WRITE(LUPRI,*) ' ' 348 END DO 349 350C 351C------------------------------------------- 352C Perform analysis for dipole strenghts. 353C------------------------------------------- 354C 355 CALL DCOPY(3*3*NEXCI*NEXCI,WORK(KOSCS2),1,WORK(KTRS),1) 356 CALL DCOPY(3*3*NEXCI*NEXCI,WORK(KVELST),1,WORK(KVELST2),1) 357C 358C----------------------------------------------- 359C Write out strength for CCS, CC2, and CCSD. 360C----------------------------------------------- 361C 362 LUOSC = LURES 363 IF (XOSCST.AND. (CCS.OR.CC2.OR.CCSD)) THEN 364C 365 WRITE(LUOSC,'(//A)') 366 * ' +==============================================' 367 * //'===============================+' 368 WRITE(LUOSC,'(1X,A26,A10,A)') 369 * '| Sym.I/F | I / F | ',MODELP,' excited sta' 370 * //'te transition properties |' 371 WRITE(LUOSC,'(A)') 372 * ' | | +-----------------------------' 373 * //'------------------------------+' 374 WRITE(LUOSC,'(1X,A)') 375 * '| | | Dipole Strength(a.u.) | Oscillator stre' 376 * //'ngth | Direction |' 377 WRITE(LUOSC,'(A)') 378 * ' +==============================================' 379 * //'===============================+' 380C 381 DO 9001 ISYMF = 1, NSYM 382 DO 9002 ISYMI = 1, ISYMF 383 DO 9003 IEXF = 1, NCCEXCI(ISYMF,1) 384 DO 9004 IEXI = 1, NCCEXCI(ISYMI,1) 385 386 ISTI = ISYOFE(ISYMI) + IEXI 387 ISTF = ISYOFE(ISYMF) + IEXF 388 IEXAD = NEXCI*(ISTF - 1) + ISTI 389 IF ((.NOT.SELQR2).AND.(ISTI.GE.ISTF)) GO TO 9004 390 391 LCALC = .FALSE. 392 LDIP = .TRUE. 393 DO IQR2 = 1, NXQR2ST 394 ISTII = IQR2STI(IQR2) 395 ISTIF = IQR2STF(IQR2) 396 IF ((ISTII.EQ.ISTI).AND.(ISTIF.EQ.ISTF)) LCALC = .TRUE. 397 END DO 398 KOSCSI = KOSCS2 + 3*3*(IEXAD-1) 399 KTRSI = KTRS + 3*3*(IEXAD-1) 400 CALL CC_XOSCPRI(WORK(KTRSI),WORK(KOSCSI), 401 * EIGVAL(ISTI),EIGVAL(ISTF), 402 * IEXI,ISYMI,IEXF,ISYMF, 403 * WORK(KEND2),LEND2,MODELP,LCALC, 404 * LDIP,LUOSC) 405 9004 CONTINUE 406 9003 CONTINUE 407 408 IF (.NOT.((ISYMI.EQ.NSYM).OR. 409 * (NCCEXCI(ISYMI,1).EQ.0).OR. 410 * (NCCEXCI(ISYMF,1).EQ.0))) THEN 411 NREST = 0 412 DO 9013 ISYM2 = ISYMI+1,NSYM 413 NREST = NREST + NCCEXCI(ISYM2,1) 414 9013 CONTINUE 415 IF (NREST.EQ.0) GOTO 9002 416 WRITE(LUOSC,'(A)') 417 * ' +----------------------------------------------' 418 * //'-------------------------------+' 419 ENDIF 420 9002 CONTINUE 421 9001 CONTINUE 422C 423 WRITE(LUOSC,'(A)') 424 * ' +==============================================' 425 * //'===============================+' 426C 427 ENDIF 428C 429 LUOSC = LURES 430 IF (XVELST.AND. (CCS.OR.CC2.OR.CCSD)) THEN 431C 432 WRITE(LUOSC,'(//A)') 433 * ' +==============================================' 434 * //'===============================+' 435 WRITE(LUOSC,'(1X,A26,A10,A)') 436 * '| Sym.I/F | I / F | ',MODELP,' excited sta' 437 * //'te transition properties |' 438 WRITE(LUOSC,'(A)') 439 * ' | | +-----------------------------' 440 * //'------------------------------+' 441 WRITE(LUOSC,'(1X,A)') 442 * '| | | Veloc. Strength(a.u.) | Oscillator stre' 443 * //'ngth | Direction |' 444 WRITE(LUOSC,'(A)') 445 * ' +==============================================' 446 * //'===============================+' 447C 448 DO 9006 ISYMF = 1, NSYM 449 DO 9007 ISYMI = 1, ISYMF 450 DO 9008 IEXF = 1, NCCEXCI(ISYMF,1) 451 DO 9009 IEXI = 1, NCCEXCI(ISYMI,1) 452 453 ISTI = ISYOFE(ISYMI) + IEXI 454 ISTF = ISYOFE(ISYMF) + IEXF 455 IEXAD = NEXCI*(ISTF - 1) + ISTI 456 IF ((.NOT.SELQR2).AND.(ISTI.GE.ISTF)) GO TO 9009 457 458 LCALC = .FALSE. 459 LDIP = .FALSE. 460 DO IQR2 = 1, NXQR2ST 461 ISTII = IQR2STI(IQR2) 462 ISTIF = IQR2STF(IQR2) 463 IF ((ISTII.EQ.ISTI).AND.(ISTIF.EQ.ISTF)) LCALC = .TRUE. 464 END DO 465 KOSCSI = KVELST + 3*3*(IEXAD-1) 466 KTRSI = KVELST2+ 3*3*(IEXAD-1) 467 CALL CC_XOSCPRI(WORK(KTRSI),WORK(KOSCSI), 468 * EIGVAL(ISTI),EIGVAL(ISTF), 469 * IEXI,ISYMI,IEXF,ISYMF, 470 * WORK(KEND2),LEND2,MODELP,LCALC, 471 * LDIP,LUOSC) 472 9009 CONTINUE 473 9008 CONTINUE 474 475 IF (.NOT.((ISYMI.EQ.NSYM).OR. 476 * (NCCEXCI(ISYMI,1).EQ.0).OR. 477 * (NCCEXCI(ISYMF,1).EQ.0))) THEN 478 NREST = 0 479 DO 9018 ISYM2 = ISYMI+1,NSYM 480 NREST = NREST + NCCEXCI(ISYM2,1) 481 9018 CONTINUE 482 IF (NREST.EQ.0) GOTO 9007 483 WRITE(LUOSC,'(A)') 484 * ' +----------------------------------------------' 485 * //'-------------------------------+' 486 ENDIF 487 9007 CONTINUE 488 9006 CONTINUE 489C 490 WRITE(LUOSC,'(A)') 491 * ' +==============================================' 492 * //'===============================+' 493C 494 ENDIF 495C 496 END 497c*DECK CC_QR2 498 SUBROUTINE CC_QR2(LABELA,ISYMA, 499 * ISTI,ISYMI,ISTF,ISYMF, 500 * RES1,WORK,LWORK) 501C 502C------------------------------------------------------------------------ 503C 504C Purpose: Calculate second residue of the quadratic response function. 505C 506C Written by Ove Christiansen 26-4-1996 507C 508C------------------------------------------------------------------------ 509C 510#include "implicit.h" 511#include "priunit.h" 512#include "maxorb.h" 513#include "ccorb.h" 514#include "iratdef.h" 515#include "cclr.h" 516#include "ccsdsym.h" 517#include "ccsdio.h" 518#include "ccsdinp.h" 519#include "ccexci.h" 520#include "ccqr2r.h" 521#include "dummy.h" 522C 523 PARAMETER( TWO = 2.0D00,TOLFRQ=1.0D-08 ) 524 525 DIMENSION WORK(LWORK) 526 CHARACTER*8 LABELA,MODEL*10 527C 528 IF ( IPRINT .GT. 10 ) THEN 529 CALL AROUND( 'IN CC_QR2: Calculating residues ') 530 ENDIF 531C 532C------------------------ 533C Allocate workspace. 534C------------------------ 535C 536 IF (ISYMA.NE.MULD2H(ISYMI,ISYMF)) 537 * CALL QUIT('Symmetry mismatch-2 in CC_QR2R') 538C 539 ISYMAI = MULD2H(ISYMA,ISYMI) 540 NTAMPAI = NT1AM(ISYMAI) + NT2AM(ISYMAI) 541 IF ( CCS ) NTAMPAI = NT1AM(ISYMAI) 542C 543 KETA = 1 544 KEND1 = KETA + NTAMPAI 545 LEND1 = LWORK - KEND1 546 IF (LEND1 .LT. 0) 547 * CALL QUIT('Insufficient space for allocation in CC_QR2-1') 548 KR1 = KEND1 549 KR11 = KEND1 550 KR12 = KEND1 + NT1AM(ISYMAI) 551 KEND2 = KR1 + NTAMPAI 552 LEND2 = LWORK - KEND2 553 IF (LEND2 .LT. 0) 554 * CALL QUIT('Insufficient space for allocation in CC_QR2-2') 555C 556C-------------------------------------- 557C Calculate Aa-matrix contribution. 558C-------------------------------------- 559C 560 CALL CC_ETAC(ISYMA,LABELA,WORK(KETA),'LE',ISTI,0, 561 * DUMMY,WORK(KEND1),LEND1) 562C 563 KR11 = KR1 564 KR12 = KR1 + NT1AM(ISYMF) 565 IOPT = 3 566 CALL CC_RDRSP('RE',ISTF,ISYMF,IOPT,MODEL,WORK(KR11), 567 * WORK(KR12)) 568C 569 EATBCN = DDOT(NTAMPAI,WORK(KETA),1,WORK(KR1),1) 570C 571 IF ( IPRINT .GT. 9 ) THEN 572 WRITE(LUPRI,*) ' Singles contribution:', 573 * DDOT(NT1AM(ISYMAI),WORK(KETA),1,WORK(KR1),1) 574 IF (.NOT. CCS) WRITE(LUPRI,*) ' Doubles contribution:', 575 * DDOT(NT2AM(ISYMAI),WORK(KETA+NT1AM(ISYMAI)),1, 576 * WORK(KR1+NT1AM(ISYMAI)),1) 577 ENDIF 578C 579C------------------------------------ 580C Add to response function array. 581C------------------------------------ 582C 583 IF (IPRINT .GT. 2 ) THEN 584 WRITE(LUPRI,'(1X,A3,A8,A3,A,F10.6)') 585 * '<i|',LABELA,'|f>',' LEi*A*REf cont. = ',EATBCN 586 ENDIF 587 RES1 = EATBCN + RES1 588C 589 NTAMPA = NT1AM(ISYMA) + NT2AM(ISYMA) 590 IF ( CCS ) NTAMPA = NT1AM(ISYMA) 591 KETA = 1 592 KEND1 = KETA + NTAMPA 593 LEND1 = LWORK - KEND1 594 KR1 = KEND1 595 KR11 = KEND1 596 KR12 = KEND1 + NT1AM(ISYMA) 597 KEND2 = KETA + NTAMPA 598 LEND2 = LWORK - KEND2 599 IF (LEND2 .LT. 0) 600 * CALL QUIT('Insufficient space for allocation in CC_QR2-2') 601C 602C------------------------------------- 603C Calculate B-matrix contribution. 604C------------------------------------- 605C 606 IF ((.NOT. CIS).AND.(.NOT.QR22N1)) THEN 607 WRITE(LUPRI,*) 'Have not been programmed this stupid way. ' 608 ENDIF 609 IF ((.NOT. CIS).AND.QR22N1) THEN 610 CALL CC_XKSI(WORK(KETA),LABELA,ISYMA,0,DUMMY,WORK(KEND1),LEND1) 611 ILSTNR = IN2AMP(ISTI,-EIGVAL(ISTI),ISYMI, 612 * ISTF,EIGVAL(ISTF),ISYMF) 613 CALL CC_RDRSP('N2',ILSTNR,ISYMA,IOPT,MODEL,WORK(KR11), 614 * WORK(KR12)) 615 ENDIF 616C 617 EATBCN = DDOT(NTAMPA,WORK(KETA),1,WORK(KR1),1) 618C 619 IF (IPRINT .GT. 9 .AND. (.NOT.CIS)) THEN 620 WRITE(LUPRI,*) ' Singles contribution:', 621 * DDOT(NT1AM(ISYMA),WORK(KETA),1,WORK(KR1),1) 622 IF (.NOT. CCS) WRITE(LUPRI,*) ' Doubles contribution:', 623 * DDOT(NT2AM(ISYMA),WORK(KETA+NT1AM(ISYMA)),1, 624 * WORK(KR1+NT1AM(ISYMA)),1) 625 ENDIF 626C 627C------------------------------------ 628C Add to response function array. 629C------------------------------------ 630C 631 IF (IPRINT .GT. 2 .AND.(.NOT.CIS)) THEN 632 IF (.NOT.QR22N1) THEN 633 WRITE(LUPRI,'(1X,A3,A8,A3,A,F10.6)') 634 * '<i|',LABELA,'|f>',' LEi*B*REf*t = ',EATBCN 635 ELSE 636 WRITE(LUPRI,'(1X,A3,A8,A3,A,F10.6)') 637 * '<i|',LABELA,'|f>',' Nif*KsiAcont. = ',EATBCN 638 ENDIF 639 ENDIF 640 IF (.NOT.CIS) RES1 = EATBCN + RES1 641C 642 RETURN 643 END 644 SUBROUTINE CC_XOSCPRI(TRS,OSC,EIGVI,EIGVF,IEXI,ISYMI,IEXF,ISYMF, 645 * WORK,LWORK,MODEL,LCALC,LDIP,LUOSC) 646C 647C------------------------------------------------------------------------------ 648C 649C Write out transition/oscillator strength between excited states. 650C 651C Written by Ove Christiansen April 1997 based on CC_OSCPRI 652C 653C------------------------------------------------------------------------------ 654C 655#include "implicit.h" 656#include "priunit.h" 657#include "dummy.h" 658#include "maxorb.h" 659 PARAMETER (TOLFRQ = 1.0D-08,ONE= 1.0D0,THR = 1.0D-08) 660C 661#include "iratdef.h" 662#include "cclr.h" 663#include "ccorb.h" 664#include "ccsdsym.h" 665#include "ccsdio.h" 666#include "ccsdinp.h" 667#include "ccqr2r.h" 668C 669 DIMENSION OSC(*),PVAL(3),PAXIS(3,3) 670 CHARACTER MODEL*10,CDIP*7 671 LOGICAL LCALC,LDIP 672C 673 IF ( IPRINT .GT. 10 ) THEN 674 CALL AROUND( 'IN CC_XOSCPRI: Output transition properties ' ) 675 ENDIF 676C 677C------------------------------------------ 678C write out transition strength matrix. 679C------------------------------------------ 680C 681 IF (LCALC) THEN 682 EIGVFI = EIGVF - EIGVI 683 WRITE(LUPRI,'(//,1X,A6,A,I2,A,I3,A,I2,A,I3,A,/,A,F11.8,1X, 684 * F11.8,A,F11.8)') 685 * MODEL(1:6), 686 * 'Transition strength matrix for (',ISYMI,',',IEXI,') -> (', 687 * ISYMF,',',IEXF,') transition', 688 * ' Excitation energies:',EIGVI,EIGVF,' and f-i energy ',EIGVFI 689 IF (LDIP) THEN 690 WRITE(LUPRI,'(1X,A)') 'Dipole gauge' 691 ELSE 692 WRITE(LUPRI,'(1X,A)') 'Velocity gauge' 693 ENDIF 694 CALL OUTPUT(TRS,1,3,1,3,3,3,1,LUPRI) 695C 696 CALL TNSRAN(TRS,PVAL,PAXIS, 697 * ALFSQ,BETSQ,ITST,ITST2, 698 * APAR1,APEN1,XKAPPA,IPAR) 699 WRITE(LUPRI,'(/,1X,A,/)') 700 * 'Principal values of diagonalized transition strength matrix:' 701 WRITE(LUPRI,'(1X,A)') ' a.u. ' 702 WRITE(LUPRI,'(1X,A,F12.8)') '1 ',PVAL(1) 703 WRITE(LUPRI,'(1X,A,F12.8)') '2 ',PVAL(2) 704 WRITE(LUPRI,'(1X,A,F12.8)') '3 ',PVAL(3) 705 WRITE(LUPRI,'(/,1X,A,/)') 706 * 'Principal axis of diagonalized transition strength matrix:' 707 CALL OUTPUT(PAXIS,1,3,1,3,3,3,1,LUPRI) 708 TRA = PVAL(1)+PVAL(2)+PVAL(3) 709C 710 IF (IPAR .EQ.1) CDIP = ' X ' 711 IF (IPAR .EQ.2) CDIP = ' Y ' 712 IF (IPAR .EQ.3) CDIP = ' Z ' 713 IF (IPAR .EQ.4) CDIP = ' (X,Y) ' 714 IF (IPAR .EQ.5) CDIP = ' (X,Z) ' 715 IF (IPAR .EQ.6) CDIP = ' (Y,Z) ' 716 IF (IPAR .EQ.7) CDIP = '(X,Y,Z)' 717 IF (IPAR .EQ.8) CDIP = ' - ' 718C 719C------------------------------------------ 720C First scale it - then 721C write out oscillator strength matrix. 722C------------------------------------------ 723C 724 IF (LDIP) THEN 725 FACT = EIGVFI*2.0D0/3.0D0 726 ELSE 727 FACT = 2.0D0/(3.0D0*EIGVFI) 728 ENDIF 729 CALL DSCAL(3*3,FACT,OSC,1) 730 WRITE(LUPRI,'(//,1X,A6,A,I2,A,I3,A,I2,A,I3,A,/,A,F11.8,1X, 731 * F11.8,A,F11.8/)') 732 * MODEL(1:6), 733 * 'Oscillator strength matrix for (',ISYMI,',',IEXI,') -> (', 734 * ISYMF,',',IEXF,') transition', 735 * ' Excitation energies:',EIGVI,EIGVF,' and f-i energy ',EIGVFI 736 CALL OUTPUT(OSC,1,3,1,3,3,3,1,LUPRI) 737 CALL TNSRAN(OSC,PVAL,PAXIS, 738 * ALFSQ,BETSQ,ITST,ITST2, 739 * APAR2,APEN2,XKAPPA,IPAR) 740 WRITE(LUPRI,'(/,1X,A,/)') 741 * 'Principal values of diagonalized oscillator strength matrix:' 742 WRITE(LUPRI,'(1X,A)') ' a.u. ' 743 WRITE(LUPRI,'(1X,A,F12.8)') '1 ',PVAL(1) 744 WRITE(LUPRI,'(1X,A,F12.8)') '2 ',PVAL(2) 745 WRITE(LUPRI,'(1X,A,F12.8)') '3 ',PVAL(3) 746 WRITE(LUPRI,'(/,1X,A,/)') 747 * 'Principal axis of diagonalized oscillator strength matrix:' 748 CALL OUTPUT(PAXIS,1,3,1,3,3,3,1,LUPRI) 749 OSCS = PVAL(1)+PVAL(2)+PVAL(3) 750c 751 NR = 2 752 IF ((.NOT.SELQR2).AND.(ISYMI.EQ.ISYMF)) NR = 3 753c IF ((IEXF+IEXI).EQ.NR) THEN 754 WRITE(LUOSC,9988) ISYMI,ISYMF,IEXI,IEXF,TRA,OSCS,CDIP 755c ELSE 756c WRITE(LUOSC,9989) IEXI,IEXF,TRA,OSCS,CDIP 757c ENDIF 758C 759 ELSE IF (.NOT.LCALC) THEN 760 CDIP = ' ? ' 761c IF ((IEXF+IEXI).EQ.NR) THEN 762 WRITE(LUOSC,9986) ISYMI,ISYMF,IEXI,IEXF,'Not calculated', 763 * 'Not calculated',CDIP 764c ELSE 765c WRITE(LUOSC,9987) IEXI,IEXF,'Not calculated', 766c * 'Not calculated',CDIP 767c ENDIF 768 ENDIF 769C 770 9986 FORMAT(1X,'|',I3,I3,' |',I3,I3,' | ',A16,4X, 771 * ' |',A15,5X,' | ',A7,' ',1X,' | ') 772 9987 FORMAT(1X,'| |',I3,I3,' | ',A16,4X, 773 * ' |',A15,5X,' | ',A7,' ',1X,' | ') 774 9988 FORMAT(1X,'|',I3,I3,' |',I3,I3,' | ',F16.7,4X, 775 * ' |',F15.7,5X,' | ',A7,' ',1X,' | ') 776 9989 FORMAT(1X,'| |',I3,I3,' | ',F16.7,4X, 777 * ' |',F15.7,5X,' | ',A7,' ',1X,' | ') 778C 779 END 780