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_EXCI 20 SUBROUTINE CC_EXCI(WORK,LWORK,LIST,APROXR12) 21C 22C---------------------------------------------------------------------- 23C 24C Purpose: Direct calculation of Coupled Cluster 25C excitation energies. 26C 27C Singlet: 28C 29C CCS, CC2, CCSD, CC3 30C 31C CCSDT-1a, CCSDT-1b 32C 33C CC(2), CCSDR(T), CCSDR(3), CCSDR(1A), CCSDR(1B) 34C 35C CCS = TDA = CIS 36C CC(2) = CIS(D) 37C 38C Triplet: 39C 40C CCS, CC2, CCSD 41C 42C The first section solves iteratively for CC excitation energies. 43C Both for left and right excitation energies. 44C The next section calculates perturbational corrections CC(2) and 45C CCSDR(). 46C 47C Written by Ove Christiansen August/November 1994 48C 49C Version 3, December 1996. 50C 51C 52C----------------------------------------------------------------------- 53C 54#include "implicit.h" 55#include "pgroup.h" 56#include "dummy.h" 57#include "priunit.h" 58#include "maxorb.h" 59#include "ccorb.h" 60#include "iratdef.h" 61#include "codata.h" 62#include "ccsections.h" 63#include "cclr.h" 64#include "ccsdsym.h" 65#include "ccsdio.h" 66#include "ccsdinp.h" 67#include "ccexci.h" 68#include "ccexgr.h" 69#include "ccfdgeo.h" 70#include "ccgr.h" 71#include "ccfro.h" 72#include "ccinftap.h" 73C 74#include "ccdeco.h" 75#include "ccdeco2.h" 76 77 LOGICAL LOCDBG 78 PARAMETER(LOCDBG=.FALSE.) 79C 80 LOGICAL KAJ,LINQCC,RSPIM2,LTRIP(4),LRST1,LRST2,LPROJECT,TRIPLET 81Cholesky pablo 82 LOGICAL LREDS,CCRSTRS,LREDS_SV 83 LOGICAL CCS_DONE 84 SAVE CCS_DONE 85 DATA CCS_DONE /.FALSE./ 86Cholesky pablo 87 DIMENSION WORK(LWORK),NCCEXSAV(8,3),CONT(28) 88 CHARACTER MODEL*24,MODELP*24,MODEL1*24,CHSYM*4,FILEX*10,FILOLD*10 89 CHARACTER MOPRPC*10 90 CHARACTER CDIP*1,CDUM*3, SPACE*32, MODELSCR*24 91 CHARACTER*8 FC1AM,FC2AM,FRHO1,FRHO2,FRHO12,FC12AM,FS12AM,FS2AM 92 CHARACTER*(3) LIST 93 CHARACTER*(3) APROXR12 94 PARAMETER (FC1AM ='CCR_C1AM',FC2AM ='CCR_C2AM') 95 PARAMETER (FRHO1 ='CCR_RHO1',FRHO2 ='CCR_RHO2') 96 PARAMETER (FRHO12='CCR_RH12',FC12AM='CCR_C12M',FS12AM='CCR_S12M') 97 PARAMETER (FS2AM ='CCR_S2DM') 98 PARAMETER (TWO = 2.0D00,TF=1.0D-04, TLOVLP = 0.5D0 ) 99 PARAMETER (ZERO = 0.0D0) 100 PARAMETER (XMONE = -1.0D00 , THRLDP = 1.0D-20) 101 PARAMETER (THRDEGEN = 1.0D-8) 102 CHARACTER*8 LABEL,LABEL1 103C 104#include "leinf.h" 105C 106 CALL QENTER('CC_EXCI') 107C 108C 109 ECURR = ZERO 110C 111 IF (CHOINT .AND. (.NOT. (CCS .OR. CC2))) THEN 112 WRITE(LUPRI,*) 'Error : ', 113 & 'Only CCS/CC2 singlet excitation energies have been ', 114 & 'implemented with Cholesky decomposed integrals' 115 CALL QUIT('Cholesky only for CCS/CC2 singlet excitations') 116 END IF 117c 118 LUFC1 = -1 119 LUFC2 = -1 120 LUFC12 = -1 121 LUFR1 = -1 122 LUFR2 = -1 123 LUFR12 = -1 124 LUFS12 = -1 125 LUFS2 = -1 126C 127 IF (CCP2) CCS = .TRUE. 128 IF (CHOINT .AND. CC2) CHOEXC = .TRUE. 129 TRIPLET = JACEXT 130C 131 LABEL = "Excita " 132 LABEL1 = "Exctot " 133C 134C---------------------------------------------------------- 135C Reinit NEXCI: Necessary if symmetry has been reduced. 136C---------------------------------------------------------- 137C 138 NEXCI = 0 139 NTRIP = 0 140 DO ISYM = 1,NSYM 141 ISYOFE(ISYM) = NEXCI 142 ITROFE(ISYM) = ISYOFE(ISYM) + NCCEXCI(ISYM,1) 143 NEXCI = ITROFE(ISYM) + NCCEXCI(ISYM,3) 144 NTRIP = NTRIP + NCCEXCI(ISYM,3) 145 DO IEX = ISYOFE(ISYM)+1, NEXCI 146 ISYEXC(IEX) = ISYM 147 END DO 148 DO IEX = ISYOFE(ISYM)+1, ITROFE(ISYM) 149 IMULTE(IEX) = 1 150 END DO 151 DO IEX = ITROFE(ISYM)+1, NEXCI 152 IMULTE(IEX) = 3 153 END DO 154 END DO 155 IF (IPRINT.GT.15) THEN 156 WRITE(LUPRI,*) 'IN CC_EXCI after Reinit' 157 WRITE(LUPRI,*) 'NEXCI: ',NEXCI 158 WRITE(LUPRI,*) 'Singlet: ',(NCCEXCI(J,1),J=1,NSYM) 159 WRITE(LUPRI,*) 'Triplet: ',(NCCEXCI(J,3),J=1,NSYM) 160 WRITE(LUPRI,*) 'ISYOFE:',(ISYOFE(J), J=1,NSYM) 161 WRITE(LUPRI,*) 'ITROFE:',(ISYOFE(J), J=1,NSYM) 162 WRITE(LUPRI,*) 'ISYEXC:',(ISYEXC(J), J=1,NEXCI) 163 WRITE(LUPRI,*) 'IMULTE:',(IMULTE(J), J=1,NEXCI) 164 WRITE(LUPRI,*) 'EIGVAL:',(EIGVAL(J), J=1,NEXCI) 165 ENDIF 166C 167C----------------------- 168C Initialize Leinfi. 169C----------------------- 170C 171 CALL CCLR_LEINFI(TRIPLET) 172 NTAMP = NT1AMX + NT2AMX 173 IF (CCS.OR.CCP2) NTAMP = NT1AMX 174 LINQCC = .FALSE. 175 NLOAD = 1 176C 177C------------------------------------------------------------------- 178C Test of Finite difference jacobian and linear transformations. 179C or test calculation by finite difference jacobian. 180C------------------------------------------------------------------- 181C 182 IF (LIST(1:2) .EQ. 'RE') ISIDE = +1 183 IF (LIST(1:2) .EQ. 'LE') ISIDE = -1 184 185 IF (JACTST.OR.FDJAC.OR.JACEXP) THEN 186 CALL CCLR_JAC(ISIDE,WORK,LWORK,APROXR12,TRIPLET) 187 IF (JACTST.OR.JACEXP) THEN 188 CALL AROUND( ' END OF JACTST ' ) 189 CALL QEXIT('CC_EXCI') 190 RETURN 191 ENDIF 192 ENDIF 193 IF ((.NOT.FDJAC).AND.(.NOT.JACTST).AND.FDEXCI) THEN 194 RSPIM2 = RSPIM 195 RSPIM = .FALSE. 196 FDJAC = .TRUE. 197 CALL CCLR_FDJAC(NT1AMX,NT2AMX,WORK,LWORK,APROXR12) 198 FDJAC = .FALSE. 199 RSPIM = RSPIM2 200 ENDIF 201C 202C----------------------------------------------------------------------- 203C Calculation of excit. energies from finite difference Jacobian: 204C----------------------------------------------------------------------- 205C 206 IF (FDEXCI) THEN 207 CALL AROUND('CC_EXCI: START OF CALCULATION OF FD EXCI.') 208 CALL CCLR_FDEXCI(WORK,LWORK) 209 CALL AROUND('CC_EXCI: END OF CALCULATION OF FD EXCI. ') 210 IF (.NOT.CCEXCI) THEN 211 CALL AROUND( ' ONLY CALCULATION FROM FD JACOBIAN' ) 212 CALL QEXIT('CC_EXCI') 213 RETURN 214 ENDIF 215 ENDIF 216C 217C--------------------------------------------- 218C Header of Excitation Energy calculation. 219C--------------------------------------------- 220C 221 CALL TIMER('START ',TIMEIN,TIMOUT) 222 WRITE (LUPRI,'(1X,A,/)') ' ' 223 WRITE (LUPRI,'(1X,A)') 224 *'*********************************************************'// 225 *'**********' 226 WRITE (LUPRI,'(1X,A)') 227 *'* '// 228 *' *' 229 WRITE (LUPRI,'(1X,A)') 230 *'*---------- OUTPUT FROM COUPLED CLUSTER LINEAR RESPONSE >'// 231 *'---------*' 232 IF ( CCEXCI ) THEN 233 WRITE (LUPRI,'(1X,A)') 234 * '* '// 235 * ' *' 236 WRITE (LUPRI,'(1X,A)') 237 * '*---------- CALCULATION OF EXCITATION ENERGIES >'// 238 * '---------*' 239 ENDIF 240 WRITE (LUPRI,'(1X,A)') 241 *'* '// 242 *' *' 243 WRITE (LUPRI,'(1X,A,/)') 244 *'*********************************************************'// 245 *'**********' 246C 247 MODEL = 'CCSD' 248 MOPRPC = 'CCSD ' 249 IF (CC2) THEN 250 IF (CHOINT) THEN 251 CALL AROUND( ' Cholesky-based CC2 Excitation Energies ') 252 ELSE 253 CALL AROUND( ' CC2 Excitation Energies ') 254 END IF 255 MODEL = 'CC2' 256 MOPRPC = 'CC2 ' 257 ENDIF 258 IF (CCP2) THEN 259 CALL AROUND( ' CCS and CC(2) Excitation Energies ') 260 MODEL = 'CC(2)' 261 MOPRPC ='CC(2) ' 262 ENDIF 263 IF (CCS.AND.(.NOT.CIS)) THEN 264 IF (CHOINT) THEN 265 CALL AROUND( ' Cholesky-based CCS Excitation Energies ') 266 CCS_DONE = .TRUE. 267 ELSE 268 CALL AROUND( ' CCS Excitation Energies ') 269 END IF 270 MODEL = 'CCS' 271 MOPRPC ='CCS ' 272 ENDIF 273 IF (CCS.AND.CIS) THEN 274 CALL AROUND( ' CIS Excitation Energies ') 275 MODEL = 'CCS' 276 MOPRPC ='CCS ' 277 ENDIF 278 IF (CC3 ) THEN 279 CALL AROUND( ' CC3 Excitation Energies ') 280 MODEL = 'CC3' 281 MOPRPC ='CC3 ' 282 ENDIF 283 IF (MLCC3 ) THEN 284 CALL AROUND( ' MLCC3 Excitation Energies ') 285 MODEL = 'CC3' 286 MOPRPC ='CC3 ' 287 ENDIF 288 IF (CC1B .AND. ( .NOT. CC1a) ) THEN 289 CALL AROUND( ' CCSDT-1b Excitation Energies ') 290 MODEL = 'CCSDT-1b' 291 MOPRPC ='CCSDT-1b ' 292 ENDIF 293 IF (CC1A ) THEN 294 CALL AROUND( ' CCSDT-1a Excitation Energies ') 295 MODEL = 'CCSDT-1a' 296 MOPRPC ='CCSDT-1a ' 297 ENDIF 298 IF (CCR3) THEN 299 CALL AROUND( ' CCSD and CCSDR(3) Excitation Energies') 300 MODEL = 'CCSD' 301 MOPRPC ='CCSD ' 302 ENDIF 303 IF (CCRT) THEN 304 CALL AROUND( ' CCSD and CCSDR(T) Excitation Energies') 305 MODEL = 'CCSD' 306 MOPRPC ='CCSD ' 307 ENDIF 308 IF (CCR1A) THEN 309 CALL AROUND( ' CCSD and CCSDR(1a) Excitation Energies') 310 MODEL = 'CCSD' 311 MOPRPC ='CCSD ' 312 ENDIF 313 IF (CCR1B) THEN 314 CALL AROUND( ' CCSD and CCSDR(1b) Excitation Energies') 315 MODEL = 'CCSD' 316 MOPRPC ='CCSD ' 317 ENDIF 318 IF (CCSD .AND. .NOT. MLCC3) THEN 319 CALL AROUND( ' CCSD Excitation Energies ') 320 MODEL = 'CCSD' 321 MOPRPC ='CCSD ' 322 ENDIF 323 IF (CIS) THEN 324 MODELP = 'CIS' 325 MOPRPC ='CIS ' 326 ELSE 327 MODELP = MODEL 328 IF (CCR12 .AND. (.NOT.CCS)) THEN 329 CALL CCSD_MODEL(MODELP,LENMOD,24,MODEL,24,APROXR12) 330 END IF 331 ENDIF 332C 333 IF (IPRINT.GT.10) WRITE(LUPRI,*) 'CC_EXCI-1: Workspace:',LWORK 334C 335C ******************************************** 336C ***** Determine CC excitation energies ***** 337C ******************************************** 338C 339C------------------------------ 340C Allocation of work space. 341C------------------------------ 342C 343 KEXCI = 1 344 KEXCP = KEXCI + NEXCI 345 KEXCP2 = KEXCP + NEXCI 346 KEXCP3 = KEXCP2 + NEXCI 347 KEXCP4 = KEXCP3 + NEXCI 348 KT1P = KEXCP4 + NEXCI 349 KWRK0 = KT1P + NEXCI 350 LWRK0 = LWORK - KWRK0 351 IF (LWRK0.LT. 0 ) CALL QUIT(' TOO LITTLE WORKSPACE IN CC_EXCI') 352 IF (IPRINT.GT.10) WRITE(LUPRI,*) 'CC_EXCI-1b: Workspace:',LWRK0 353C 354 KEXCPS = KEXCP 355 KEXCPS2 = KEXCP2 356 KEXCPS3 = KEXCP3 357 KEXCPS4 = KEXCP4 358 KEXCIS = KEXCI 359 KT1PS = KT1P 360C 361C------------------------------------- 362C Start loop over symmetry classes. 363C------------------------------------- 364C 365 DO 9000 ISYM = 1, NSYM 366C 367C 368C-------------------------------------- 369C Start loop over multiplicities: 370C-------------------------------------- 371C 372 DO 8000 IMULT = 1, 3, 2 373C 374C----------------------------- 375C Initialize variables. 376C----------------------------- 377C 378 IF (IMULT.EQ.1) THEN 379 TRIPLET = .FALSE. 380 ELSE IF (IMULT.EQ.3) THEN 381 TRIPLET = .TRUE. 382 ELSE 383 CALL QUIT('Illegal multiplicity in CC_EXCI.') 384 END IF 385C 386 IF (TRIPLET .AND. NCCEXCI(ISYM,IMULT).GT.0) THEN 387 IF (.NOT. (CCS.OR.CC2.OR.CCSD.OR.CC3)) THEN 388 WRITE (LUPRI,*) 'WARNING: ',MODELP, 389 & ' TRIPLET EXCITATION ENERGIES', 390 & ' ARE NOT (YET) AVAILABLE.' 391 END IF 392 IF (CHOINT) THEN 393 WRITE(LUPRI,*) 'Error : ', 394 & 'Only singlet excitation energies have been ', 395 & 'implemented with Cholesky decomposed integrals' 396 WRITE(LUPRI,*) 'Triplet excitations ignored' 397 END IF 398 END IF 399C 400 LREDS = CCR12 .OR. (CCSDT .AND. CCSDT_DIIS) 401 & .OR. (CHOINT .AND. CC2 .AND. CHEXDI) 402 LINQCC = .FALSE. 403C 404 IF (LIST(1:2) .EQ. 'RE') ISIDE = +1 405 IF (LIST(1:2) .EQ. 'LE') ISIDE = -1 406 IF ((LIST(1:2) .EQ. 'LE').AND.CCSDT) THEN 407 IF (.NOT. CC3) THEN 408 CALL QUIT('Left excitation energies and '// 409 * '(non cc3) iterative triples not allowed') 410 ENDIF 411 ENDIF 412C 413 ISYMTR = ISYM 414c IF (CCS) THEN 415 IF (CCS .OR. (CC2 .AND. CHOINT)) THEN 416 NCCVAR = NT1AM(ISYMTR) 417 ELSE 418 NCCVAR = NT1AM(ISYMTR) + NT2AM(ISYMTR) 419 IF (CCR12) NCCVAR = NCCVAR + NTR12AM(ISYMTR) 420 IF (TRIPLET) NCCVAR = NCCVAR + NT2AMA(ISYMTR) 421 END IF 422c 423 IF (TRIPLET .AND. CHOINT) THEN 424 NCCEXSAV(ISYM,IMULT) = 0 425 NCCEXCI(ISYM,IMULT) = 0 426 GOTO 8000 427 ELSE 428 CALL CCLR_LEINFI(TRIPLET) 429 END IF 430C 431C--------------- 432C OUTPUT. 433C--------------- 434C 435C 436 IF (ISYM .GT. 1) WRITE(LUPRI,'(//A)') 437 * ' ***********************************************' 438 * //'********************************' 439 WRITE (LUPRI,'(/A)') ' --------------------------' 440 WRITE (LUPRI,'(A,I5)') ' Symmetry class Nr.: ',ISYM 441 WRITE (LUPRI,'(A,I5)') ' Multiplicity : ',IMULT 442 WRITE (LUPRI,'(A,/)') ' --------------------------' 443 WRITE (LUPRI,'(A,I10)') 444 * ' Length of Excitation vectors in this class is:',NCCVAR 445C 446 IF (NCCVAR .LT.NCCEXCI(ISYM,IMULT)) THEN 447 WRITE(LUPRI,*) 'Demand for more excitation energies than ' 448 * //'amplitudes with this symmetry/multiplicity!!!!' 449 WRITE(LUPRI,*) 'Nr. of amplitudes with symmetry ',ISYM, 450 * ' and multiplicity ',IMULT,' is ',NCCVAR 451 WRITE (LUPRI,*) 'Nr. of exci.e. requested with sym ',ISYM, 452 * ' IS ',NCCEXCI(ISYM,IMULT) 453 NCCEXCI(ISYM,IMULT) = NCCVAR 454 ENDIF 455 NCCEXSAV(ISYM,IMULT) = NCCEXCI(ISYM,IMULT) 456C NCC(ISYM,IMULT) = NCCEXCI(ISYM,IMULT) 457C 458 IF (NCCEXCI(ISYMTR,IMULT).EQ.0) THEN 459 WRITE (LUPRI,'(A)') 460 * ' No Excitation Energies calculated in this symmetry class' 461 GOTO 8000 462 ENDIF 463C 464 IF (IMULT.EQ.1) THEN 465 IBOFF = ISYOFE(ISYM) 466 ELSE 467 IBOFF = ITROFE(ISYM) 468 END IF 469 NSIM = NCCEXCI(ISYM,IMULT) 470 IRST = IBOFF + 1 471 NSIMR = NSIM 472 IREND = IRST + NSIMR - 1 473C 474C--------------------------------------------------------------- 475C Allocation of work space. 476C NB: The gradient vectors is not in memory at this time. 477C--------------------------------------------------------------- 478C 479C 480 KIPLAC = KWRK0 481 KREDH = KIPLAC + MAXRED 482 KREDGD = KREDH + MAXRED*MAXRED 483 KEIVAL = KREDGD + MAXRED*NSIMR 484 KSOLEQ = KEIVAL + MAXRED 485 KWRK1 = KSOLEQ + MAXRED*MAXRED 486 487 KREDS = KWRK1 488 IF (LREDS) KWRK1 = KREDS + MAXRED*MAXRED 489 490casm IF (CCSDT .AND. CCSDT_DIIS) THEN 491 IF ((CCSDT .AND. CCSDT_DIIS) .OR. 492 & (CHOINT .AND. CC2 .AND. CHEXDI)) THEN 493 KAMAT = KWRK1 494 KITRAN = KAMAT + (MXDIIS+1)*(MXDIIS+1)*NSIM 495 KCONV = KITRAN + NSIM 496 KRNORM = KCONV + NSIM 497 KWRK1 = KRNORM + NSIM 498 END IF 499 500 LWRK1 = LWORK - KWRK1 501C 502 IF (LWRK1.LT. 0 ) 503 * CALL QUIT(' TOO LITTLE WORKSPACE IN CC_EXCI') 504 CALL DZERO(WORK(KEIVAL),MAXRED) 505C 506C-------------------------------------------------------------- 507C If CCR3 do spaghetti goto to pert. corrections section. 508C-------------------------------------------------------------- 509C 510 IF ( CCR3 ) GOTO 1200 511C 512 WRITE (LUPRI,'(A,I3,A)') 513 * ' Converging for',NCCEXCI(ISYM,IMULT),' roots.' 514C 515C----------------------------------------------------------------------- 516C Choose start omega for CC3: EOMINP contains input choice, 517C CCSD or CCSDR excitations energies. 518C Choose the first, for omesc and mxtomn this means the highest, 519C and thus probably the hightest shift. 520C----------------------------------------------------------------------- 521C 522 IF (CCSDT .OR. MLCC3) THEN 523 ECURR = EOMINP(1,ISYMTR,IMULT) 524 ENDIF 525C 526 IF ( STSD.AND.CCSDT ) THEN 527 CCSDT = .FALSE. 528 KAJ = .TRUE. 529 ENDIF 530C 531C------------------- 532C Open files. 533C------------------- 534C 535 CALL CC_FILOP(FRHO1,LUFR1,FRHO2,LUFR2,FRHO12,LUFR12, 536 * FC1AM,LUFC1,FC2AM,LUFC2,FC12AM,LUFC12, 537 * FS12AM,LUFS12,FS2AM,LUFS2) 538C 539C----------------------------- 540C Create start vectors. 541C----------------------------- 542C 543 LPROJECT = .FALSE. 544 CALL CCEQ_STR(FC1AM,LUFC1,FC2AM,LUFC2,FC12AM,LUFC12, 545 * FS12AM,LUFS12,FS2AM,LUFS2,LPROJECT,ISTATPRJ, 546 * TRIPLET,ISIDE,IRST,NSIMR,NUPVEC, 547 * NREDH,WORK(KEIVAL),WORK(KIPLAC), 548 * WORK(KWRK1),LWRK1,LIST) 549 CALL FLSHFO(LUPRI) 550c 551C 552C------------------------------------------ 553C Solve equations by call to solver. 554C------------------------------------------ 555C 556casm IF ((OMESC .AND. CCSDT .AND. CCSDT_DIIS)) THEN 557 IF ((OMESC .AND. CCSDT .AND. CCSDT_DIIS) .OR. 558 & (CHOINT .AND. CC2 .AND. CHEXDI)) THEN 559 560 IF (CHOINT) THEN 561 DO ICOUNT = 1,NSIM 562 IF (CCS_DONE) THEN 563 WORK(KEIVAL-1+ICOUNT) = EIGVAL(ISYOFE(ISYM)+ICOUNT) 564 ELSE 565 WORK(KEIVAL-1+ICOUNT) = 0.0D0 566 END IF 567 END DO 568C 569 IF (DV4DIS) THEN 570 NUPVEC_SV = NUPVEC 571 LREDS_SV = LREDS 572 LREDS = .FALSE. 573 THRLE_SV = THRLE 574 THRLE = 1.0D-3 575 ECURR = 0.0D0 576 CALL CCEQ_SOL(LIST,LPROJECT,ISTATPRJ,ECURR, 577 * FRHO1,LUFR1,FRHO2,LUFR2,FRHO12,LUFR12, 578 * FC1AM,LUFC1,FC2AM,LUFC2,FC12AM,LUFC12, 579 * LREDS,WORK(KREDS),FS12AM,LUFS12,FS2AM,LUFS2, 580 * LINQCC,TRIPLET,ISIDE,IRST,NSIMR,NUPVEC, 581 * NREDH,WORK(KREDH),WORK(KEIVAL), 582 * WORK(KSOLEQ),WORK(KWRK1),LWRK1,APROXR12) 583C 584C Discard all trial vectors unless last one 585C 586 CCRSTRS = CCRSTR 587 CCRSTR = .TRUE. 588C 589 KIPLAC = KWRK1 590 KWRK2 = KIPLAC + MAXRED 591 LWRK2 = LWORK - KWRK2 592 IF (LWRK2 .LT. 0) 593 & CALL QUIT('Not enough space for DV4DIS') 594C 595 CALL CCEQ_STR(FC1AM,LUFC1,FC2AM,LUFC2,FC12AM,LUFC12, 596 * FS12AM,LUFS12,FS2AM,LUFS2,LPROJECT,ISTATPRJ, 597 * TRIPLET,ISIDE,IRST,NSIMR,NUPVEC, 598 * NREDH,WORK(KEIVAL),WORK(KIPLAC), 599 * WORK(KWRK2),LWRK2,LIST) 600c NUPVEC = NUPVEC_SV 601c NREDH = NUPVEC 602C 603 CCRSTR = CCRSTRS 604 LREDS = LREDS_SV 605 THRLE = THRLE_SV 606 END IF 607C 608 ELSE 609 DO ICOUNT = 1, NSIM 610 IF (TRIPLET) THEN 611 WORK(KEIVAL-1+ICOUNT) = EIGVAL(ITROFE(ISYM)+ICOUNT) 612 ELSE 613 WORK(KEIVAL-1+ICOUNT) = EIGVAL(ISYOFE(ISYM)+ICOUNT) 614 END IF 615 END DO 616 END IF 617C 618 CALL CCDIIS_SOL(FRHO1,LUFR1,FRHO2,LUFR2, 619 * FC1AM,LUFC1,FC2AM,LUFC2,LIST,IRST, 620 * TRIPLET,ISIDE,NSIMR,NUPVEC, 621 * NREDH,WORK(KREDH),WORK(KREDS),WORK(KEIVAL), 622 * WORK(KSOLEQ),WORK(KAMAT), 623 * WORK(KITRAN),WORK(KCONV),WORK(KRNORM), 624 * WORK(KWRK1),LWRK1) 625C 626 ELSE 627 CALL CCEQ_SOL(LIST,LPROJECT,ISTATPRJ,ECURR, 628 * FRHO1,LUFR1,FRHO2,LUFR2,FRHO12,LUFR12, 629 * FC1AM,LUFC1,FC2AM,LUFC2,FC12AM,LUFC12, 630 * LREDS,WORK(KREDS),FS12AM,LUFS12,FS2AM,LUFS2, 631 * LINQCC,TRIPLET,ISIDE,IRST,NSIMR,NUPVEC, 632 * NREDH,WORK(KREDH),WORK(KEIVAL), 633 * WORK(KSOLEQ),WORK(KWRK1),LWRK1,APROXR12) 634 END IF 635 636 ! test eigenvector: in case of degeneracies orthogonalize 637 ! eigenvectors to the same eigenvalue 638 CALL CCTSTSOL(WORK(KSOLEQ),WORK(KEIVAL),THRDEGEN,NREDH,NSIMR) 639 640 IF ( STSD.AND.KAJ ) THEN 641 CCSDT = .TRUE. 642 ENDIF 643C 644C------------------------------------- 645C Write out Excitation energies. 646C------------------------------------- 647C 648 WRITE (LUPRI,'(//A,I3)') ' SYMMETRY CLASS NR. ',ISYM 649 WRITE (LUPRI,'(A,I3)') ' MULTIPLICITY ',IMULT 650 IF (.NOT. CCSDT) THEN 651 IF (LIST(1:2) .EQ. 'RE') THEN 652 WRITE (LUPRI,'(//1X,A10,1X,A26)') 653 * MODELP,'right excitation energies:' 654 ELSE IF (LIST(1:2) .EQ. 'LE') THEN 655 WRITE (LUPRI,'(//1X,A10,1X,A25)') 656 * MODELP,'left excitation energies:' 657 ENDIF 658 WRITE (LUPRI,'(A)') 659 * ' ====================================' 660 ELSE 661 IF ( STSD ) THEN 662 WRITE (LUPRI,'(//A11,A10,A12)') 663 * ' CCSD with ',MODELP,' amplitudes:' 664 WRITE (LUPRI,'(A)') 665 * ' ===============================' 666 ELSE 667 WRITE (LUPRI,'(//,1X,A10,A33,F10.6)') 668 * MODELP,' excitation energies with Omega= ',ECURR 669 WRITE (LUPRI,'(A)') 670 * ' =====================================================' 671 ENDIF 672C 673 ENDIF 674C 675 WRITE (LUPRI,'(A//A/A)') 676 * ' (conversion factor used: 1 au = 27.2113957 eV)', 677 * ' Excitation no. Hartree eV ', 678 * ' -------------- ------- --' 679 DO 400 I = 1,NCCEXCI(ISYM,IMULT) 680 WRITE (LUPRI,'(I10,4X,2F20.10)') 681 * I,WORK(KEIVAL-1+I),WORK(KEIVAL-1+I)*XTEV 682 400 CONTINUE 683C 684 WRITE (LUPRI,'(//A,2I5,/A/A)') 685 * ' Total excited state energies for states of symmetry/spin', 686 * ISYM,IMULT, 687 * ' Excitation no. Energy (Hartree) ', 688 * ' ------------------------------------- ' 689 DO 401 I = 1,NCCEXCI(ISYM,IMULT) 690 WRITE (LUPRI,'(A2,2I5,F30.15)') 691 * '@@',ISYM,I,WORK(KEIVAL-1+I)+ECCGRS 692 401 CONTINUE 693C 694C------------------------------------ 695C Analysis of solution vectors. 696C------------------------------------ 697C 698 NVARPT = NCCVAR+ 2*NALLAI(ISYM) 699 KWRK2 = KWRK1 + NVARPT 700 LWRK2 = LWORK - KWRK2 701 NSIMUL = MIN(NSIMR,LWRK2/NCCVAR) 702 NBATCH = (NSIMR -1)/NSIMUL + 1 703 IOFF1 = 1 704 ICOUNT = 0 705 DO 500 I = 1,NBATCH 706 IOFF2 = MIN(NSIMUL,NCCEXCI(ISYM,IMULT) - (I-1)*NSIMUL) 707 CALL CCCONV(LUFC1,FC1AM,LUFC2,FC2AM,LUFC12,FC12AM, 708 * TRIPLET,CCR12,NREDH,IOFF1,IOFF2,WORK(KSOLEQ), 709 * WORK(KWRK2),WORK(KWRK1)) 710 711 IF (LOCDBG) THEN 712 WRITE(LUPRI,*) 'Overlap matrix for solution vectors:' 713 CALL PROVLP(WORK(KWRK2),WORK(KWRK2),NCCVAR,IOFF2, 714 * WORK(KWRK2+IOFF2*NCCVAR),LUPRI) 715 END IF 716 IF ( IPRINT .GT. 30 ) THEN 717 CALL AROUND('CC_EXCI: VECTORS IN MO BASIS' ) 718 CALL OUTPUT(WORK(KWRK2),1,NCCVAR,1,NCCEXCI(ISYM,IMULT), 719 * NCCVAR,NCCEXCI(ISYM,IMULT),1,LUPRI) 720 ENDIF 721C 722 DO 510 J = 1,IOFF2 723 ICOUNT = ICOUNT + 1 724 IF (TRIPLET) THEN 725 ILSTNR = ITROFE(ISYM) + ICOUNT 726 ELSE 727 ILSTNR = ISYOFE(ISYM) + ICOUNT 728 END IF 729 730 WRITE(LUPRI,'(//1X,2A,I3)') 'Analysis of the Coupled ', 731 * 'Cluster Excitation Vector Number :', ICOUNT 732 WRITE(LUPRI,'(1X,2A)') '-----------------------------', 733 * '--------------------------------' 734 WRITE(LUPRI,'(/10X,A,23X,F10.4,A)') 735 * 'Excitation Energy : ',WORK(KEIVAL+ICOUNT-1)*XTEV,' eV' 736 WORK(KEXCIS+(ICOUNT-1)) = WORK(KEIVAL+ICOUNT-1) 737 IF (TRIPLET) THEN 738 KCAM1 = KWRK2 + (J-1)*NCCVAR 739 KCAMP = KCAM1 + NT1AM(ISYMTR) 740 KCAMM = KCAMP + NT2AM(ISYMTR) 741 CALL CC_PRAM3(WORK(KCAM1), WORK(KCAMP), WORK(KCAMM), 742 * WORK(KT1PS+(ICOUNT-1)),PTP,PTM,ISYMTR, 743 * .FALSE.) 744 ELSE 745 CALL CC_PRAM(WORK(KWRK2+(J-1)*NCCVAR), 746 * WORK(KT1PS+(ICOUNT-1)),ISYMTR, 747 * .FALSE.) 748 END IF 749 CALL FLSHFO(LUPRI) 750 751C ----------------------------------------------------- 752C Save solution vectors on file and excitation energies 753C in the result array EIGVAL: 754C (if we compute both left and right eigenvector, 755C don't overwrite with left eigenvalues the right ones) 756C ----------------------------------------------------- 757 IOPT = 3 758 CALL CC_WRRSP(LIST,ILSTNR,ISYM,IOPT,MODELP,DUMMY, 759 * WORK(KWRK2+(J-1)*NCCVAR), 760 * WORK(KWRK2+(J-1)*NCCVAR+NT1AM(ISYM)), 761 * WORK(KWRK1),NVARPT) 762 IF (CCR12) THEN 763 IOPT = 32 764 CALL CC_WRRSP(LIST,ILSTNR,ISYM,IOPT,MODELP,DUMMY,DUMMY, 765 * WORK(KWRK2+(J-1)*NCCVAR+NT1AM(ISYM)+NT2AM(ISYM)), 766 * WORK(KWRK1),NVARPT) 767 END IF 768 769 IF (LIST(1:2).EQ.'RE'.OR.LHTR) THEN 770 EIGVAL(ILSTNR) = WORK(KEIVAL+ICOUNT-1) 771 ENDIF 772 773 510 CONTINUE 774 IOFF1 = IOFF1 + NSIMUL 775 500 CONTINUE 776C 777 IF (EXCI_CONT .AND. LIST(1:2).EQ.'RE') THEN 778 DO I = 1, NSIMR 779 KC1AM = KWRK0 780 KWRK1 = KC1AM + NT1AM(ISYM) 781 IF (.NOT.CCS) THEN 782 KC2AM = KWRK1 783 KWRK1 = KC2AM + NT2AM(ISYM) 784 END IF 785 IF (CCR12) THEN 786 KC12AM = KWRK1 787 KWRK1 = KC12AM + NTR12AM(ISYM) 788 END IF 789 LWRK1 = LWORK - KWRK1 790 IF(LWRK1.LT.0) CALL QUIT('OUT OF WORKSPACE IN CC_EXCI') 791 792 IF (TRIPLET) THEN 793 ILSTNR = ITROFE(ISYM) + I 794 ELSE 795 ILSTNR = ISYOFE(ISYM) + I 796 END IF 797 798 IOPT = 3 799 CALL CC_RDRSP(LIST,ILSTNR,ISYM,IOPT,MODEL, 800 * WORK(KC1AM),WORK(KC2AM)) 801 IF (CCR12) THEN 802 IOPT = 32 803 CALL CC_RDRSP(LIST,ILSTNR,ISYM,IOPT,MODEL, 804 * DUMMY,WORK(KC12AM)) 805 END IF 806 807 CALL CC_ATRR(ECURR,ISYM,1,WORK(KWRK0),LWRK0,.TRUE.,CONT, 808 & APROXR12,.FALSE.) 809 CALL CC_EXCIT_CONT(I,ISYM,TRIPLET,MODELP, 810 * CONT(4),CONT(1)) 811 END DO 812 END IF 813C 814 CALL FLSHFO(LUPRI) 815C 816C---------------------------------------------------------- 817C========================================================== 818C IF CCSDT we may solve until self consistency. 819C And we have to if Cholesky CC2 820C========================================================== 821C---------------------------------------------------------- 822C 823C 824 IF (OMESC .AND. ((CCSDT.AND.(.NOT.CCSDT_DIIS)) 825 & .OR. (CHOINT .AND. CC2 .AND. (.NOT. CHEXDI))) 826 & .OR. MLCC3) THEN 827C 828 CALL DZERO(WORK(KEXCIS),NCCEXCI(ISYMTR,IMULT)) 829 CALL DZERO(WORK(KT1PS),NCCEXCI(ISYMTR,IMULT)) 830C 831C----------------------------------------------------------- 832C Solve for all NOMINP eigenvalues self-consistently. 833C----------------------------------------------------------- 834C 835 DO 1000 IOM = 1, NOMINP(ISYMTR,IMULT) 836 837C 838C----------------------- 839C Solve equtions. 840C----------------------- 841C 842 ISTATE = IOMINP(IOM,ISYMTR,IMULT) 843C 844 NCCEXCI(ISYM,IMULT) = ISTATE 845C NCC(ISYM,IMULT) = ISTATE 846 NSIMR = ISTATE 847 NSIDE = 1 848C 849C-------------------------------- 850C Create start vectors. 851C-------------------------------- 852C 853 LPROJECT = .FALSE. 854Cholesky pablo 855 IF (CHOINT .AND. CC2) THEN 856 CCRSTRS = CCRSTR 857 CCRSTR = .TRUE. 858 CALL CCEQ_STR(FC1AM,LUFC1,FC2AM,LUFC2,FC12AM,LUFC12, 859 * FS12AM,LUFS12,FS2AM,LUFS2,LPROJECT,ISTATPRJ, 860 * TRIPLET,ISIDE,IRST,NSIMR,NUPVEC, 861 * NREDH,WORK(KEIVAL),WORK(KIPLAC), 862 * WORK(KWRK1),LWRK1,LIST) 863 CCRSTR = CCRSTRS 864 ELSE 865 CALL CCEQ_STR(FC1AM,LUFC1,FC2AM,LUFC2,FC12AM,LUFC12, 866 * FS12AM,LUFS12,FS2AM,LUFS2,LPROJECT,ISTATPRJ, 867 * TRIPLET,ISIDE,IRST,NSIMR,NUPVEC, 868 * NREDH,WORK(KEIVAL),WORK(KIPLAC), 869 * WORK(KWRK1),LWRK1,LIST) 870 END IF 871Cholesky pablo 872 CALL FLSHFO(LUPRI) 873C 874C------------------------------------------------------------------- 875C If EOMINP read in is different from zero, use this 876C omega. 877C------------------------------------------------------------------- 878C 879 ITSC = 1 880 IF ( STSD ) ITSC = 0 881C 882 IF (TRIPLET) THEN 883 ILSTNR = ITROFE(ISYM) + ISTATE 884 ELSE 885 ILSTNR = ISYOFE(ISYM) + ISTATE 886 ENDIF 887C ECURR1 = EIGVAL(ISTATE) 888 ECURR1 = EIGVAL(ILSTNR) 889 IF (ITSC .EQ. 0 ) THEN 890 EITOL = 1.0D-3 891 IF (ABS(EOMINP(IOM,ISYMTR,IMULT)).GT.EITOL) THEN 892 ECURR1 = EOMINP(IOM,ISYMTR,IMULT) 893 WRITE(LUPRI,'(/,1X,A24,F10.6)') 894 * 'Omega put equal to input',ECURR1 895 ENDIF 896 ENDIF 897C 898C------------------------------------------------- 899C 900C LOOP OMESC 901C 902C Loop until solution is selfconsistent. 903C 904C------------------------------------------------- 905C 906 2000 CONTINUE 907C 908 ITSC = ITSC + 1 909 ECURR = ECURR1 910C 911 IF (CHOINT) THEN 912 WRITE(LUPRI,'(A,I3,A,F8.5)') 'Calculating state :', 913 & ISTATE, ' with present energy :',ecurr 914 CALL FLSHFO(LUPRI) 915 END IF 916C 917C--------------------------------------------- 918C Solve equations by call to solver. 919C--------------------------------------------- 920C 921 CALL CCEQ_SOL(LIST,LPROJECT,ISTATPRJ,ECURR, 922 * FRHO1,LUFR1,FRHO2,LUFR2,FRHO12,LUFR12, 923 * FC1AM,LUFC1,FC2AM,LUFC2,FC12AM,LUFC12, 924 * LREDS,WORK(KREDS),FS12AM,LUFS12,FS2AM,LUFS2, 925 * LINQCC,TRIPLET,ISIDE,IRST,NSIMR,NUPVEC, 926 * NREDH,WORK(KREDH),WORK(KEIVAL), 927 * WORK(KSOLEQ),WORK(KWRK1),LWRK1,APROXR12) 928C 929C------------------------------------- 930C Write out Excitation energies. 931C------------------------------------- 932C 933 WRITE (LUPRI,'(//,1X,A10,A33,F10.6)') 934 * MODELP,' excitation energies with Omega= ',ECURR 935 WRITE (LUPRI,'(A/A//A/A/)') 936 * ' =====================================================', 937 * ' (conversion factor used: 1 au = 27.2113957 eV)', 938 * ' Excitation no. Hartree eV', 939 * ' -------------- ------- --' 940C 941 DO 900 I = 1,NCCEXCI(ISYM,IMULT) 942C 943 WRITE (LUPRI,'(I10,2F20.6)') 944 * I,WORK(KEIVAL-1+I),WORK(KEIVAL-1+I)*XTEV 945C 946 900 CONTINUE 947C 948 ECURR1 = WORK(KEIVAL-1+ISTATE) 949C 950 IF (IPRINT.GT. 1) THEN 951 WRITE(LUPRI,'(/1x,A13,I2,A15,2X,F10.6)') 952 * 'Omega in the ',ITSC,'-th iteration: ',ECURR1 953 WRITE(LUPRI,'(1x,A32,F10.6)') 954 * 'Omega in previous iteration: ',ECURR 955 WRITE(LUPRI,'(1x,A32,F10.6,/)') 956 * 'Tolerance for self consistency: ',TOLSC 957 ENDIF 958C 959 IF (ABS(ECURR1-ECURR).LT. TOLSC) THEN 960C 961 WRITE(LUPRI,'(/,1x,A23,F10.6,A11,I3,A4,I3,A11/)') 962 * 'Converged root to diff.' 963 * ,ECURR-ECURR1,' for root ',ISTATE 964 * ,' in ', ITSC ,' iterations' 965 EOMINP(IOM,ISYMTR,IMULT) = ECURR1 966C 967 WRITE (LUPRI,'(//A,2I3,/A/A)') 968 * ' Total excited state energies for states of symmetry/spin ', 969 * ISYM,IMULT, 970 * ' Excitation no. Energy (Hartree) ', 971 * ' ------------------------------------- ' 972 WRITE (LUPRI,'(A3,2I5,F30.15)') 973C 974CKH * '@@@',ISYM,I,WORK(KEIVAL-1+I)+ECCGRS 975C 976 * '@@@',ISYM,I-1,EOMINP(IOM,ISYMTR,IMULT)+ECCGRS 977C 978C 979 ENDIF 980C 981C---------------------------------------------------------------------- 982C If equations not are selfconsistent with respect to omega, 983C then take solution vectors and energy and start again. 984C Else write out the analysis of the vectors. 985C---------------------------------------------------------------------- 986C 987 NVARPT = NCCVAR+ 2*NALLAI(ISYM) 988 KWRK2 = KWRK1 + NVARPT 989 LWRK2 = LWORK - KWRK2 990 991 THRESH = 0.05 992 MAXLIN = 100 993 NSIMUL = MIN(NCCEXCI(ISYM,IMULT),LWRK2/NCCVAR -1 ) 994 NBATCH = (NCCEXCI(ISYM,IMULT)-1)/NSIMUL + 1 995 IOFF1 = 1 996 ICOUNT = 0 997C 998 NSIMA = 1 999 NSIMB = 0 1000 NBLE3 = 0 1001C 1002 DO 600 I = 1,NBATCH 1003 IOFF2 = MIN(NSIMUL,NCCEXCI(ISYM,IMULT) - (I-1)*NSIMUL) 1004 CALL CCCONV(LUFC1,FC1AM,LUFC2,FC2AM,LUFC12,FC12AM, 1005 * TRIPLET,CCR12,NREDH,IOFF1,IOFF2,WORK(KSOLEQ), 1006 * WORK(KWRK2),WORK(KWRK1)) 1007C 1008 IF ( IPRINT .GT. 30 ) THEN 1009 CALL AROUND('CC_EXCI: VECTORS IN MO BASIS' ) 1010 CALL OUTPUT(WORK(KWRK2),1,NCCVAR,1, 1011 * NCCEXCI(ISYM,IMULT), 1012 * NCCVAR,NCCEXCI(ISYM,IMULT),1,LUPRI) 1013 ENDIF 1014C 1015C----------------------------------------------------------------- 1016C If converged write out the analysis of the vectors. 1017C----------------------------------------------------------------- 1018C 1019 DO 610 J = 1,IOFF2 1020 ICOUNT = ICOUNT + 1 1021C ------------------------------- 1022C save vectors on standard files: 1023C ------------------------------- 1024 IF (TRIPLET) THEN 1025 ILSTNR = ITROFE(ISYM) + ICOUNT 1026 ELSE 1027 ILSTNR = ISYOFE(ISYM) + ICOUNT 1028 END IF 1029 IOPT = 3 1030 CALL CC_WRRSP(LIST,ILSTNR,ISYM,IOPT,MODEL,DUMMY, 1031 * WORK(KWRK2+(ICOUNT-1)*NCCVAR), 1032 * WORK(KWRK2+(ICOUNT-1)*NCCVAR+NT1AM(ISYM)), 1033 * WORK(KWRK1),NVARPT) 1034C 1035 IF (ABS(ECURR1-ECURR).LE. TOLSC) THEN 1036C 1037 WRITE(LUPRI,'(//1X,A,I2)') 1038 *'Analysis of the Coupled Cluster Excitation Vector Number : ', 1039 * ICOUNT 1040 WRITE(LUPRI,'(1X,A)') 1041 *'-------------------------------------------------------------' 1042 WRITE(LUPRI,'(/10X,A,23X,F10.4,A)') 1043 * 'Excitation Energy : ',WORK(KEIVAL+ICOUNT-1)*XTEV, 1044 * ' eV' 1045C 1046 IF (TRIPLET) THEN 1047 KCAM1 = KWRK2 + (ICOUNT-1)*NCCVAR 1048 KCAMP = KCAM1 + NT1AM(ISYMTR) 1049 KCAMM = KCAMP + NT2AM(ISYMTR) 1050 CALL CC_PRAM3(WORK(KCAM1), WORK(KCAMP), 1051 * WORK(KCAMM),PT1, 1052 * PTP,PTM,ISYMTR,.FALSE.) 1053 ELSE 1054 CALL CC_PRAM(WORK(KWRK2 + (ICOUNT-1)*NCCVAR), 1055 * PT1,ISYMTR,.FALSE.) 1056 ENDIF 1057 IF (ICOUNT .EQ. ISTATE) THEN 1058 WORK(KEXCIS+(ICOUNT-1)) = WORK(KEIVAL+ICOUNT-1) 1059 WORK(KT1PS+(ICOUNT-1)) = PT1 1060 ENDIF 1061 1062C ------------------------------- 1063C save / check excitation energy: 1064C ------------------------------- 1065 IF ((LIST(1:2) .EQ. 'LE').AND.(.NOT.LHTR)) THEN 1066 IF(ABS(EIGVAL(ILSTNR)-WORK(KEIVAL+ICOUNT-1)).GT.TF) 1067 & WRITE(LUPRI,'(////,1X,A,//)') 1068 & 'Warning - Large difference between '// 1069 & 'left and right excitation energies' 1070 ELSE 1071 EIGVAL(ILSTNR) = WORK(KEIVAL+ICOUNT-1) 1072 END IF 1073C 1074Cholesky pablo: Call to CCEQ_STR moved after the loop 1075C in order to call it just once. 1076C 1077 ENDIF 1078C 1079 610 CONTINUE 1080C 1081 IOFF1 = IOFF1 + NSIMUL 1082 600 CONTINUE 1083C 1084C--------------------------------------------------------- 1085C If not self-consistent then return to 2000 and 1086C solved eigenvalue equations with new omega. 1087C--------------------------------------------------------- 1088C 1089 IF (ABS(ECURR1-ECURR).GT. TOLSC) THEN 1090C 1091Cholesky pablo 1092 LPROJECT = .FALSE. 1093 IF (CHOINT .AND. CC2) THEN 1094 CCRSTRS = CCRSTR 1095 CCRSTR = .TRUE. 1096 CALL CCEQ_STR(FC1AM,LUFC1,FC2AM,LUFC2,FC12AM, 1097 * LUFC12,FS12AM,LUFS12,FS2AM,LUFS2, 1098 * LPROJECT,ISTATPRJ, 1099 * TRIPLET,ISIDE,IRST,NSIMR,NUPVEC, 1100 * NREDH,WORK(KEIVAL),WORK(KIPLAC), 1101 * WORK(KWRK1),LWRK1,LIST) 1102 CCRSTR = CCRSTRS 1103 ELSE 1104 CALL CCEQ_STR(FC1AM,LUFC1,FC2AM,LUFC2,FC12AM, 1105 * LUFC12,FS12AM,LUFS12,FS2AM,LUFS2, 1106 * LPROJECT,ISTATPRJ, 1107 * TRIPLET,ISIDE,IRST,NSIMR,NUPVEC, 1108 * NREDH,WORK(KEIVAL),WORK(KIPLAC), 1109 * WORK(KWRK1),LWRK1,LIST) 1110 END IF 1111Cholesky pablo 1112 CALL FLSHFO(LUPRI) 1113 GOTO 2000 1114C 1115 ENDIF 1116C 1117 1000 CONTINUE 1118C 1119 ENDIF 1120C 1121C------------------------------- 1122C Close and delete files. 1123C------------------------------- 1124C 1125 CALL CC_FILCL(FRHO1,LUFR1,FRHO2,LUFR2,FRHO12,LUFR12, 1126 * FC1AM,LUFC1,FC2AM,LUFC2,FC12AM,LUFC12, 1127 * FS12AM,LUFS12,FS2AM,LUFS2) 1128 1129C----------------------------------------------------------- 1130C restore number of excitation within symmetry class: 1131C----------------------------------------------------------- 1132 1133 NCCEXCI(ISYM,IMULT) = NCCEXSAV(ISYM,IMULT) 1134 1135C--------------------------------------------------------------------- 1136C First, check if left and right eigenvalues agree. Then 1137C normalize Eigenvectors such that (RE|RE) = 1 and (LE|RE) = 1 1138C--------------------------------------------------------------------- 1139 IF (LIST(1:2).EQ.'LE' .AND. (.NOT.LHTR)) THEN 1140 ! set threshold within which left and right eigenvalues 1141 ! should agree, else complain 1142 THREIG = 10 * THREXC 1143 IF (OMESC.AND.CCSDT.AND.(.NOT.CCSDT_DIIS)) THREIG = TOLSC 1144 1145 DO ICOUNT = 1, NCCEXCI(ISYM,IMULT) 1146 IF (TRIPLET) THEN 1147 ILSTNR = ITROFE(ISYM) + ICOUNT 1148 ELSE 1149 ILSTNR = ISYOFE(ISYM) + ICOUNT 1150 END IF 1151 IF(ABS(EIGVAL(ILSTNR)-WORK(KEXCIS+ICOUNT-1)).GT.THREIG)THEN 1152 WRITE(LUPRI,'(//,1X,A)') 'Warning - Large differ' 1153 * //'ence between left and right excitation energies' 1154 WRITE(LUPRI,'(1X,A,I5)') 'Symmetry class:',ISYM 1155 WRITE(LUPRI,'(1X,A,I5)') 'Multiplicity :',IMULT 1156 WRITE(LUPRI,'(1X,A,I5)') 'State number :',ICOUNT 1157 WRITE(LUPRI,'(1X,A,2F20.10)') 'Eigenvalues :', 1158 * EIGVAL(ILSTNR),WORK(KEXCIS+ICOUNT-1) 1159 WRITE(LUPRI,'(1X,A,2F20.10)') 'Diff/Threshold:', 1160 * EIGVAL(ILSTNR)-WORK(KEXCIS+ICOUNT-1),THREIG 1161 WRITE(LUPRI,'(1X,A,//)') 1162 * 'Warning - Check your input and output carefully!!' 1163 END IF 1164 END DO 1165 1166 CALL CCEXNORM(ISYM,IMULT,THREIG,WORK,LWORK) 1167 1168 END IF 1169C 1170C==================================================================== 1171C-------------------------------------------------------------------- 1172C 1173C PERTURBATIVE CORRECTIONS TO CCS OR CCSD EXCITATION ENERGIES. 1174C 1175C-------------------------------------------------------------------- 1176C==================================================================== 1177C 1178 1200 CONTINUE 1179C 1180 IF ((LIST(1:2).EQ.'LE').AND. (.NOT.TRIPLET) .AND. 1181 * (.NOT.CCR12) .AND. 1182 * (CCP2.OR.((CCR3.OR.CCRT).OR.(CCR1A.OR.CCR1B)))) THEN 1183C 1184C--------------------------------------------------------------- 1185C If CCSDR(3) and first time (only calculating left CCSD) 1186C then we skip the rest. 1187C--------------------------------------------------------------- 1188C 1189 IF ( CCR3 .AND. ( .NOT. CCT)) GOTO 8000 1190C 1191 KE1 = KWRK1 1192 KE2 = KE1 + NCCEXCI(ISYM,IMULT) 1193 KE3 = KE2 + NCCEXCI(ISYM,IMULT) 1194 KE4 = KE3 + NCCEXCI(ISYM,IMULT) 1195 KWRK2 = KE4 + NCCEXCI(ISYM,IMULT) 1196 LWRK2 = LWORK - KWRK2 1197C 1198 CALL CC_PCEXCI(WORK(KEIVAL),WORK(KE1),WORK(KE2),WORK(KE3), 1199 * WORK(KE4),WORK(KWRK2),LWRK2) 1200C 1201 IF (CCP2 ) THEN 1202 WRITE (LUPRI,'(//,1X,A)') 1203 * 'Excitation energies with doubles corrections' 1204 ELSE IF (CCRT.OR.CCR3.OR.CCR1A.OR.CCR1B) THEN 1205 WRITE (LUPRI,'(//,1X,A)') 1206 * 'Excitation energies with triples corrections' 1207 ENDIF 1208C 1209 IF (CCP2) THEN 1210 WRITE (LUPRI,'(//A)') 1211 * ' CC-model Excitation no. Hartree eV' 1212 ENDIF 1213C 1214 IF (CCR3.OR.CCRT.OR.CCR1A.OR.CCR1B) THEN 1215 WRITE (LUPRI,'(//A)') 1216 * ' CC-model Excitation no. Hartree' 1217 * //' eV' 1218 ENDIF 1219C 1220C------------------------------------------------ 1221C Loop over excitations in symmetry class. 1222C------------------------------------------------ 1223C 1224 DO 1300 IST = 1,NCCEXCI(ISYM,IMULT) 1225C 1226 OME1 = WORK(KE1-1+IST) 1227 OME2 = WORK(KE2-1+IST) 1228 OME3 = WORK(KE3-1+IST) 1229 OME4 = WORK(KE4-1+IST) 1230C 1231C----------------------------------------- 1232C Write out Excitation energies. 1233C----------------------------------------- 1234C 1235 IF (CCP2) THEN 1236 WRITE (LUPRI,'(A)') 1237 * ' ---------------------------- -------- --' 1238 * //'------' 1239 WRITE (LUPRI,'(A16,I10,2F20.6)') ' CCS ', 1240 * IST,WORK(KEIVAL-1+IST),WORK(KEIVAL-1+IST)*XTEV 1241 WRITE (LUPRI,'(A)') 1242 * ' ---------------------------- -------- --' 1243 * //'------' 1244C 1245 WRITE (LUPRI,'(A16,I10,2F20.6)') ' C(0)A(0+2)C(0) ', 1246 * IST,OME1,OME1*XTEV 1247 WRITE (LUPRI,'(A16,I10,2F20.6)') ' C(0)A(1)C(1) ', 1248 * IST,OME2,OME2*XTEV 1249C 1250 OMECOR = OME1 + OME2 1251C 1252 WRITE (LUPRI,'(A)') 1253 * ' ---------------------------- -------- --' 1254 * //'------' 1255 WRITE (LUPRI,'(A16,I10,2F20.6)') ' CC(2) ', 1256 * IST,OMECOR,OMECOR*XTEV 1257 WRITE (LUPRI,'(A/)') 1258 * ' ============================ ======== ==' 1259 * //'======' 1260C 1261 ENDIF 1262C 1263 IF (CCR3.OR.CCRT.OR.CCR1A.OR.CCR1B) THEN 1264 ITST = 0 1265 DO 1310 IOM = 1, NOMINP(ISYMTR,IMULT) 1266 IF (IEX .EQ. IOMINP(IOM,ISYMTR,IMULT)) ITST = ITST + 1 1267 1310 CONTINUE 1268 IF ((ITST .EQ. 0 ).AND.CCSDT) GOTO 1300 1269C 1270 WRITE (LUPRI,'(A)') 1271 * ' -------- -------------- ----' 1272 * //'---- ---------' 1273 IF (.NOT. CCR3 ) THEN 1274 WRITE (LUPRI,'(A16,5X,I10,5X,2F20.6)')' CCSD ', 1275 * IST,WORK(KEIVAL-1+IST),WORK(KEIVAL-1+IST)*XTEV 1276 WRITE (LUPRI,'(A)') 1277 * ' -------- -------------- ----' 1278 * //'---- ---------' 1279 ELSE IF (IPRINT .GT. 10) THEN 1280 WRITE (LUPRI,'(A16,5X,I10,5X,2F20.6)')' CCSD*(T*ampl.) ', 1281 * IST,WORK(KEIVAL-1+IST),WORK(KEIVAL-1+IST)*XTEV 1282 WRITE (LUPRI,'(A)') 1283 * ' -------- -------------- ----' 1284 * //'---- ---------' 1285 ENDIF 1286C 1287 IF (CCR3) WRITE (LUPRI,'(A16,5X,I10,5X,2F20.6)') 1288 * ' CCSDR(3) ',IST,OME1,OME1*XTEV 1289 IF (CCRT) WRITE (LUPRI,'(A16,5X,I10,5X,2F20.6)') 1290 * ' CCSDR(T) ',IST,OME2,OME2*XTEV 1291 IF (CCR1A) WRITE (LUPRI,'(A16,5X,I10,5X,2F20.6)') 1292 * ' CCSDR(1a) ',IST,OME3,OME3*XTEV 1293 IF (CCR1B) WRITE (LUPRI,'(A16,5X,I10,5X,2F20.6)') 1294 * ' CCSDR(1b) ',IST,OME4,OME4*XTEV 1295C 1296 OMECOR = OME1 1297C 1298 WRITE (LUPRI,'(A/)') 1299 * ' -------- -------------- ----' 1300 * //'---- ---------' 1301C 1302 ENDIF 1303C 1304 IF (CCP2) WORK(KEXCPS+IST-1) = OMECOR 1305 IF (CCR3) WORK(KEXCPS+IST-1) = OME1 1306 IF (CCRT) WORK(KEXCPS2+IST-1) = OME2 1307 IF (CCR1A) WORK(KEXCPS3+IST-1) = OME3 1308 IF (CCR1B) WORK(KEXCPS4+IST-1) = OME4 1309C 1310C----------------------------------- 1311C End of loop over states. 1312C----------------------------------- 1313C 1314 1300 CONTINUE 1315C 1316 ENDIF 1317C 1318C---------------------------------------------- 1319C End of symmetry and multiplicity loop. 1320C---------------------------------------------- 1321C 1322 KEXCPS = KEXCPS + NCCEXSAV(ISYM,IMULT) 1323 KEXCPS2 = KEXCPS2 + NCCEXSAV(ISYM,IMULT) 1324 KEXCPS3 = KEXCPS3 + NCCEXSAV(ISYM,IMULT) 1325 KEXCPS4 = KEXCPS4 + NCCEXSAV(ISYM,IMULT) 1326 KEXCIS = KEXCIS + NCCEXSAV(ISYM,IMULT) 1327 KT1PS = KT1PS + NCCEXSAV(ISYM,IMULT) 1328C 1329 8000 CONTINUE 1330C 1331 9000 CONTINUE 1332C 1333 KEXCIS = KEXCI 1334 KT1PS = KT1P 1335C 1336 IF ((.NOT.CCR3).AND.(LHTR .OR. (.NOT.LIST(1:2).EQ.'LE'))) THEN 1337C 1338 WRITE(LURES,'(//A)') 1339 * ' +==============================================' 1340 * //'===============================+' 1341 WRITE(LURES,'(1X,A26,A10,A)') 1342 * '| sym. | Exci. | ',MODELP,' Excitation energies' 1343 * //' | ||T1|| |' 1344 WRITE(LURES,'(A)') 1345 * ' |(spin, | +-----------------------------' 1346 * //'-------------------------------+' 1347 WRITE(LURES,'(1X,A)') 1348 * '| spat) | | Hartree | eV. |' 1349 * //' cm-1 | % |' 1350 WRITE(LURES,'(A)') 1351 * ' +==============================================' 1352 * //'===============================+' 1353C 1354 DO 9100 ISYM = 1, NSYM 1355C 1356 DO 9050 IMULT = 1, 3, 2 1357C 1358 DO 9200 IEX = 1, NCCEXSAV(ISYM,IMULT) 1359C 1360 IF ( OMESC ) THEN 1361 ITST = 0 1362 DO 9210 IOM = 1, NOMINP(ISYM,IMULT) 1363 IF ((IOMINP(IOM,ISYM,IMULT).EQ.IEX).AND. 1364 * (.NOT.(OMEINP))) 1365 * EOMINP(IOM,ISYM,IMULT)=WORK(KEXCIS+IEX-1) 1366 IF (IEX .EQ. IOMINP(IOM,ISYM,IMULT)) ITST = ITST + 1 1367 9210 CONTINUE 1368 IF ((ITST .EQ. 0 ).AND.CCSDT) GOTO 9200 1369 ENDIF 1370C 1371c IF ( IEX .EQ. 1) THEN 1372 WRITE(LURES,9990) IMULT,REP(ISYM-1),IEX, 1373 * WORK(KEXCIS+IEX-1), 1374 * WORK(KEXCIS+IEX-1)*XTEV,WORK(KEXCIS+IEX-1)*XTKAYS, 1375 * WORK(KT1PS+IEX-1) 1376c ELSE 1377c WRITE(LURES,9991) IEX,WORK(KEXCIS+IEX-1), 1378c * WORK(KEXCIS+IEX-1)*XTEV,WORK(KEXCIS+IEX-1)*XTKAYS, 1379c * WORK(KT1PS+IEX-1) 1380c ENDIF 1381 OME = WORK(KEXCIS+IEX-1) 1382 ETOT = ECCGRS+OME 1383 CALL WRIPRO(OME,MOPRPC,0, 1384 * LABEL,LABEL,LABEL,LABEL, 1385 * ECCGRS,OME,ETOT,1,ISYM,IMULT,IEX) 1386 CALL WRIPRO(ETOT,MOPRPC,0, 1387 * LABEL1,LABEL1,LABEL1,LABEL1, 1388 * ECCGRS,OME,ETOT,1,ISYM,IMULT,IEX) 1389 9200 CONTINUE 1390C 1391 DO 9220 IEX = 1, NCCEXSAV(ISYM,IMULT) 1392 IF ((ISYM.EQ.IXSTSY).AND.(IEX.EQ.IXSTAT)) THEN 1393 OMECCX = WORK(KEXCIS+IEX-1) 1394 ECCXST = ECCGRS + OMECCX 1395 ENDIF 1396 9220 CONTINUE 1397C 1398 IF (.NOT.(NCCEXSAV(ISYM,IMULT).EQ.0)) THEN 1399 NREST = 0 1400 IF (IMULT.EQ.1) NREST = NCCEXCI(ISYM,3) 1401 DO ISYM2 = ISYM+1,NSYM 1402 DO IMULT2 = 1, 3, 2 1403 NREST = NREST + NCCEXCI(ISYM2,IMULT2) 1404 END DO 1405 END DO 1406 IF (NREST.EQ.0) GOTO 9100 1407 WRITE(LURES,'(A)') 1408 * ' +----------------------------------------------' 1409 * //'-------------------------------+' 1410 ENDIF 1411 KEXCIS = KEXCIS + NCCEXSAV(ISYM,IMULT) 1412 KT1PS = KT1PS + NCCEXSAV(ISYM,IMULT) 1413C 1414 9050 CONTINUE 1415 9100 CONTINUE 1416C 1417 WRITE(LURES,'(A)') 1418 * ' +==============================================' 1419 * //'===============================+' 1420C 1421 WRITE(LURES,'(//5X,A)') 'Total energies in Hartree:' 1422 KEXCIS = KEXCI 1423 KT1PS = KT1P 1424 DO ISYM = 1, NSYM 1425 DO IMULT = 1, 3, 2 1426 DO IEX = 1, NCCEXSAV(ISYM,IMULT) 1427 ITST = 1 1428 IF ( OMESC ) THEN 1429 ITST = 0 1430 DO IOM = 1, NOMINP(ISYM,IMULT) 1431 IF (IEX .EQ. IOMINP(IOM,ISYM,IMULT)) ITST = ITST + 1 1432 END DO 1433 ENDIF 1434 IF (.NOT.(ITST.EQ.0.AND.CCSDT)) THEN 1435 WRITE(LURES,9989) IEX,IMULT,REP(ISYM-1), 1436 * WORK(KEXCIS+IEX-1)+ECCGRS 1437 END IF 1438 END DO 1439 KEXCIS = KEXCIS + NCCEXSAV(ISYM,IMULT) 1440 KT1PS = KT1PS + NCCEXSAV(ISYM,IMULT) 1441 END DO 1442 END DO 1443 1444 ENDIF 1445C 1446 9989 FORMAT(10X,I4,' ^',I1,A3,F20.10) 1447 9990 FORMAT(1X,'| ^',I1,A3,' | ',I4,' | ',F13.7, 1448 * ' | ',F13.5,' | ',F13.3,' | ',F6.2,' |') 1449 9991 FORMAT(1X,'| | ',I4,' | ',F13.7, 1450 * ' | ',F13.5,' | ',F13.3,' | ',F6.2,' |') 1451C 1452 IF ((CCP2.OR.((CCR3.AND.CCT).OR.CCRT.OR.CCR1A.OR.CCR1B)) 1453 * .AND.(LIST(1:2).EQ.'LE')) THEN 1454C 1455 KEXCPS = KEXCP 1456 KEXCPS2 = KEXCP2 1457 KEXCPS3 = KEXCP3 1458 KEXCPS4 = KEXCP4 1459C 1460 WRITE(LURES,'(//A/)') 1461 * ' Excitation energies with perturbational corrections ' 1462C 1463C-------------------------------------- 1464C Loop over triples corrections. 1465C-------------------------------------- 1466C 1467 LTRIP(1) = CCR3 1468 LTRIP(2) = CCRT 1469 LTRIP(3) = CCR1A 1470 LTRIP(4) = CCR1B 1471C 1472 CCR3 = .FALSE. 1473 CCRT = .FALSE. 1474 CCR1A = .FALSE. 1475 CCR1B = .FALSE. 1476C 1477 IF ( CCP2 ) THEN 1478 NTRIP = 1 1479 ELSE 1480 NTRIP = 4 1481 ENDIF 1482C 1483 DO 9299 ITRIP = 1, NTRIP 1484C 1485 IF ((.NOT. LTRIP(ITRIP)).AND.(.NOT.CCP2)) GOTO 9299 1486C 1487 IF ((ITRIP.EQ.1).AND.(.NOT.CCP2)) CCR3 = .TRUE. 1488 IF (ITRIP.EQ.2) CCRT = .TRUE. 1489 IF (ITRIP.EQ.3) CCR1A = .TRUE. 1490 IF (ITRIP.EQ.4) CCR1B = .TRUE. 1491C 1492 IF (CCP2 ) MODEL1 = 'CIS(D)' 1493 IF (CCR3 ) MODEL1 = 'CCSDR(3)' 1494 IF (CCRT ) MODEL1 = 'CCSDR(T)' 1495 IF (CCR1A) MODEL1 = 'CCSDR(1a)' 1496 IF (CCR1B) MODEL1 = 'CCSDR(1b)' 1497 1498 IF (CCR12) THEN 1499 MODELSCR = MODEL1 1500 CALL CCSD_MODEL(MODEL1,LENMOD,24,MODELSCR,24,APROXR12) 1501 END IF 1502C 1503 WRITE(LURES,'(/A)') 1504 * ' +==============================================' 1505 * //'=====================+' 1506 WRITE(LURES,'(1X,A26,A10,A)') 1507 * '| sym. | Exci. | ',MODEL1,' Excitation energies' 1508 * //' |' 1509 WRITE(LURES,'(A)') 1510 * ' |(spin, | +-----------------------------' 1511 * //'---------------------+' 1512 WRITE(LURES,'(1X,A)') 1513 * '| spat) | | Hartree | eV. |' 1514 * //' cm-1 |' 1515 WRITE(LURES,'(A)') 1516 * ' +==============================================' 1517 * //'=====================+' 1518C 1519 DO 9300 ISYM = 1, NSYM 1520C 1521C DO 9310 IMULT = 1, 3, 2 1522 DO 9310 IMULT = 1, 1, 2 1523C 1524 DO 9400 IEX = 1, NCCEXCI(ISYM,IMULT) 1525C 1526 IF (CCP2) OME = WORK(KEXCPS +IEX-1) 1527 IF (CCR3) OME = WORK(KEXCPS +IEX-1) 1528 IF (CCRT) OME = WORK(KEXCPS2+IEX-1) 1529 IF (CCR1A) OME = WORK(KEXCPS3+IEX-1) 1530 IF (CCR1B) OME = WORK(KEXCPS4+IEX-1) 1531C 1532 ITST = 0 1533 DO 9410 IOM = 1, NOMINP(ISYM,IMULT) 1534 IF ((IOMINP(IOM,ISYM,IMULT).EQ.IEX).AND. 1535 * (.NOT.(OMEINP)).AND.OMESC) 1536 * EOMINP(IOM,ISYM,IMULT)=OME 1537 IF (IEX .EQ. IOMINP(IOM,ISYM,IMULT)) ITST=ITST+1 1538 9410 CONTINUE 1539C 1540 IF ((ITST .EQ. 0 ).AND.(.NOT.CCP2)) GOTO 9400 1541C 1542c IF ( IEX .EQ. 1) THEN 1543 WRITE(LURES,9992) IMULT,REP(ISYM-1),IEX,OME, 1544 * OME*XTEV,OME*XTKAYS 1545c ELSE 1546c WRITE(LURES,9993) IEX,OME, 1547c * OME*XTEV,OME*XTKAYS 1548c ENDIF 1549 ETOT = ECCGRS+OME 1550 CALL WRIPRO(OME,MODEL1,0, 1551 * LABEL,LABEL,LABEL,LABEL, 1552 * 0.0D0,0.0D0,0.0D0,1,ISYM,IMULT,IEX) 1553 CALL WRIPRO(ETOT,MODEL1,0, 1554 * LABEL1,LABEL1,LABEL1,LABEL1, 1555 * 0.0D0,0.0D0,0.0D0,1,ISYM,IMULT,IEX) 1556 1557 IF ((ISYM .EQ. IXSTSY).AND.(IEX.EQ.IXSTAT)) THEN 1558 OMECCX = OME 1559 ECCXST = ECCGRS + OMECCX 1560 ENDIF 1561C 1562 9400 CONTINUE 1563C 1564 IF (.NOT.((ISYM .EQ. NSYM).OR. 1565 * (NCCEXCI(ISYM,IMULT).EQ.0))) THEN 1566 NREST = 0 1567 DO 9350 ISYM2 = ISYM+1,NSYM 1568 NREST = NREST + NCCEXCI(ISYM2,IMULT) 1569 9350 CONTINUE 1570 IF (NREST.EQ.0) GOTO 9300 1571 WRITE(LURES,'(A)') 1572 * ' +----------------------------------------------' 1573 * //'---------------------+' 1574 ENDIF 1575C 1576 IF (CCP2) KEXCPS = KEXCPS + NCCEXCI(ISYM,IMULT) 1577 IF (CCR3) KEXCPS = KEXCPS + NCCEXCI(ISYM,IMULT) 1578 IF (CCRT) KEXCPS2 = KEXCPS2 + NCCEXCI(ISYM,IMULT) 1579 IF (CCR1A) KEXCPS3 = KEXCPS3 + NCCEXCI(ISYM,IMULT) 1580 IF (CCR1B) KEXCPS4 = KEXCPS4 + NCCEXCI(ISYM,IMULT) 1581C 1582 9310 CONTINUE 1583 9300 CONTINUE 1584C 1585 WRITE(LURES,'(A)') 1586 * ' +==============================================' 1587 * //'=====================+' 1588C 1589 KEXCPS = KEXCP 1590 KEXCPS2 = KEXCP2 1591 KEXCPS3 = KEXCP3 1592 KEXCPS4 = KEXCP4 1593 WRITE(LURES,'(//5X,A)') 'Total energies in Hartree:' 1594 DO ISYM = 1, NSYM 1595C DO IMULT = 1, 3, 2 1596 DO IMULT = 1, 1, 2 1597 DO IEX = 1, NCCEXCI(ISYM,IMULT) 1598 IF (CCP2) OME = WORK(KEXCPS +IEX-1) 1599 IF (CCR3) OME = WORK(KEXCPS +IEX-1) 1600 IF (CCRT) OME = WORK(KEXCPS2+IEX-1) 1601 IF (CCR1A) OME = WORK(KEXCPS3+IEX-1) 1602 IF (CCR1B) OME = WORK(KEXCPS4+IEX-1) 1603 ITST = 0 1604 DO IOM = 1, NOMINP(ISYM,IMULT) 1605 IF (IEX .EQ. IOMINP(IOM,ISYM,IMULT)) ITST=ITST+1 1606 END DO 1607 IF (.NOT.(ITST.EQ.0 .AND. .NOT.CCP2)) THEN 1608 WRITE(LURES,9989) IEX,IMULT,REP(ISYM-1), 1609 * OME+ECCGRS 1610 END IF 1611 END DO 1612 IF (CCP2) KEXCPS = KEXCPS + NCCEXCI(ISYM,IMULT) 1613 IF (CCR3) KEXCPS = KEXCPS + NCCEXCI(ISYM,IMULT) 1614 IF (CCRT) KEXCPS2 = KEXCPS2 + NCCEXCI(ISYM,IMULT) 1615 IF (CCR1A) KEXCPS3 = KEXCPS3 + NCCEXCI(ISYM,IMULT) 1616 IF (CCR1B) KEXCPS4 = KEXCPS4 + NCCEXCI(ISYM,IMULT) 1617 END DO 1618 END DO 1619C 1620 IF (CCR3) CCR3 = .FALSE. 1621 IF (CCRT) CCRT = .FALSE. 1622 IF (CCR1A) CCR1A = .FALSE. 1623 IF (CCR1B) CCR1B = .FALSE. 1624C 1625 9299 CONTINUE 1626C 1627 CCR3 = LTRIP(1) 1628 CCRT = LTRIP(2) 1629 CCR1A = LTRIP(3) 1630 CCR1B = LTRIP(4) 1631C 1632 9992 FORMAT(1X,'| ^',I1,A3,' | ',I4,' | ',F13.7, 1633 * ' | ',F13.5,' | ',F13.3,' |') 1634 9993 FORMAT(1X,'| | ',I4,' | ',F13.7, 1635 * ' | ',F13.5,' | ',F13.3,' |') 1636C 1637 ENDIF 1638C 1639 IF (IPRINT .GT.10) THEN 1640 CALL AROUND( ' END OF CC_EXCI ' ) 1641 ENDIF 1642C 1643 CALL QEXIT('CC_EXCI') 1644C 1645 RETURN 1646 END 1647c*DECK CCLR_LEINFI 1648 SUBROUTINE CCLR_LEINFI(TRIPLET) 1649C 1650C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 1651C Purpose: set common Leinf for calculation of 1652C Coupled Cluster excitation energies. 1653C 1654C Written by Ove Christiansen November 1994 1655C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 1656C 1657#include "implicit.h" 1658#include "priunit.h" 1659#include "ccorb.h" 1660#include "ccsdinp.h" 1661#include "cclr.h" 1662#include "ccsdsym.h" 1663#include "leinf.h" 1664Cholesky 1665#include "maxorb.h" 1666#include "ccdeco.h" 1667Cholesky 1668C 1669 LOGICAL TRIPLET 1670C 1671 CALL QENTER('CCLR_LEINFI') 1672C 1673C-------------------------- 1674C Set COMMON /LEINF/ 1675C-------------------------- 1676C 1677 IF (IPRINT .GT.10) THEN 1678 CALL AROUND( ' START OF CCLR_LEINFI ' ) 1679 ENDIF 1680C 1681 THRLE = THREXC 1682 MAXLE = MAXITE 1683 IF (CHOINT .AND. CHEXDI) MAXLE = MAXLE / 2 1684 IPRLE = IPRINT 1685C 1686Cholesky 1687C 1688 IF (CHOINT) THEN 1689 IF (TRIPLET) 1690 & CALL QUIT('CCLR_LEINFI: Triplet Cholesky not implemented') 1691 IF (CC2 .OR. CCS) THEN 1692 LETYPA = NT1AM(ISYMTR) 1693 ELSE 1694 CALL QUIT('CCLR_LEINFI: Only Cholesky CC2 implemented') 1695 ENDIF 1696 ELSE 1697 LETYPA = NT1AM(ISYMTR) + NT2AM(ISYMTR) 1698 IF ( TRIPLET ) LETYPA = LETYPA + NT2AMA(ISYMTR) 1699 IF ( CCR12 ) LETYPA = LETYPA + NTR12AM(ISYMTR) 1700 IF (CCS.OR.CCP2) LETYPA = NT1AM(ISYMTR) 1701 END IF 1702C 1703 LETYPB = 0 1704 LETOT = LETYPA + LETYPB 1705 NBASIS = 1 1706 NCREF = 0 1707 LULEA = 45 1708 LULEB = 46 1709 LULEC = 47 1710C 1711C-------------------------------------------------------------------- 1712C B space has zero dimension. 1713C A space : integral distribution, and vectors and intermediates. 1714C-------------------------------------------------------------------- 1715C 1716 LLEWB = 0 1717 LLEWA = 0 1718 NCCEX = 0 1719C 1720 DO 10 ISYM = 1, NSYM 1721 IF (NDISAO(ISYM).GT. LLEWA) LLEWA = NDISAO(ISYM) 1722 10 CONTINUE 1723C 1724 NRHO2 = 2*NT2ORT(1) + NT2SQ(1) 1725 IF (CCS) NRHO2 = 0 1726C 1727 LLEWA = LLEWA + NRHO2 1728C 1729 LLEWA = LLEWA + (10+NRHFT)*N2BST(1) 1730 LLEWA = LLEWA + 2*NT2BGD(1) 1731 LLEWA = LLEWA + NDSRHF(1) 1732 LLEWA = LLEWA + 2*NT2MAO(1,1) 1733C 1734 IF (CCSDT) LLEWA = LLEWA + NT2SQ(1) 1735C 1736Cholesky 1737 IF (CHOINT .AND. CC2) THEN 1738 NRHO2 = 0 1739 LLEWA = 0 1740 END IF 1741Cholesky 1742C 1743C--------------------------------------------------------------- 1744C t2tcor is taken care of in routine, that is; it is forced 1745C to be false if there is problems. 1746C--------------------------------------------------------------- 1747C 1748C IF (T2TCOR) LLEWA = LLEWA + NT2SQ(1) 1749C 1750 IF (IPRINT .GT.60) THEN 1751 CALL AROUND('COMMON CCLR in CCEXCI') 1752 WRITE (LUPRI,'(A32,F24.12)') 'IN CCLR_LEINFI: THREXC', THREXC 1753 IMULT = 1 1754 IF (TRIPLET) IMULT = 3 1755 WRITE(LUPRI,1) 'NCCEXCI:',(NCCEXCI(I,IMULT), I=1,NSYM) 1756 ENDIF 1757 IF (IPRINT .GT.60) THEN 1758 CALL AROUND('COMMON /LEINF/ in CCEXCI') 1759 WRITE (LUPRI,'(A32,F24.12)') 'IN CCLR_LEINFI: THRLE ', THRLE 1760 WRITE (LUPRI,'(A32,I12)') 'IN CCLR_LEINFI: MAXLE ', MAXLE 1761 WRITE (LUPRI,'(A32,I12)') 'IN CCLR_LEINFI: IPRLE ', IPRLE 1762 WRITE (LUPRI,'(A32,I12)') 'IN CCLR_LEINFI: LETYPA', LETYPA 1763 WRITE (LUPRI,'(A32,I12)') 'IN CCLR_LEINFI: LETYPB', LETYPB 1764 WRITE (LUPRI,'(A32,I12)') 'IN CCLR_LEINFI: LLEWA ', LLEWA 1765 WRITE (LUPRI,'(A32,I12)') 'IN CCLR_LEINFI: LLEWB ', LLEWB 1766 WRITE (LUPRI,'(A32,I12)') 'IN CCLR_LEINFI: LETOT ', LETOT 1767 WRITE (LUPRI,'(A32,I12)') 'IN CCLR_LEINFI: NSIDE ', NSIDE 1768 WRITE (LUPRI,'(A32,I12)') 'IN CCLR_LEINFI: NCREF ', NCREF 1769 WRITE (LUPRI,'(A32,I12)') 'IN CCLR_LEINFI: LULEA ', LULEA 1770 WRITE (LUPRI,'(A32,I12)') 'IN CCLR_LEINFI: LULEB ', LULEB 1771 WRITE (LUPRI,'(A32,I12)') 'IN CCLR_LEINFI: LULEC ', LULEC 1772 END IF 1773C 1774 IF (IPRINT .GT.10) THEN 1775 CALL AROUND( ' END OF CCLR_LEINFI ' ) 1776 ENDIF 1777C 1778 1 FORMAT(3X,A8,8I8) 1779C 1780 CALL QEXIT('CCLR_LEINFI') 1781C 1782 END 1783c*DECK CC_PCEXCI 1784 SUBROUTINE CC_PCEXCI(OMES,OME1,OME2,OME3,OME4,WORK,LWORK) 1785C 1786C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 1787C 1788C Purpose: 1789C Direct calculation of perturbative corrections 1790C to Coupled Cluster excitation energies. 1791C 1792C CC(2) correction to CCS - identical to CIS(D) 1793C (CCS =TDA=CIS ),(CC(2)=CIS(D)) 1794C 1795C CCSDR(3),CCSDR(1a),CCSDR(1b),CCSDR(T) 1796C 1797C Written by Ove Christiansen 10-3-1995 1798C 8-10-1995 1799C 26-2-1996 1800C 10-12-1996 1801C 3-2-1997 1802C Version 3, February 1996 1803C 1804C 1805C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 1806C 1807#include "implicit.h" 1808#include "priunit.h" 1809#include "maxorb.h" 1810#include "ccorb.h" 1811#include "iratdef.h" 1812#include "cclr.h" 1813#include "ccsdsym.h" 1814#include "ccsdio.h" 1815#include "ccsdinp.h" 1816#include "ccexci.h" 1817C 1818 PARAMETER(NTRIP = 4,MXISPT = 3) 1819C 1820 LOGICAL TRIPLET 1821 LOGICAL CCSSAV, CC2SAV, CC3SAV, LINQCC, LTRIP(NTRIP),LHSOLD 1822 DIMENSION OME1(*),OME2(*),OME3(*),OME4(*),WORK(LWORK),OMES(*) 1823 DIMENSION IADR(MXISPT) 1824C 1825 CHARACTER CHSYM*4,MODEL1*10,MODEL*10,APROXR12*3 1826 CHARACTER LBLPT(MXISPT)*8 1827 CHARACTER*8 FC1AM,FC2AM,FRHO1,FRHO2 1828 PARAMETER (TWO = 2.0D00,TF=1.0D-04, TLOVLP = 0.5D0 ) 1829 PARAMETER (FC1AM='CCR_C1AM',FC2AM='CCR_C2AM') 1830 PARAMETER (FRHO1='CCR_RHO1',FRHO2='CCR_RHO2') 1831C 1832#include "leinf.h" 1833C 1834 CALL QENTER('CC_PCEXCI') 1835C 1836 LUFC1 = -1 1837 LUFC2 = -1 1838 LUFR1 = -1 1839 LUFR2 = -1 1840C 1841 WRITE (LUPRI,'(1X,A,/)') ' ' 1842 WRITE (LUPRI,'(1X,A)') 1843 *'*********************************************************'// 1844 *'**********' 1845 WRITE (LUPRI,'(1X,A)') 1846 *'* '// 1847 *' *' 1848 IF ( CCR1A.OR.CCR1B.OR.CCRT.OR.CCR3.OR.CCP2) THEN 1849 WRITE (LUPRI,'(1X,A)') 1850 * '* '// 1851 * ' *' 1852 WRITE (LUPRI,'(1X,A)') 1853 * '*---------- CALCULATION OF PERTURBATIONAL CORRECTIONS >'// 1854 * '---------*' 1855 ENDIF 1856 WRITE (LUPRI,'(1X,A)') 1857 *'* '// 1858 *' *' 1859 WRITE (LUPRI,'(1X,A,/)') 1860 *'*********************************************************'// 1861 *'**********' 1862C 1863 TRIPLET = .FALSE. 1864 IMULT = 1 1865C 1866 MODEL = 'CCSD ' 1867 IF (CCS) MODEL = 'CC2 ' 1868 IF (CCP2) MODEL = 'CC(2) ' 1869 IF (CCR1A) MODEL = 'CCSDR(1A) ' 1870 IF (CCR1B) MODEL = 'CCSDR(1B) ' 1871 IF (CCRT) MODEL = 'CCSDR(T) ' 1872 IF (CCR3) MODEL = 'CCSDR(3) ' 1873C 1874 IF ( CCR3.OR.CCRT.OR.CCR1A.OR.CCR1B ) THEN 1875 CC3SAV = CCSDT 1876 CCSDT = .FALSE. 1877 ENDIF 1878C 1879 IF ( CCP2 ) THEN 1880 CCSSAV = CCS 1881 CC2SAV = CC2 1882 CCS = .FALSE. 1883 CC2 = .TRUE. 1884 ENDIF 1885C 1886 IF (CCP2.OR.((CCR3.OR.CCRT).OR.(CCR1A.OR.CCR1B))) THEN 1887C 1888 WRITE (LUPRI,'(//A)') 1889 * ' ========================================================' 1890 * //'==============' 1891 WRITE (LUPRI,'(A)') 1892 * ' ### ' 1893 * //' ###' 1894 WRITE (LUPRI,'(A)') 1895 * ' ### ' 1896 * //' ###' 1897 WRITE (LUPRI,'(A)') 1898 * ' ### Perturbational Corrections to Excitation ' 1899 * //' energies. ###' 1900 WRITE (LUPRI,'(A)') 1901 * ' ### ' 1902 * //' ###' 1903 WRITE (LUPRI,'(A)') 1904 * ' ### ' 1905 * //' ###' 1906 IF (CCR3.OR.CCRT.OR.CCR1A.OR.CCR1B) THEN 1907 WRITE(LUPRI,'(A)') 1908 * ' ### Calculating triples corrected CCSD excitation e' 1909 * //'nergies. ###' 1910 ENDIF 1911 IF (CCP2) THEN 1912 WRITE(LUPRI,'(A)') 1913 * ' ### Calculating doubles corrected CCS excitation e' 1914 * //'nergies. ###' 1915 ENDIF 1916 WRITE (LUPRI,'(A)') 1917 * ' ### ' 1918 * //' ###' 1919 WRITE (LUPRI,'(A)') 1920 * ' ### ' 1921 * //' ###' 1922 WRITE (LUPRI,'(A)') 1923 * ' ========================================================' 1924 * //'==============' 1925C 1926 ENDIF 1927C 1928 IF (CCR3.OR.CCRT.OR.CCR1A.OR.CCR1B ) THEN 1929C 1930C-------------------------- 1931C Work allocation 2. 1932C-------------------------- 1933C 1934 NTAMP = NT1AM(ISYMTR) + NT2AM(ISYMTR) 1935C 1936 KRV1 = 1 1937 KEND1 = KRV1 + NTAMP 1938 LEND1 = LWORK - KEND1 1939C 1940 IF (LEND1.LE.0) CALL QUIT('Too little work space in cc_pc') 1941C 1942 DO 50 IV = 1, NCCEXCI(ISYMTR,IMULT) 1943C 1944 ISTATE = ISYOFE(ISYMTR) + IV 1945 EIGV = EIGVAL(ISTATE) 1946 IF (IPRINT .GT. 5) THEN 1947 WRITE(LUPRI,'(/,1x,A,I3,/1X,A,I3,A,F16.8)') 1948 * 'Calculating triples corrections for state ',ISTATE, 1949 * 'of symmetry ',ISYMTR,' and eigenvalue: ',EIGV 1950 ENDIF 1951C 1952C------------------------------------------------------------ 1953C Calculate (A(0)+A(1)*(A(0)-E(ccsd)1)-1*A(1))*C(0) 1954C Loop over various triple correction approaches. 1955C------------------------------------------------------------ 1956C 1957 ECURR = EIGV 1958C 1959 ITST = 0 1960 DO 55 IOM = 1, NOMINP(ISYMTR,IMULT) 1961 IF (IOMINP(IOM,ISYMTR,IMULT).EQ.IV) ITST = ITST + 1 1962 55 CONTINUE 1963C 1964 IF (ITST .EQ. 0 ) THEN 1965 IF (CCR3) OME1(IV) = 0.0 1966 IF (CCRT) OME2(IV) = 0.0 1967 IF (CCR1A) OME3(IV) = 0.0 1968 IF (CCR1B) OME4(IV) = 0.0 1969 GOTO 50 1970 ENDIF 1971C 1972 IF ( IPRINT. GT. 5) THEN 1973 WRITE(LUPRI,'(A,F10.6)') ' Doing partioned triples ' 1974 * //'linear transformation with ECURR= ',ECURR 1975 ENDIF 1976C 1977C-------------------------------------------------- 1978C loop over different triple corrections. 1979C-------------------------------------------------- 1980C 1981 LTRIP(1) = CCR3 1982 LTRIP(2) = CCRT 1983 LTRIP(3) = CCR1A 1984 LTRIP(4) = CCR1B 1985C 1986 CCR3 = .FALSE. 1987 CCRT = .FALSE. 1988 CCR1A = .FALSE. 1989 CCR1B = .FALSE. 1990C 1991 DO 60 ITRIP = 1, NTRIP 1992C 1993 IF (.NOT. LTRIP(ITRIP)) GOTO 60 1994C 1995 CCSDT = .TRUE. 1996 IF (ITRIP.EQ.1) CC3 = .TRUE. 1997 IF (ITRIP.EQ.2) CCRT = .TRUE. 1998 IF (ITRIP.EQ.3) CC1A = .TRUE. 1999 IF (ITRIP.EQ.4) CC1B = .TRUE. 2000C 2001 IF ( IPRINT .GT. 0 ) THEN 2002 WRITE(LUPRI,'(/,A,F10.6)') ' ECURR: ',ECURR 2003 WRITE(LUPRI,'(A,F10.6)') ' L*R Norm is ', 2004 & XNORM(ISTATE)**2 2005 ENDIF 2006 IF ( IPRINT .GT. 5 ) THEN 2007 WRITE(LUPRI,*) 'CC3,CC1A,CC1B,CCRT,CC3LR: ', 2008 * CC3,CC1A,CC1B,CCRT,CC3LR 2009 ENDIF 2010C 2011C------------------------------------------------------------ 2012C Readin right solution and find excitation energy. 2013C------------------------------------------------------------ 2014C 2015 IOPT = 3 2016 CALL CC_RDRSP('RE',ISTATE,ISYMTR,IOPT,MODEL,WORK(KRV1), 2017 * WORK(KRV1+NT1AM(ISYMTR))) 2018 IF (IPRINT .GT. 40 ) THEN 2019 CALL AROUND( 'In CC_PCEXCI: right eigen vector ' ) 2020 CALL CC_PRP(WORK(KRV1),WORK(KRV1+NT1AM(ISYMTR)), 2021 * ISYMTR,1,1) 2022 ENDIF 2023C 2024C-------------------------------------------------------------------- 2025C Calculate linear transformation. 2026C Input vector is first elements in work - 2027C Output vector replaces vector as first elements in work. 2028C-------------------------------------------------------------------- 2029C 2030 CALL CC_ATRR(ECURR,ISYMTR,1,WORK,LWORK,.FALSE.,DUMMY, 2031 & APROXR12,.FALSE.) 2032C 2033 IF (IPRINT .GT.25) THEN 2034 CALL AROUND('CC_PCEXCI: partioned A*R(0) vector. ') 2035 CALL CC_PRP(WORK(KRV1),WORK(KRV1+NT1AM(ISYMTR)), 2036 * ISYMTR,1,1) 2037 ENDIF 2038C 2039C-------------------------------- 2040C Readin left solution. 2041C-------------------------------- 2042C 2043 KLV1 = KEND1 2044 KEND2 = KLV1 + NTAMP 2045 LEND2 = LWORK - KEND2 2046 IF (LEND2.LE.0) CALL QUIT('Too little work '// 2047 & 'space in cc_pc-2') 2048C 2049 IOPT = 3 2050 CALL CC_RDRSP('LE',ISTATE,ISYMTR,IOPT,MODEL,WORK(KLV1), 2051 * WORK(KLV1+NT1AM(ISYMTR))) 2052 IF (IPRINT .GT. 40 ) THEN 2053 CALL AROUND( 'In CC_PCEXCI: Left Eigen vector ' ) 2054 CALL CC_PRP(WORK(KLV1),WORK(KLV1+NT1AM(ISYMTR)), 2055 * ISYMTR,1,1) 2056 ENDIF 2057C 2058 IF (CC3) OME1(IV) = 2059 * DDOT(NTAMP,WORK(KLV1),1,WORK(KRV1),1)/XNORM(ISTATE)**2 2060 IF (CCRT) OME2(IV) = 2061 * DDOT(NTAMP,WORK(KLV1),1,WORK(KRV1),1)/XNORM(ISTATE)**2 2062 IF (CC1A) OME3(IV) = 2063 * DDOT(NTAMP,WORK(KLV1),1,WORK(KRV1),1)/XNORM(ISTATE)**2 2064 IF (CC1B) OME4(IV) = 2065 * DDOT(NTAMP,WORK(KLV1),1,WORK(KRV1),1)/XNORM(ISTATE)**2 2066C 2067 IF ( IPRINT .GT. 1) THEN 2068 IF (CC3) WRITE(LUPRI,'(1X,A9,I3,1X,A15,1X,F10.6)') 2069 * 'Exci. nr.',IV,'in CCSDR(3) is',OME1(IV) 2070 IF (CCRT) WRITE(LUPRI,'(1X,A9,I3,1X,A15,1X,F10.6)') 2071 * 'Exci. nr.',IV,'in CCSDR(T) is',OME2(IV) 2072 IF (CC1A) WRITE(LUPRI,'(1X,A9,I3,1X,A15,1X,F10.6)') 2073 * 'Exci. nr.',IV,'in CCSDR(1A) is',OME3(IV) 2074 IF (CC1B) WRITE(LUPRI,'(1X,A9,I3,1X,A15,1X,F10.6)') 2075 * 'Exci. nr.',IV,'in CCSDR(1B) is',OME4(IV) 2076 ENDIF 2077 CCSDT = .FALSE. 2078 IF (CC3) CC3 = .FALSE. 2079 IF (CCRT) CCRT = .FALSE. 2080 IF (CC1A) CC1A = .FALSE. 2081 IF (CC1B) CC1B = .FALSE. 2082C 2083 60 CONTINUE 2084C 2085 CCR3 = LTRIP(1) 2086 CCRT = LTRIP(2) 2087 CCR1A = LTRIP(3) 2088 CCR1B = LTRIP(4) 2089C 2090 50 CONTINUE 2091C 2092 ENDIF 2093C 2094C------------------- 2095C CC(2) section. 2096C NB!! kan laves smatere ved at lave logisk flag 2097C saa kun C1 transformeres altsaa spare allokering 2098C til C2SQ. 2099C------------------- 2100C 2101 IF (CCP2) THEN 2102C 2103C-------------------------- 2104C Work allocation 3. 2105C-------------------------- 2106C 2107 NTAMP = NT1AM(ISYMTR) + NT2AM(ISYMTR) 2108 NT = NT1AM(ISYMTR) 2109C 2110 DO 100 IV = 1, NCCEXCI(ISYMTR,IMULT) 2111C 2112 ISTATE = ISYOFE(ISYMTR) + IV 2113 EIGV = EIGVAL(ISTATE) 2114 IF (IPRINT .GT. 5) THEN 2115 WRITE(LUPRI,'(/,1x,A,I3,/1X,A,I3,A,F16.8)') 2116 * 'Calculating doubles corrections for state ',ISTATE, 2117 * 'of symmetry ',ISYMTR,' and eigenvalue: ',EIGV 2118 ENDIF 2119C 2120 KRV1 = 1 2121 KEND1 = KRV1 + NTAMP 2122 LEND1 = LWORK - KEND1 2123 IF (LEND1.LE.0) CALL QUIT('Too little work space in cc_pc') 2124 CALL DZERO(WORK(KRV1),NTAMP) 2125C 2126 IOPT = 1 2127 CALL CC_RDRSP('RE',ISTATE,ISYMTR,IOPT,MODEL,WORK(KRV1), 2128 * WORK(KRV1+NT1AM(ISYMTR))) 2129 IF (IPRINT .GT.15) THEN 2130 CALL AROUND('CC_PCEXCI: C(0) vector. ') 2131 CALL CC_PRP(WORK(KRV1),WORK(KRV1+NT1AM(ISYMTR)), 2132 * ISYMTR,1,1) 2133 ENDIF 2134C 2135 CALL CC_ATRR(ECURR,ISYMTR,1,WORK,LWORK,.FALSE.,DUMMY, 2136 & APROXR12,.FALSE.) 2137 IF (IPRINT .GT.15) THEN 2138 CALL AROUND('CC_PCEXCI: A*C(0) vector. ') 2139 CALL CC_PRP(WORK(KRV1),WORK(KRV1+NT1AM(ISYMTR)), 2140 * ISYMTR,1,1) 2141 ENDIF 2142C 2143C----------------------------------------------------------------------- 2144C Scale Vector with orbital energy diff. and exci. energy. 2145C----------------------------------------------------------------------- 2146C 2147 KLV1 = KEND1 2148 KEND2 = KLV1 + NT 2149 LEND2 = LWORK - KEND2 2150 IF (LEND2.LE.0) CALL QUIT('Too little work space in cc_pc-2') 2151 IOPT = 1 2152 CALL CC_RDRSP('LE',ISTATE,ISYMTR,IOPT,MODEL,WORK(KLV1), 2153 * WORK(KRV1+NT1AM(ISYMTR))) 2154 OME1(IV) = DDOT(NT1AM(ISYMTR),WORK(KLV1),1,WORK(KRV1),1) 2155 CALL CC_OMEC(OME2(IV),WORK(KRV1+NT1AM(ISYMTR)),EIGV, 2156 * WORK(KEND1),LEND1,ISYMTR) 2157 100 CONTINUE 2158C 2159 ENDIF 2160C 2161 IF (CCP2) THEN 2162C 2163 CCS = CCSSAV 2164 CC2 = CC2SAV 2165C 2166 ENDIF 2167C 2168 IF ( CCR3.OR.CCRT.OR.CCR1A.OR.CCR1B ) THEN 2169 CC3SAV = .FALSE. 2170 CCSDT = CC3SAV 2171 ENDIF 2172C 2173 CALL QEXIT('CC_PCEXCI') 2174C 2175 END 2176C /* Deck cc_vscal */ 2177 SUBROUTINE CC_VSCAL(OMEGA1,OMEGA2,OME,WORK,LWORK,ISYMTR) 2178C 2179C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 2180C Ove Christiansen 10-5-1995 2181C 2182C Purpose: Scale OMEGA with diagonals. 2183C eps(a) - eps(i) 2184C eps(a) + eps(b) - eps(i) - eps(j) 2185C Used for calculating pert. corr. amplitudes in 2186C CCSDR(3) 2187C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 2188C 2189#include "implicit.h" 2190#include "dummy.h" 2191C 2192 DIMENSION OMEGA1(*),OMEGA2(*),WORK(LWORK) 2193C 2194#include "priunit.h" 2195#include "inftap.h" 2196#include "ccorb.h" 2197#include "ccsdsym.h" 2198#include "ccsdinp.h" 2199Cholesky 2200#include "maxorb.h" 2201#include "ccdeco.h" 2202Cholesky 2203C 2204 INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 2205C 2206 CALL QENTER('CC_VSCAL') 2207C 2208C----------------------- 2209C Memory allocation. 2210C----------------------- 2211C 2212 KSCR1 = 1 2213 KEND = KSCR1 + NORBTS 2214 LWRK = LWORK - KEND 2215C 2216 IF (LWRK .LT. 0) THEN 2217 CALL QUIT('Insufficient space in CC2_FCK') 2218 ENDIF 2219C 2220C------------------------------------- 2221C Read canonical orbital energies. 2222C------------------------------------- 2223C 2224 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY, 2225 & .FALSE.) 2226 REWIND LUSIFC 2227C 2228 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 2229 READ (LUSIFC) 2230 READ (LUSIFC) (WORK(I), I=1,NORBTS) 2231C 2232 CALL GPCLOSE(LUSIFC,'KEEP') 2233C 2234 IF (FROIMP .OR. FROEXP) 2235 * CALL CCSD_DELFRO(WORK(KSCR1),WORK(KEND),LWRK) 2236C 2237 IF (IPRINT .GT. 80) THEN 2238 CALL AROUND('CC_VSCAL - Orbital energies. ') 2239 WRITE (LUPRI,*) (WORK(I), I=1,NORBTS) 2240 CALL AROUND('CC_VSCAL - start - : RHO1,RHO2 ') 2241 CALL CC_PRP(OMEGA1,OMEGA2,ISYMTR,1,1) 2242 ENDIF 2243C 2244C---------------------- 2245C Transform vector. 2246C---------------------- 2247C 2248 DO 10 ISYMI = 1, NSYM 2249C 2250 ISYMA = MULD2H(ISYMI,ISYMTR) 2251C 2252 DO 20 I = 1, NRHF(ISYMI) 2253C 2254 MI = IORB(ISYMI) + I 2255C 2256 DO 30 A = 1, NVIR(ISYMA) 2257C 2258 MA = IORB(ISYMA) + NRHF(ISYMA) + A 2259C 2260 NAI = IT1AM(ISYMA,ISYMI) 2261 * + NVIR(ISYMA)*(I - 1) + A 2262C 2263 DEN = (WORK(MA) - WORK(MI) - OME ) 2264C 2265 XIDEN = 1/DEN 2266C 2267 OMEGA1(NAI) = OMEGA1(NAI)*XIDEN 2268C 2269 30 CONTINUE 2270 20 CONTINUE 2271 10 CONTINUE 2272C 2273C 2274C Skip doubles part if CCS or Cholesky CC2 2275C 2276 IF (.NOT. (CCS .OR. (CHOINT .AND. CC2))) THEN 2277 DO 100 ISYMBJ = 1,NSYM 2278C 2279 ISYMAI = MULD2H(ISYMBJ,ISYMTR) 2280C 2281 DO 110 ISYMJ = 1,NSYM 2282C 2283 ISYMB = MULD2H(ISYMJ,ISYMBJ) 2284C 2285 DO 120 ISYMI = 1,NSYM 2286C 2287 ISYMA = MULD2H(ISYMI,ISYMAI) 2288C 2289 DO 130 J = 1,NRHF(ISYMJ) 2290C 2291 MJ = IORB(ISYMJ) + J 2292C 2293 DO 140 B = 1,NVIR(ISYMB) 2294C 2295 NBJ = IT1AM(ISYMB,ISYMJ) 2296 * + NVIR(ISYMB)*(J - 1) + B 2297C 2298 MB = IORB(ISYMB) + NRHF(ISYMB) + B 2299C 2300 DO 150 I = 1,NRHF(ISYMI) 2301C 2302 MI = IORB(ISYMI) + I 2303C 2304 DO 160 A = 1,NVIR(ISYMA) 2305C 2306 NAI = IT1AM(ISYMA,ISYMI) 2307 * + NVIR(ISYMA)*(I - 1) + A 2308C 2309 MA = IORB(ISYMA) + NRHF(ISYMA) + A 2310C 2311 IF (((ISYMAI.EQ.ISYMBJ).AND. 2312 * (NAI .GT. NBJ)).OR.(ISYMAI.GT.ISYMBJ)) 2313 * GOTO 160 2314C 2315 IF (ISYMAI.EQ.ISYMBJ) THEN 2316 NAIBJ = IT2AM(ISYMAI,ISYMBJ) 2317 * + INDEX(NAI,NBJ) 2318 ELSE 2319 NAIBJ = IT2AM(ISYMAI,ISYMBJ) 2320 * + NT1AM(ISYMAI)*(NBJ-1) + NAI 2321 ENDIF 2322C 2323 DEN = (WORK(MA) + WORK(MB) 2324 * - WORK(MI) - WORK(MJ) - OME ) 2325C 2326 XIDEN = 1/DEN 2327C 2328 OMEGA2(NAIBJ) = OMEGA2(NAIBJ)*XIDEN 2329C 2330 160 CONTINUE 2331 150 CONTINUE 2332 140 CONTINUE 2333 130 CONTINUE 2334 120 CONTINUE 2335 110 CONTINUE 2336 100 CONTINUE 2337 ENDIF 2338C 2339 IF (IPRINT .GT. 80) THEN 2340 IF (CCS .OR. (CHOINT.AND.CC2)) THEN 2341 I = 0 2342 CALL AROUND('CC_VSCAL - end - : RHO1') 2343 ELSE 2344 I = 1 2345 CALL AROUND('CC_VSCAL - end - : RHO1,RHO2') 2346 ENDIF 2347 CALL CC_PRP(OMEGA1,OMEGA2,ISYMTR,1,I) 2348 ENDIF 2349C 2350 CALL QEXIT('CC_VSCAL') 2351C 2352 RETURN 2353 END 2354C /* Deck cc_omec */ 2355 SUBROUTINE CC_OMEC(OMEC,OMEGA2,OME,WORK,LWORK,ISYMTR) 2356C 2357C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 2358C Ove Christiansen 10-5-1995 2359C 2360C Purpose: Used in CC(2) for scaling 2361C in this case OME is different 2362C from zero. 2363C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 2364C 2365#include "implicit.h" 2366#include "priunit.h" 2367#include "dummy.h" 2368C 2369 DIMENSION OMEGA2(*),WORK(LWORK) 2370C 2371#include "inftap.h" 2372#include "ccorb.h" 2373#include "ccsdsym.h" 2374#include "ccsdinp.h" 2375C 2376 INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 2377C 2378 CALL QENTER('CC_OMEC') 2379C 2380C----------------------- 2381C Memory allocation. 2382C----------------------- 2383C 2384 KSCR1 = 1 2385 KEND = KSCR1 + NORBTS 2386 LWRK = LWORK - KEND 2387C 2388 IF (LWRK .LT. 0) THEN 2389 CALL QUIT('Insufficient space in CC_OMEC') 2390 ENDIF 2391C 2392C------------------------------------- 2393C Read canonical orbital energies. 2394C------------------------------------- 2395C 2396 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY, 2397 & .FALSE.) 2398 REWIND LUSIFC 2399C 2400 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 2401 READ (LUSIFC) 2402 READ (LUSIFC) (WORK(I), I=1,NORBTS) 2403C 2404 CALL GPCLOSE(LUSIFC,'KEEP') 2405C 2406 IF (FROIMP .OR. FROEXP) 2407 * CALL CCSD_DELFRO(WORK(KSCR1),WORK(KEND),LWRK) 2408C 2409C---------------------------- 2410C Calculate contribution. 2411C---------------------------- 2412C 2413 OMEC = 0.0D00 2414C 2415 DO 100 ISYMBJ = 1,NSYM 2416C 2417 ISYMAI = MULD2H(ISYMBJ,ISYMTR) 2418C 2419 DO 110 ISYMJ = 1,NSYM 2420C 2421 ISYMB = MULD2H(ISYMJ,ISYMBJ) 2422C 2423 DO 120 ISYMI = 1,NSYM 2424C 2425 ISYMA = MULD2H(ISYMI,ISYMAI) 2426 ISYMAJ = MULD2H(ISYMA,ISYMJ) 2427 ISYMBI = MULD2H(ISYMB,ISYMI) 2428C 2429 DO 130 J = 1,NRHF(ISYMJ) 2430C 2431 MJ = IORB(ISYMJ) + J 2432C 2433 DO 140 B = 1,NVIR(ISYMB) 2434C 2435 NBJ = IT1AM(ISYMB,ISYMJ) 2436 * + NVIR(ISYMB)*(J - 1) + B 2437C 2438 MB = IORB(ISYMB) + NRHF(ISYMB) + B 2439C 2440 DO 150 I = 1,NRHF(ISYMI) 2441C 2442 NBI = IT1AM(ISYMB,ISYMI) 2443 * + NVIR(ISYMB)*(I - 1) + B 2444C 2445 MI = IORB(ISYMI) + I 2446C 2447 DO 160 A = 1,NVIR(ISYMA) 2448C 2449 NAI = IT1AM(ISYMA,ISYMI) 2450 * + NVIR(ISYMA)*(I - 1) + A 2451C 2452 NAJ = IT1AM(ISYMA,ISYMJ) 2453 * + NVIR(ISYMA)*(J - 1) + A 2454C 2455 MA = IORB(ISYMA) + NRHF(ISYMA) + A 2456C 2457 IF (((ISYMAI.EQ.ISYMBJ).AND. 2458 * (NAI .GT. NBJ)).OR.(ISYMAI.GT.ISYMBJ)) 2459 * GOTO 160 2460C 2461 IF (ISYMAI.EQ.ISYMBJ) THEN 2462 NAIBJ = IT2AM(ISYMAI,ISYMBJ) 2463 * + INDEX(NAI,NBJ) 2464 ELSE 2465 NAIBJ = IT2AM(ISYMAI,ISYMBJ) 2466 * + NT1AM(ISYMAI)*(NBJ-1) + NAI 2467 ENDIF 2468C 2469 IF (ISYMAJ .EQ. ISYMBI) THEN 2470C 2471 NAJBI = IT2AM(ISYMAJ,ISYMBI) 2472 * + INDEX(NAJ,NBI) 2473C 2474 ELSE IF (ISYMAJ .LT. ISYMBI) THEN 2475C 2476 NAJBI = IT2AM(ISYMAJ,ISYMBI) 2477 * + NT1AM(ISYMAJ)*(NBI - 1) + NAJ 2478C 2479 ELSE IF (ISYMBI .LT. ISYMAJ) THEN 2480C 2481 NAJBI = IT2AM(ISYMAJ,ISYMBI) 2482 * + NT1AM(ISYMBI)*(NAJ - 1) + NBI 2483C 2484 ENDIF 2485C 2486 DEN = - ( WORK(MA) + WORK(MB) 2487 * - WORK(MI) - WORK(MJ) - OME ) 2488C 2489 XIDEN = 1/DEN 2490 XIDEN2 = XIDEN 2491C 2492 IF (ISYMAI.EQ.ISYMBJ) THEN 2493 IF (NAI.EQ.NBJ) XIDEN2 = XIDEN2*2 2494 ENDIF 2495C 2496 OMEC = OMEC + 2497 * (2*OMEGA2(NAIBJ)-OMEGA2(NAJBI))* 2498 * OMEGA2(NAIBJ)*XIDEN2 2499C 2500 160 CONTINUE 2501 150 CONTINUE 2502 140 CONTINUE 2503 130 CONTINUE 2504 120 CONTINUE 2505 110 CONTINUE 2506 100 CONTINUE 2507C 2508 CALL QEXIT('CC_OMEC') 2509C 2510 RETURN 2511 END 2512c*DECK CC_ATRR 2513 SUBROUTINE CC_ATRR(ECURR,ISYMV,ISIDE,WORK,LWORK,LCONT,CONT, 2514 & APROXR12,TRIPLET) 2515C 2516C---------------------------------------------------------------------- 2517C Ove Christiansen December 1996. 2518C 2519C Jacobian transformation with ONE vector. 2520C Vector is first element in WORK and on output is replaced 2521C by its linear transformed. 2522C For CCR12 on input R_12 is passed after the conventional trial vector 2523C and while on output A * R_12 and S * R_12 passed at this place 2524C 2525C removed problem with CCS left transformation, C.H., October 1997 2526C---------------------------------------------------------------------- 2527C 2528#include "implicit.h" 2529#include "priunit.h" 2530#include "maxorb.h" 2531#include "ccorb.h" 2532#include "iratdef.h" 2533#include "cclr.h" 2534#include "ccsdsym.h" 2535#include "ccsdio.h" 2536#include "ccsdinp.h" 2537#include "r12int.h" 2538C 2539 LOGICAL LCONT,LOCDBG, TRIPLET 2540 CHARACTER*8 FC1AM,FC2AM,FRHO1,FRHO2,FR2SD,FRHO12,FC12AM,FS12AM, 2541 & FS2AM 2542 PARAMETER (FC1AM ='CCR_C1AM',FC2AM ='CCR_C2AM',FC12AM='CCR_C12M') 2543 PARAMETER (FRHO1 ='CCR_RHO1',FRHO2 ='CCR_RHO2',FRHO12='CCR_RH12') 2544 PARAMETER ( FS2AM ='CCR_S2DM',FS12AM='CCR_S12M') 2545 PARAMETER (TWO = 2.0D00, ZERO = 0.0D0) 2546 PARAMETER (LOCDBG = .FALSE.) 2547 CHARACTER*3 APROXR12 2548C 2549 DIMENSION WORK(LWORK), CONT(28) 2550C 2551 CALL QENTER('CC_ATRR') 2552 IF (CCS) THEN 2553 NCCVAR = NT1AM(ISYMV) 2554 ELSE 2555 NCCVAR = NT1AM(ISYMV) + NT2AM(ISYMV) 2556 IF (TRIPLET) NCCVAR = NCCVAR + NT2AM(ISYMV) 2557 IF (CCR12) NCCVAR = NCCVAR + NTR12AM(ISYMV) 2558 END IF 2559C 2560 NSIMTR = 1 2561 ISYMTR = ISYMV 2562C 2563 LUFR1 = -1 2564 LUFR2 = -1 2565 LUFR12 = -1 2566 LUFC1 = -1 2567 LUFC2 = -1 2568 LUFC12 = -1 2569 LUFS12 = -1 2570 LUFS2 = -1 2571C 2572C---------------- 2573C Open files. 2574C---------------- 2575C 2576 CALL CC_FILOP(FRHO1,LUFR1,FRHO2,LUFR2,FRHO12,LUFR12, 2577 * FC1AM,LUFC1,FC2AM,LUFC2,FC12AM,LUFC12, 2578 * FS12AM,LUFS12,FS2AM,LUFS2) 2579C 2580C---------------------------------------------------------------------- 2581C Make rho2 file name. 2582C For CCSD rho2 has to be stored on different file 2583C due to different length. 2584C---------------------------------------------------------------------- 2585C 2586 IF (.NOT. (CCS.OR.CC2)) THEN 2587 LUFSD = -1 2588 FR2SD = 'CC_TRA2_' 2589 ELSE 2590 LUFSD = LUFR2 2591 FR2SD = FRHO2 2592 ENDIF 2593C 2594C----------------------------------------------------------------------- 2595C Cheat and do CCS left transformation by right hand transformation. 2596C----------------------------------------------------------------------- 2597C 2598 MYSIDE = ISIDE 2599 IF ((ISIDE .EQ. -1 ) .AND. CCS ) MYSIDE = 1 2600 2601C 2602C----------------------------------------------------------------------- 2603C 2604 K1 = 1 2605 IVEC = K1 2606 IF (.NOT. (CCS .OR. CC2)) THEN 2607 ITR = 1 2608 ELSE 2609 ITR = K1 2610 ENDIF 2611C 2612C 2613C-------------------- 2614C Allocations. 2615C-------------------- 2616C 2617 NRHO2 = MAX(NT2AM(ISYMTR),NT2AM(1),2*NT2ORT(ISYMTR)) 2618CCH IF (ISIDE .EQ. -1) NRHO2 = MAX(NRHO2,2*NT2ORT(1)) 2619 IF (MYSIDE .EQ. -1) NRHO2 = MAX(NRHO2,2*NT2ORT(1)) 2620 2621 IF ( CC2 ) NRHO2 = MAX(NT2AM(ISYMTR),NT2AM(1)) 2622 IF ( CCS ) NRHO2 = 2 2623C 2624 IF (TRIPLET) NRHO2 = MAX(NRHO2,2*NT2AM(ISYMTR),NT2SQ(ISYMTR), 2625 & NT2ORT(ISYMTR)+NT2ORT3(ISYMTR)) 2626 NC2AM = MAX(NT2SQ(ISYMTR),NT2SQ(1), 2627 * NT2AM(ISYMTR)+2*NT2ORT(1),NT2R12(1)) 2628 IF ( CC2 ) NC2AM = MAX(NT2SQ(ISYMTR),NT2SQ(1)) 2629 IF ( CCS ) NC2AM = 2 2630C 2631 NRHO1 = NT1AM(ISYMTR)*NSIMTR 2632C 2633 KRHO1 = 1 2634 KRHO2 = KRHO1 + NRHO1 2635 KC1AM = KRHO2 + NRHO2 2636 KC2AM = KC1AM + NRHO1 2637 KEND1 = KC2AM + NC2AM 2638 LWRK1 = LWORK - KEND1 2639 IF ( LWRK1 .LE. 0 ) CALL QUIT('Insufficient workspace in CC_ATRR') 2640C 2641C 2642C--------------------------------------------------------------------- 2643C Prepare the C-amplitudes. 2644C--------------------------------------------------------------------- 2645C 2646 IF ( DEBUG .OR. LOCDBG ) THEN 2647 RHO1N = DDOT(NT1AM(ISYMTR),WORK(KRHO1),1,WORK(KRHO1),1) 2648 WRITE(LUPRI,1) 'Norm of C1AM -first in CC_ATRR: ',RHO1N 2649 IF (.NOT. CCS) THEN 2650 RHO2N = DDOT(NT2AM(ISYMTR),WORK(KRHO2),1,WORK(KRHO2),1) 2651 WRITE(LUPRI,1) 'Norm of C2AM -first in CC_ATRR: ',RHO2N 2652 ENDIF 2653 ENDIF 2654C 2655 IF (.NOT. CCS) THEN 2656 NRHO2 = NT2AM(ISYMTR) 2657 IF (TRIPLET) NRHO2 = 2*NRHO2 2658 END IF 2659C 2660 DO 70 IV = 1, NSIMTR ! Note that NSIMTR is always one! 2661 NR1 = IV + K1 - 1 2662 IF (.NOT. (CCS.OR.CC2)) THEN 2663 NR2 = IV 2664 ELSE 2665 NR2 = NR1 2666 ENDIF 2667 CALL CC_WVEC(LUFC1,FC1AM,NT1AM(ISYMTR),NT1AM(ISYMTR), 2668 * NR1,WORK(KRHO1)) 2669 IF (.NOT.CCS) THEN 2670 CALL CC_WVEC(LUFC2,FC2AM,NRHO2,NRHO2,NR2,WORK(KRHO2)) 2671 ENDIF 2672 IF (LCONT) THEN 2673 CONT(1) = DDOT(NT1AM(ISYMTR),WORK(KRHO1),1,WORK(KRHO1),1) 2674 IF (CCS) THEN 2675 CONT(2) = 0.0D0 2676 ELSE IF (CCR12.AND.IANR12.EQ.2) THEN 2677 CONT(2) = 0.0D0 2678 ELSE 2679 LRHO2 = NT2AM(ISYMTR) 2680 IF (TRIPLET) LRHO2 = 2*LRHO2 2681 CONT(2) = DDOT(LRHO2,WORK(KRHO2),1,WORK(KRHO2),1) 2682 END IF 2683 CONT(3) = 0.0D0 2684 END IF 2685 IF (CCR12) THEN 2686 KRHO12 = NT1AM(ISYMTR)+NRHO2+1 2687 CALL CC_WVEC(LUFC12,FC12AM,NTR12AM(ISYMTR),NTR12AM(ISYMTR), 2688 * NR1,WORK(KRHO12)) 2689 END IF 2690 70 CONTINUE 2691C 2692 CALL DCOPY(NT1AM(ISYMTR),WORK(KRHO1),1,WORK(KC1AM),1) 2693 IF (.NOT.TRIPLET) THEN 2694 IF ( .NOT. CCS ) THEN 2695 IF ( MYSIDE .GE. 1) THEN 2696 CALL CCLR_DIASCL(WORK(KRHO2),TWO,ISYMTR) 2697 ENDIF 2698 CALL CC_T2SQ(WORK(KRHO2),WORK(KC2AM),ISYMTR) 2699 ENDIF 2700 ENDIF 2701C 2702C--------------- 2703C File open. 2704C--------------- 2705C 2706C 2707 IF (.NOT. (CCS.OR.CC2)) THEN 2708 CALL WOPEN2(LUFSD,FR2SD,64,0) 2709 ENDIF 2710C 2711C--------------------- 2712C Zero rho vector. 2713C--------------------- 2714C 2715 NRHO2 = MAX(NT2AM(ISYMTR),2*NT2ORT(ISYMTR)) 2716 IF (TRIPLET.AND.ISIZE.EQ.-1) 2717 & NRHO2 = MAX(NT2SQ(ISYMTR),NT2ORT(ISYMTR)+NT2ORT3(ISYMTR)) 2718 IF ( CC2 ) NRHO2 = NT2AM(ISYMTR) 2719 IF ( CCS ) NRHO2 = 2 2720 CALL DZERO(WORK(KRHO1),NRHO1) 2721 CALL DZERO(WORK(KRHO2),NRHO2) 2722 NRHO2 = NT2AM(ISYMTR) 2723 IF(TRIPLET) NRHO2 = 2*NRHO2 2724 DO 80 IV = 1, NSIMTR 2725 NR1 = IV + K1 - 1 2726 IF (.NOT. (CCS.OR.CC2)) THEN 2727 NR2 = IV 2728 ELSE 2729 NR2 = NR1 2730 ENDIF 2731 CALL CC_WVEC(LUFR1,FRHO1,NT1AM(ISYMTR),NT1AM(ISYMTR), 2732 * NR1,WORK(KRHO1)) 2733 IF (.NOT.CCS) THEN 2734 CALL CC_WVEC(LUFSD,FR2SD,NRHO2,NRHO2,NR2,WORK(KRHO2)) 2735 ENDIF 2736 80 CONTINUE 2737C 2738C---------------------------------- 2739C Calculate transformed vectors. 2740C----------------------------------- 2741C 2742 LRHO1 = NT1AM(ISYMTR) 2743C 2744 IF (MYSIDE .EQ. 1) THEN 2745c get A*R 2746 IF (.NOT.TRIPLET) THEN 2747 CALL CC_RHTR(ECURR,FRHO1,LUFR1,FR2SD,LUFSD,FRHO12,LUFR12, 2748 * FC1AM,LUFC1,FC2AM,LUFC2,FC12AM,LUFC12, 2749 * WORK(KRHO1),WORK(KRHO2), 2750 * WORK(KC1AM),WORK(KC2AM), 2751 * WORK(KEND1),LWRK1,NSIMTR, 2752 * IVEC,ITR,LRHO1,LCONT,CONT(7),APROXR12) 2753 ELSE 2754 CALL CC_RHTR3(ECURR,FRHO1,LUFR1,FR2SD,LUFSD, 2755 & FC1AM,LUFC1,FC2AM,LUFC2, 2756 & WORK,LWORK,NSIMTR, 2757 & 1,1) 2758 END IF 2759 ELSE IF (MYSIDE .EQ. -1) THEN 2760c get R*A 2761 IF (.NOT.TRIPLET) THEN 2762 CALL CC_LHTR(ECURR, 2763 * FRHO1,LUFR1,FR2SD,LUFSD,FRHO12,LUFR12, 2764 * FC1AM,LUFC1,FC2AM,LUFC2,FC12AM,LUFC12, 2765 * WORK(KRHO1),WORK(KRHO2), 2766 * WORK(KC1AM),WORK(KC2AM), 2767 * WORK(KEND1),LWRK1,NSIMTR, 2768 * IVEC,ITR,LRHO1,APROXR12) 2769 ELSE 2770 CALL CC_LHTR3(ECURR, 2771 * FRHO1,LUFR1,FR2SD,LUFSD, 2772 * FC1AM,LUFC1,FC2AM,LUFC2, 2773 * WORK(KRHO1),WORK(KRHO2), 2774 * WORK(KC1AM),WORK(KC2AM), 2775 * WORK(KEND1),LWRK1,NSIMTR, 2776 * 1,1,LRHO1) 2777 END IF 2778 ELSE 2779 CALL QUIT('CC_ATRR; ISIDE should be -1 or +1 ') 2780 ENDIF 2781 2782 IF ( CCR12 ) THEN 2783 IF (MYSIDE .EQ. 1) THEN 2784 CALL CC_R12METRIC(ISYMTR,BRASCL,KETSCL,WORK(KEND1),LWRK1, 2785 * FC2AM,LUFC2,FC12AM,LUFC12,FS12AM,LUFS12, 2786 * FS2AM,LUFS2,IVEC,.FALSE.,DUMMY) 2787 ELSE IF (MYSIDE .EQ. -1) THEN 2788 CALL CC_R12METRIC(ISYMTR,0.5D0*KETSCL,2.0D0*BRASCL, 2789 * WORK(KEND1),LWRK1,FC2AM,LUFC2, 2790 * FC12AM,LUFC12,FS12AM,LUFS12,FS2AM,LUFS2, 2791 * IVEC,.FALSE.,DUMMY) 2792 END IF 2793 2794 KS12AM = NT1AM(ISYMTR)+NT2AM(ISYMTR)+NTR12AM(ISYMTR)+1 2795 CALL CC_RVEC(LUFS12,FS12AM,NTR12AM(ISYMTR),NTR12AM(ISYMTR), 2796 * IVEC,WORK(KS12AM)) 2797 IF (IANR12.EQ.2) THEN 2798 KS2AM = KS12AM + NTR12AM(ISYMTR) 2799 KEND1 = KS2AM + NT2AM(ISYMTR) 2800 LWRK1 = LWORK - KEND1 2801 IF (LWRK1.LT.0) THEN 2802 CALL QUIT('Insufficient work space in cc_atrr') 2803 END IF 2804 CALL CC_RVEC(LUFS2,FS2AM,NT2AM(ISYMTR),NT2AM(ISYMTR), 2805 * IVEC,WORK(KS2AM)) 2806 END IF 2807 END IF 2808 2809 NRHO2 = MAX(NT2AM(ISYMTR),2*NT2ORT(ISYMTR)) 2810 NRHO22 = NT2AM(ISYMTR) 2811 IF (CC2 ) NRHO2 = NT2AM(ISYMTR) 2812 IF (TRIPLET) THEN 2813 NRHO2 = NT2AM(ISYMTR) + NT2AMA(ISYMTR) 2814 NRHO22 = NT2AM(ISYMTR) + NT2AMA(ISYMTR) 2815 ENDIF 2816 DO 90 IV = 1, NSIMTR 2817 NR1 = IV + K1 - 1 2818 IF (.NOT. (CCS.OR.CC2)) THEN 2819 NR2 = IV 2820 ELSE 2821 NR2 = NR1 2822 ENDIF 2823 CALL CC_RVEC(LUFR1,FRHO1,NT1AM(ISYMTR),NT1AM(ISYMTR), 2824 * NR1,WORK(KRHO1)) 2825 IF (.NOT.CCS) THEN 2826 CALL CC_RVEC(LUFSD,FR2SD,NRHO2,NRHO22, 2827 * NR2,WORK(KRHO2)) 2828 END IF 2829 IF (CCR12) THEN 2830 KRHO12 = NT1AM(ISYMTR)+NT2AM(ISYMTR)+1 2831 CALL CC_RVEC(LUFR12,FRHO12,NTR12AM(ISYMTR),NTR12AM(ISYMTR), 2832 * NR1,WORK(KRHO12)) 2833 END IF 2834 2835 IF (LCONT) THEN 2836 IF (CCS) THEN 2837 KSCR12 = NT1AM(ISYMTR)+NT2AM(ISYMTR)+1 2838 KEND2 = KSCR12 + NT1AM(ISYMTR) 2839 ELSE IF (CCR12) THEN 2840 IF (IANR12.EQ.2) THEN 2841 KSCR12 = KS2AM + NT2AM(ISYMTR) 2842 ELSE 2843 KSCR12 = KS12AM + NTR12AM(ISYMTR) 2844 END IF 2845 KEND2 = KSCR12 + 2846 & MAX(NT1AM(ISYMTR),NT2AM(ISYMTR),NTR12AM(ISYMTR)) 2847 ELSE 2848 KSCR12 = NT1AM(ISYMTR)+NT2AM(ISYMTR)+1 2849 KEND2 = KSCR12 + MAX(NT1AM(ISYMTR),NT2AM(ISYMTR)) 2850 END IF 2851 LWRK2 = LWORK - KEND2 2852 IF(LWRK2.LE.0)CALL QUIT('Insufficient workspace in CC_ATRR') 2853 2854 CALL CC_RVEC(LUFC1,FC1AM,NT1AM(ISYMTR),NT1AM(ISYMTR), 2855 * 1,WORK(KSCR12)) 2856 CONT(4) = DDOT(NT1AM(ISYMTR),WORK(KRHO1),1,WORK(KSCR12),1) 2857 2858 IF (CCS) THEN 2859 CONT(5) = 0.0d0 2860 ELSE 2861 CALL CC_RVEC(LUFC2,FC2AM,NT2AM(ISYMTR),NT2AM(ISYMTR), 2862 & 1,WORK(KSCR12)) 2863 CONT(5) = DDOT(NT2AM(ISYMTR),WORK(KRHO2),1,WORK(KSCR12),1) 2864chf 2865 IF (CCR12.AND.IANR12.EQ.2) THEN 2866 CONT(2)= CONT(2)+ 2867 & DDOT(NT2AM(ISYMTR),WORK(KS2AM),1,WORK(KSCR12),1) 2868 END IF 2869chf 2870 END IF 2871 2872 IF (CCR12) THEN 2873 CALL CC_RVEC(LUFC12,FC12AM,NTR12AM(ISYMTR),NTR12AM(ISYMTR), 2874 * NR1,WORK(KSCR12)) 2875 CONT(6)=DDOT(NTR12AM(ISYMTR),WORK(KSCR12),1,WORK(KRHO12),1) 2876 CONT(3)=DDOT(NTR12AM(ISYMTR),WORK(KSCR12),1,WORK(KS12AM),1) 2877 ELSE 2878 CONT(6) = 0.0D0 2879 CONT(3) = 0.0D0 2880 END IF 2881 2882 END IF 2883 2884 IF ((IPRINT.GT.45) .OR. LOCDBG) THEN 2885 CALL AROUND('CC_ATRR: RHO = trans. Vector ') 2886 CALL CC_PRP(WORK(KRHO1),WORK(KRHO2),ISYMTR,1,1) 2887 IF (CCR12) THEN 2888 CALL CC_PRPR12(WORK(KRHO12),ISYMTR,1,.TRUE.) 2889 END IF 2890 ENDIF 2891 IF (.NOT.(CCS.OR.CC2)) THEN 2892 CALL CC_WVEC(LUFR2,FRHO2,NT2AM(ISYMTR), 2893 * NT2AM(ISYMTR),NR1,WORK(KRHO2)) 2894 ENDIF 2895 90 CONTINUE 2896C 2897C---------------------------- 2898C Close and delete files. 2899C---------------------------- 2900C 2901 CALL CC_FILCL(FRHO1,LUFR1,FRHO2,LUFR2,FRHO12,LUFR12, 2902 * FC1AM,LUFC1,FC2AM,LUFC2,FC12AM,LUFC12, 2903 * FS12AM,LUFS12,FS2AM,LUFS2) 2904 IF (.NOT. (CCS.OR.CC2)) THEN 2905 CALL WCLOSE2(LUFSD,FR2SD,'DELETE') 2906 ENDIF 2907C 2908 1 FORMAT(1x,A35,1X,E20.10) 2909C 2910 CALL QEXIT('CC_ATRR') 2911 RETURN 2912 END 2913c*DECK CC_REDEIG 2914 SUBROUTINE CC_REDEIG(WORK,LWORK,OMEHIT) 2915C 2916C----------------------------------------------------------------------------- 2917C 2918C Purpose: Find exci with OMEHIT and skip further calculation of 2919C as many of the other as possible. 2920C 2921C Written by Ove Christiansen 230899 2922C 2923C----------------------------------------------------------------------------- 2924 2925#include "implicit.h" 2926#include "priunit.h" 2927#include "cclr.h" 2928#include "ccorb.h" 2929#include "ccsdsym.h" 2930#include "ccsdinp.h" 2931#include "ccsections.h" 2932#include "ccexci.h" 2933#include "ccexgr.h" 2934#include "ccfdgeo.h" 2935#include "ccgr.h" 2936#include "cclres.h" 2937C 2938 DIMENSION WORK(LWORK) 2939C 2940 CALL QENTER('CC_REDEIG') 2941C 2942 WRITE (LUPRI,'(/,1X,A)') 2943 *'*********************************************************'// 2944 *'**********' 2945 WRITE(LUPRI,'(/,1X,A,F20.10)') 2946 * 'Search for reference excitation energy: ',OMEHIT 2947 WRITE(LUPRI,'(A)') ' List of excitation energies:' 2948 WRITE(LUPRI,'(A)') ' Nr. Sym. Multi Nr.in sym/mul Exci(au.)' 2949 DO IEXCI=1,NEXCI 2950 NR = IEXCI 2951 IS = ILSTSYM('RE ',IEXCI) 2952 IMUL = IMULTE(IEXCI) 2953 IF (IMUL.EQ.1) IST = IEXCI-ISYOFE(ILSTSYM('RE ',IEXCI)) 2954 IF (IMUL.EQ.3) IST = IEXCI-ITROFE(ILSTSYM('RE ',IEXCI)) 2955 OM = EIGVAL(IEXCI) 2956 WRITE(LUPRI,'(4I5,F20.10)') NR,IS,IMUL,IST,OM 2957 ENDDO 2958C 2959 OMEREF = 1.0D10 2960 IOMREF = 1 2961 DO IEXCI=1,NEXCI 2962 IF (ABS(OMEHIT-EIGVAL(IEXCI)).LT.OMEREF) THEN 2963 IOMREF = IEXCI 2964 OMEREF = ABS(OMEHIT-EIGVAL(IEXCI)) 2965 ENDIF 2966 ENDDO 2967C 2968 ISYHIT = ILSTSYM('RE ',IOMREF) 2969 IMUHIT = IMULTE(IOMREF) 2970 IF (IMUHIT.EQ.1) ISTHIT = IOMREF - ISYOFE(ISYHIT) 2971 IF (IMUHIT.EQ.3) ISTHIT = IOMREF - ITROFE(ISYHIT) 2972C 2973C We will, we will, MUH It!!! 2974C Let me hear it again 2975C We will, we will, MUH It!!! 2976C 2977 WRITE(LUPRI,'(/,F15.10,4(A,I3))') 2978 * OMEHIT,' is closest to exci. nr.',IOMREF, 2979 * ' which is nr. ',ISTHIT, 2980 * ' of symmetry ',ISYHIT, 2981 * ' and multiplicity ',IMUHIT 2982C 2983C Test if closest correspondence is ok. 2984C 2985 IF (MARGIN) THEN 2986 WRITE(LUPRI,'(2(A,F15.6))') ' Margin: ',XMARGIN, 2987 * ' correspondence:',OMEREF 2988 IF (OMEREF.LT.XMARGIN) THEN 2989 WRITE(LUPRI,'(/,A)') ' Closest correspondence is acceptable' 2990 ELSE 2991 WRITE(LUPRI,'(/,A)') 2992 * ' Closest correspondence is NOT acceptable' 2993 CALL QUIT( 2994 * ' Search for specific excitation energy was not satisfactory') 2995 ENDIF 2996 ENDIF 2997C 2998 DO IMULT = 1, 3, 2 2999 DO ISYM = 1, NSYM 3000 IF ((ISYM.EQ.ISYHIT).AND.(IMULT.EQ.IMUHIT)) THEN 3001 NCCEXCI(ISYM,IMULT) = ISTHIT 3002 ELSE 3003 NCCEXCI(ISYM,IMULT) = 0 3004 ENDIF 3005 ENDDO 3006 ENDDO 3007C 3008 DO I=1,ISTHIT 3009 IF (IMUHIT.EQ.1) EIGVAL(I) = EIGVAL(ISYOFE(ISYHIT)+I) 3010 IF (IMUHIT.EQ.3) EIGVAL(I) = EIGVAL(ITROFE(ISYHIT)+I) 3011 ENDDO 3012C ---------------------------- 3013C set up NEXCI + sym into et. 3014C ---------------------------- 3015 NEXCI = 0 3016 NTRIP = 0 3017 DO ISYM = 1,NSYM 3018 ISYOFE(ISYM) = NEXCI 3019 ITROFE(ISYM) = ISYOFE(ISYM) + NCCEXCI(ISYM,1) 3020 NEXCI = ITROFE(ISYM) + NCCEXCI(ISYM,3) 3021 NTRIP = NTRIP + NCCEXCI(ISYM,3) 3022 DO IEX = ISYOFE(ISYM)+1, NEXCI 3023 ISYEXC(IEX) = ISYM 3024 END DO 3025 DO IEX = ISYOFE(ISYM)+1, ITROFE(ISYM) 3026 IMULTE(IEX) = 1 3027 END DO 3028 DO IEX = ITROFE(ISYM)+1, NEXCI 3029 IMULTE(IEX) = 3 3030 END DO 3031 END DO 3032C 3033 CALL FLSHFO(LUPRI) 3034C 3035 IF (IPRINT.GT.15) THEN 3036 WRITE(LUPRI,*) 'IN CC_REDEIG after Reinit' 3037 WRITE(LUPRI,*) 'NEXCI: ',NEXCI 3038 WRITE(LUPRI,*) 'Singlet: ',(NCCEXCI(J,1),J=1,NSYM) 3039 WRITE(LUPRI,*) 'Triplet: ',(NCCEXCI(J,3),J=1,NSYM) 3040 WRITE(LUPRI,*) 'ISYOFE:',(ISYOFE(J), J=1,NSYM) 3041 WRITE(LUPRI,*) 'ITROFE:',(ISYOFE(J), J=1,NSYM) 3042 WRITE(LUPRI,*) 'ISYEXC:',(ISYEXC(J), J=1,NEXCI) 3043 WRITE(LUPRI,*) 'IMULTE:',(IMULTE(J), J=1,NEXCI) 3044 WRITE(LUPRI,*) 'EIGVAL:',(EIGVAL(J), J=1,NEXCI) 3045 ENDIF 3046C 3047C------------------------------- 3048C Overwriting IXSTAT IXSTSY. 3049C------------------------------- 3050C 3051 IXSTSY = ISYHIT 3052 IXSTAT = ISTHIT 3053 IXSTMU = IMUHIT 3054 OMECCX = EIGVAL(ISTHIT) 3055 ECCXST = ECCGRS + OMECCX 3056C 3057C----------------------------------------------------- 3058C Calculate only residues for this specific state. 3059C Input overwritten. 3060C----------------------------------------------------- 3061C 3062 IF (CCLRSD) THEN 3063 SELLRS =.TRUE. 3064 NSELRS = 1 3065 ISELRS(NSELRS,1) = ISYHIT 3066 ISELRS(NSELRS,2) = ISTHIT 3067 ENDIF 3068C 3069 WRITE(LUPRI,*) ' Note: Optimization and residue calculation '// 3070 * ' only carried out for this state' 3071C 3072 WRITE (LUPRI,'(/,1X,A)') 3073 *'*********************************************************'// 3074 *'**********' 3075C 3076 CALL QEXIT('CC_REDEIG') 3077C 3078 END 3079*=====================================================================* 3080 subroutine cctstsol(solvec,eigval,thrld,nred,nvec) 3081*---------------------------------------------------------------------* 3082* Purpose: test for degenerate eigenvalues/vectors. 3083* for degenerate solutions orthogonalize eigenvectors 3084* all eigenvectors are normalized to one 3085* 3086* Written by Christof Haettig, Februar 2003, Aarhus 3087*---------------------------------------------------------------------* 3088 implicit none 3089#include "priunit.h" 3090#include "cclr.h" 3091 3092 logical locdbg 3093 parameter (locdbg = .false.) 3094 3095 integer ivec, jvec, ndeg, nvec, nred 3096#if defined (SYS_CRAY) 3097 real fac, ddot, thrld, solvec(maxred,*), eigval(*) 3098#else 3099 double precision fac, ddot, thrld, solvec(maxred,*), eigval(*) 3100#endif 3101 3102 if (locdbg) then 3103 write(lupri,*) 'Matrix with solution vectors in reduced space:' 3104 call output(solvec,1,maxred,1,maxred,maxred,maxred,1,lupri) 3105 end if 3106 3107 do ivec = 1, nvec 3108 ndeg = 1 3109 do jvec = 1, nvec 3110 if (ivec.ne.jvec. and. 3111 & abs(eigval(ivec)-eigval(jvec)).lt.thrld) then 3112 ndeg = ndeg+1 3113 fac = -ddot(nred,solvec(1,jvec),1,solvec(1,ivec),1) 3114 call daxpy(nred,fac,solvec(1,jvec),1,solvec(1,ivec),1) 3115 end if 3116 end do 3117 if (ndeg.gt.1) then 3118 write(lupri,'(/1x,a,i3,f20.10)') 3119 & 'degenerate eigenvector found:',ivec,eigval(ivec) 3120 write(lupri,'(1x,a,i2,1x,a,/)') 3121 & 'eigenvector was orthogonalized against ',ndeg-1, 3122 & 'degenerate vectors' 3123 end if 3124 fac = dsqrt(ddot(nred,solvec(1,ivec),1,solvec(1,ivec),1)) 3125 if (dabs(fac-1.0d0).gt.1.0d-14) then 3126 call dscal(nred,1.0d0/fac,solvec(1,ivec),1) 3127 write(lupri,'(/1x,a,i3,f16.8,a,g16.8)') 3128 & 'eigenvector for state',ivec,eigval(ivec), 3129 & 'has been renormalized... scaling factor:',fac 3130 end if 3131 end do 3132 3133 return 3134 end 3135*=====================================================================* 3136 SUBROUTINE CC_EXCIT_CONT(ISTATE,ISYM,TRIPLET,MODEL,CONT,SCONT) 3137 IMPLICIT NONE 3138#include "priunit.h" 3139#include "ccsdinp.h" 3140#include "r12int.h" 3141 3142 CHARACTER*(*) MODEL 3143 CHARACTER*(8) MULTIPLICITY 3144 LOGICAL TRIPLET 3145 INTEGER ISTATE,ISYM 3146 3147#if defined (SYS_CRAY) 3148 REAL CONT(*), SCONT(*), EAE, ESE 3149#else 3150 DOUBLE PRECISION CONT(*), SCONT(*), EAE, ESE 3151#endif 3152 3153 MULTIPLICITY = 'singlet ' 3154 IF (TRIPLET) MULTIPLICITY = 'triplet ' 3155 3156 WRITE(LUPRI,'(//A)') 'Analysis of excitation energy:' 3157 WRITE(LUPRI,'(30("="),/)') 3158 WRITE(LUPRI,'(A,I3,5X,2A,5X,A,I3)')'Symmetry:',ISYM, 3159 * 'Multiplicity: ',MULTIPLICITY,'State nr.:',ISTATE 3160 WRITE(LUPRI,'(/72("-")/)') 3161 WRITE(LUPRI,'(3(A,F10.6,5X))') 'C1 * R1 :',CONT(1), 3162 * 'C1 * C1 :',SCONT(1),'Ratio:',CONT(1)/SCONT(1) 3163 IF (.NOT.CCS) THEN 3164 IF (CCR12.AND.(IANR12.EQ.1.OR.IANR12.EQ.3)) THEN 3165 WRITE(LUPRI,'(3(A,F10.6,5X))') 'C2 * R2 :',CONT(2), 3166 * 'C2 * C2 :',SCONT(2),'Ratio:', 3167 * CONT(2)/SCONT(2) 3168 ELSE IF (CCR12.AND.IANR12.EQ.2) THEN 3169 WRITE(LUPRI,'(3(A,F10.6,5X))') 'C2 * R2 :',CONT(2), 3170 * 'C2*C2 + C2*S23*C12 :',SCONT(2),'Ratio:', 3171 * CONT(2)/SCONT(2) 3172 ELSE 3173 WRITE(LUPRI,'(3(A,F10.6,5X))') 'C2 * R2 :',CONT(2), 3174 * 'C2 * C2 :',SCONT(2),'Ratio:', 3175 * CONT(2)/SCONT(2) 3176 END IF 3177 END IF 3178 IF (CCR12) THEN 3179 IF (IANR12.EQ.1.OR.IANR12.EQ.3) THEN 3180 WRITE(LUPRI,'(3(A,F10.6,5X))') 'C12 * R12:',CONT(3), 3181 * 'C12 * S * C12 :',SCONT(3),'Ratio:',CONT(3)/SCONT(3) 3182 ELSE IF (IANR12.EQ.2) THEN 3183 WRITE(LUPRI,'(3(A,F10.6,5X))') 'C12 * R12:',CONT(3), 3184 * 'C12*S32*C2 + C12*S12*C12:',SCONT(3),'Ratio:', 3185 * CONT(3)/SCONT(3) 3186 END IF 3187 END IF 3188 EAE = CONT(1)+CONT(2)+CONT(3) 3189 ESE = SCONT(1)+SCONT(2)+SCONT(3) 3190 WRITE(LUPRI,'(72("."))') 3191 WRITE(LUPRI,'(3(A,F10.6,5X),/)') 'C * R :', EAE, 3192 * ' C * C :',ESE, 'Ratio:',EAE/ESE 3193 3194 WRITE(LUPRI,'(A)') 'singles only contributions:' 3195 WRITE(LUPRI,'(13X,A,F10.6)') "J(T1) :",CONT(4)/ESE 3196 WRITE(LUPRI,'(13X,A,F10.6)') "J(E1) :",CONT(5)/ESE 3197 WRITE(LUPRI,'(A)') "...................................." 3198 WRITE(LUPRI,'(A,2F10.6,/)') "<R1|[Hhat,R1]|HF> :", 3199 * (CONT(4)+CONT(5))/ESE, (CONT(4)+CONT(5))/ESE 3200 3201 IF (.NOT.CCS) THEN 3202 WRITE(LUPRI,'(A)') 'doubles contributions:' 3203 WRITE(LUPRI,'(13X,A,F10.6)') "G(T2)+H(T2) :",CONT(7)/ESE 3204 WRITE(LUPRI,'(13X,A,F10.6)') "I(T2,E1) :",CONT(9)/ESE 3205 WRITE(LUPRI,'(A)') "...................................." 3206 WRITE(LUPRI,'(A,2F10.6)') "<R1|[[H,T2],R1]|HF> :", 3207 * (CONT(7)+CONT(9))/ESE, (CONT(7)+CONT(9))/ESE 3208 WRITE(LUPRI,'(13X,A,F10.6)') "G(E2)+H(E2) :",CONT(6)/ESE 3209 WRITE(LUPRI,'(13X,A,F10.6)') "I(E2,T1) :",CONT(8)/ESE 3210 WRITE(LUPRI,'(A)') "...................................." 3211 WRITE(LUPRI,'(A,F10.6)') "<R1|[Hhat,R2]|HF> :", 3212 * (CONT(6)+CONT(8))/ESE 3213 3214 IF (CC2) THEN 3215 WRITE(LUPRI,'(3X,A,F10.6)') "<R2|[Hhat,R1]|HF> :", 3216 * CONT(10)/ESE 3217 WRITE(LUPRI,'(3X,A,F10.6)') "<R2|[F,R2]|HF> :", 3218 * CONT(11)/ESE 3219 WRITE(LUPRI,'(A)') "...................................." 3220 WRITE(LUPRI,'(A,2F10.6)') "sum of doubles blocks :", 3221 * (CONT(6)+CONT(8)+CONT(10)+CONT(11))/ESE, 3222 * (CONT(6)+CONT(8)+CONT(10)+CONT(11))/ESE 3223 WRITE(LUPRI,'(A,F10.6)') "doubles total :", 3224 * (CONT(7)+CONT(9)+CONT(6)+CONT(8)+CONT(10)+CONT(11))/ESE 3225 END IF 3226 END IF 3227 3228 IF (CCSD) THEN 3229 WRITE(LUPRI,'(A)') 'Analysis incomplete for CCSD!' 3230 END IF 3231 3232 IF (CCR12) THEN 3233 WRITE(LUPRI,'(/A)') 'R12 contributions:' 3234 WRITE(LUPRI,'(4X,A,2F10.6)') "<R1|[[H,T2'],R1]|HF> :", 3235 * CONT(18)/ESE, CONT(18)/ESE 3236 IF (IANR12.EQ.2.OR.IANR12.EQ.3) THEN 3237 WRITE(LUPRI,'(4X,A,2F10.6)') "1HP and 1IP contrib :", 3238 * CONT(22)/ESE 3239 END IF 3240 WRITE(LUPRI,'(4X,A,F10.6)') "<R1|[Hhat,R2']|HF> :", 3241 * CONT(19)/ESE 3242 IF (IANR12.EQ.2.OR.IANR12.EQ.3) THEN 3243 WRITE(LUPRI,'(4X,A,F10.6)') "2HP and 2IP contrib :", 3244 * CONT(23)/ESE 3245 END IF 3246 WRITE(LUPRI,'(4X,A,F10.6)') "<R2'|[Hhat,R1]|HF> :", 3247 * CONT(20)/ESE 3248 WRITE(LUPRI,'(4X,A,F10.6)') "<R2'|[F,R2']|HF> :", 3249 * CONT(21)/ESE 3250 IF (IANR12.EQ.2.OR.IANR12.EQ.3) THEN 3251 WRITE(LUPRI,'(4X,A,F10.6)') "<R2'|[F,R2]|HF> :", 3252 * CONT(12)/ESE 3253 WRITE(LUPRI,'(4X,A,F10.6)') "<R2 |[F,R2']|HF> :", 3254 * CONT(13)/ESE 3255 WRITE(LUPRI,'(A)') "...................................." 3256 WRITE(LUPRI,'(A,2F10.6)') "sum of R12 doubles blocks:", 3257 * (CONT(19)+CONT(20)+CONT(21)+CONT(12)+CONT(13))/ESE, 3258 * (CONT(19)+CONT(20)+CONT(21)+CONT(12)+CONT(13))/ESE 3259 WRITE(LUPRI,'(A,F10.6)') "R12 doubles total :", 3260 * (CONT(18)+CONT(19)+CONT(20)+CONT(21)+CONT(12)+CONT(13))/ESE 3261 WRITE(LUPRI,'(A)') "...................................." 3262 WRITE(LUPRI,'(A,10X,F10.6)') "total excitation energy :", 3263 * (CONT(4)+CONT(5)+CONT(6)+CONT(7)+CONT(8)+CONT(9)+ 3264 * CONT(10)+CONT(11)+CONT(18)+CONT(19)+CONT(20)+ 3265 * CONT(21)+CONT(12)+CONT(13))/ESE 3266 ELSE IF (IANR12.EQ.1) THEN 3267 WRITE(LUPRI,'(A)') "...................................." 3268 WRITE(LUPRI,'(A,2F10.6)') "sum of R12 doubles blocks:", 3269 * (CONT(19)+CONT(20)+CONT(21))/ESE, 3270 * (CONT(19)+CONT(20)+CONT(21))/ESE 3271 3272 WRITE(LUPRI,'(A,F10.6)') "R12 doubles total :", 3273 * (CONT(18)+CONT(19)+CONT(20)+CONT(21))/ESE 3274 WRITE(LUPRI,'(A)') "...................................." 3275 WRITE(LUPRI,'(A,10X,F10.6)') "total excitation energy :", 3276 * (CONT(4)+CONT(5)+CONT(6)+CONT(7)+CONT(8)+CONT(9)+ 3277 * CONT(10)+CONT(11)+CONT(18)+CONT(19)+CONT(20)+ 3278 * CONT(21))/ESE 3279 END IF 3280 END IF 3281 3282 WRITE(LUPRI,'(72("-"))') 3283 WRITE(LUPRI,'(//)') 3284 3285 RETURN 3286 END 3287*=====================================================================* 3288