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_slv */ 20 SUBROUTINE CC_SLV(AODEN,ETLM,DIELCONV,WORK,LWORK) 21C 22C----------------------------------------------------------------------------- 23C 24C Purpose: Direct calculation of Coupled Cluster solvent effects. 25C 26C CCS(CIS/HF)(nci), MP2(nci), CC2(nci), CCSD, CC3(nci), MCC2(nci) 27C 28C SLV98,OC 29C Ove Christiansen, April 1998. 30C 31C----------------------------------------------------------------------------- 32C 33#include "implicit.h" 34#include "maxorb.h" 35#include "mxcent.h" 36 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, IZERO = 0 , TWO = 2.0D0) 37#include "dummy.h" 38#include "iratdef.h" 39#include "priunit.h" 40#include "ccorb.h" 41#include "ccsdsym.h" 42#include "ccsdio.h" 43#include "ccsdinp.h" 44#include "ccfield.h" 45#include "exeinf.h" 46#include "ccfdgeo.h" 47#include "ccslvinf.h" 48#include "ccinftap.h" 49C 50 LOGICAL DIELCONV 51 DIMENSION WORK(LWORK),AODEN(*),ETLM(NLMCU,2) 52 CHARACTER MODEL*10 53 CHARACTER MODELPRI*4 54C 55C----------- 56C Header 57C----------- 58C 59 WRITE (LUPRI,'(1X,A,/)') ' ' 60 WRITE (LUPRI,'(1X,A)') 61 *'*********************************************************'// 62 *'**********' 63 WRITE (LUPRI,'(1X,A)') 64 *'* Output from coupled cluster solvent program '// 65 *' *' 66 WRITE (LUPRI,'(1X,A,/)') 67 *'*********************************************************'// 68 *'**********' 69C 70 WRITE(LUPRI,'(/,1X,A)') 71 *'+=====================================================+' 72 WRITE(LUPRI,'(1X,A)') 73 *'| Lmax Cavity(a.u.) Epsilon_st Epsilon_op |' 74 WRITE(LUPRI,'(1X,A)') 75 *'+-----------------------------------------------------+' 76 WRITE(LUPRI,'(1X,A,I3,4X,F13.8,2X,F13.8,3X,F13.8,A)') 77 * '| ',LMAXCU,RCAVCU,EPSTCU,EPOPCU,' |' 78 WRITE(LUPRI,'(1X,A)') 79 * '+=====================================================+' 80C 81 MODEL = 'CCSD' 82 IF (CC2.AND.(.NOT.MCC2)) THEN 83 CALL AROUND('Coupled Cluster model is: CC2') 84 MODEL = 'CC2' 85 MODELPRI = ' CC2' 86 ENDIF 87 IF (MCC2) THEN 88 CALL AROUND('Coupled Cluster model is: MCC2') 89 MODEL = 'MCC2' 90 MODELPRI = 'MCC2' 91 ENDIF 92 IF (MP2) THEN 93 CALL AROUND('Model is second order pert. theory: MP2 ') 94 MODEL = 'MP2' 95 MODELPRI = ' MP2' 96 ENDIF 97 IF (CCS.AND.(.NOT.CIS)) THEN 98 CALL AROUND('Coupled Cluster model is: CCS') 99 MODEL = 'CCS' 100 MODELPRI = ' CCS' 101 ENDIF 102 IF (CCS.AND.CIS) THEN 103 CALL AROUND('CIS model in use ') 104 MODEL = 'CCS' 105 MODELPRI = ' CIS' 106 ENDIF 107 IF (CCD) THEN 108 CALL AROUND('Coupled Cluster model is: CCD') 109 MODEL = 'CCD' 110 MODELPRI = ' CCD' 111 ENDIF 112 IF (CC3 ) THEN 113 CALL AROUND('Coupled Cluster model is: CC3') 114 MODEL = 'CC3' 115 MODELPRI = ' CC3' 116 CALL QUIT('CC3 first order properties not implemented') 117 ENDIF 118 IF (CC1A) THEN 119 CALL AROUND('Coupled Cluster model is: CCSDT-1a') 120 MODEL = 'CCSDT-1a' 121 CALL QUIT( 'CCSDT-1a first order properties not implemented') 122 ENDIF 123 IF (CC1B) THEN 124 CALL AROUND('Coupled Cluster model is: CCSDT-1b') 125 MODEL = 'CCSDT-1b' 126 CALL QUIT( 'CCSDT-1b first order properties not implemented') 127 ENDIF 128 IF (CCSD) THEN 129 CALL AROUND('Coupled Cluster model is: CCSD') 130 MODEL = 'CCSD' 131 MODELPRI = 'CCSD' 132 ENDIF 133C 134c IF ((.NOT. CCSD).OR.(.NOT.CC2)) THEN 135c CALL QUIT( ' Solvent not implemented for other than CC2 and CCSD') 136c ENDIF 137C 138 IF (IPRINT.GT.10) WRITE(LUPRI,*) 'CC_SLV-1: Workspace:',LWORK 139C 140C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 141C SLV98,OC 142C Solvent section 143C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 144C 145C 146C------------------------------------------------------ 147C Calculate the solvent contribution to the energy. 148C------------------------------------------------------ 149C 150 CALL CC_SLVE(ECCSLCON,AODEN,ETLM,WORK,LWORK) 151C 152C-------------------------------------------------------------- 153C Calculate new solvent energy. 154C-------------------------------------------------------------- 155C 156 ECCCU = ECCGRS + ECCSLCON 157 IF (ABS(ECCCU-ECCPR).LT.CVGESOL) LSLECVG = .TRUE. 158 WRITE(LUPRI,'(1X,A,I3,A,F20.10)') 159 * 'Solvent energy contribution in iteration',ICCSLIT,': ', 160 * ECCSLCON 161 WRITE(LUPRI,'(1X,A,F20.10)') 162 * 'CC energy in the current solvent iteration: ',ECCCU 163 WRITE(LUPRI,'(1X,A,F20.10)') 164 * 'CC energy in the previous solvent iteration: ',ECCPR 165 WRITE(LUPRI,*)'LSLECVG: ',LSLECVG 166 WRITE(LUPRI,*) 167 *' Change in Total energy in this solvent it.:', 168 * ECCCU - ECCPR 169 ECCPR = ECCCU 170 DIELCONV = .FALSE. 171 IF (LSLECVG.AND.LSLTCVG.AND.LSLLCVG) THEN 172 DIELCONV = .TRUE. 173 WRITE(LUPRI,'(/,1X,A,I3,A)') 174 *'Coupled cluster solvent equations are converged in ',ICCSLIT, 175 *' solvent iterations' 176 WRITE(LUPRI,'(/,1X,A8,A,F30.16,/)') 177 * MODEL,'converged energy in solvent: ',ECCCU 178 WRITE(LUPRI,'(/,1X,A8,A,F30.16,/)') 179 * MODEL,'solvation energy: ',ECCSLCON 180 WRITE(LURES,'(/,1X,A,I2,A,F8.4,A,F8.4,A,F8.4,A)') 181 * 'Solvent: L_max=',LMAXCU,', R_cav=',RCAVCU,', Eps_st =',EPSTCU, 182 * ', Eps_op =',EPOPCU,': ' 183 WRITE(LURES,'(12X,A8,A,F20.10)') 184 * MODEL,' Total energy: ',ECCCU 185 WRITE(LURES,'(12X,A8,A,F20.10,/)') 186 * MODEL,' Solvation energy: ',ECCSLCON 187 ECCGRS = ECCCU 188 ELSE 189 ICCSLIT = ICCSLIT + 1 190 IF (ICCSLIT.GT.MXCCSLIT) THEN 191 WRITE(LUPRI,*) 'Maximum number of solvent iterations', 192 * MXCCSLIT,' is reached' 193 CALL QUIT( 'Maximum number of solvent iterations reached') 194 ENDIF 195 ENDIF 196C 197C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 198C SLV98,OC 199C End of solvent section. 200C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 201C 202 WRITE (LUPRI,'(1X,A)') 203 *'*********************************************************'// 204 *'**********' 205 WRITE (LUPRI,'(1X,A)') 206 *'* End of coupled cluster solvent program '// 207 *' *' 208 WRITE (LUPRI,'(1X,A)') 209 *'*********************************************************'// 210 *'**********' 211C 212 END 213C /* Deck cc_slve */ 214 SUBROUTINE CC_SLVE(ESOLT,AODEN,ETLM,WORK,LWORK) 215C 216C----------------------------------------------------------------------------- 217C 218C Purpose: Calculation of Coupled Cluster solvent energy. 219C for given solvent defined by RCAVCU,EPSTCU,LMAXCU. 220C 221C Based on SOLGRD from Sirius. 222C 223C SLV98,OC 224C Ove Christiansen, April 1998. 225C 226C----------------------------------------------------------------------------- 227C 228#include "implicit.h" 229#include "maxorb.h" 230#include "mxcent.h" 231 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, IZERO = 0 , TWO = 2.0D0) 232#include "dummy.h" 233#include "iratdef.h" 234#include "priunit.h" 235#include "ccorb.h" 236#include "ccsdsym.h" 237#include "ccsdio.h" 238#include "ccsdinp.h" 239#include "ccfield.h" 240#include "exeinf.h" 241#include "ccfdgeo.h" 242#include "ccslvinf.h" 243#include "ccinftap.h" 244#include "thrzer.h" 245C 246 DIMENSION WORK(LWORK),AODEN(*),ETLM(NLMCU,2) 247 CHARACTER MODEL*10 248C 249 LOGICAL FIRST 250 SAVE FIRST 251 DATA FIRST /.TRUE./ 252C 253 DOUBLE PRECISION GL 254 EXTERNAL GL 255C 256C----------------------- 257C Check g_l factors. 258C----------------------- 259C 260 IF (IPRINT.GT.10) THEN 261 WRITE(LUPRI,*) 'L G_l(Eps) for Eps,Rcav=',EPSTCU,RCAVCU 262 DO L = 0, LMAXCU 263 WRITE(LUPRI,*) L,GL(L,EPSTCU,RCAVCU) 264 ENDDO 265 ENDIF 266C 267C------------------------------------------------------- 268C Check integrals and readin nuclear contributions. 269C------------------------------------------------------- 270C 271 LUCSOL = -1 272 CALL GPOPEN(LUCSOL,'AOSOLINT','OLD',' ','UNFORMATTED',IDUMMY, 273 & .FALSE.) 274 REWIND(LUCSOL) 275 CALL MOLLAB('SOLVRLM ',LUCSOL,LUERR) 276C 277 IF (FIRST) THEN 278 READ (LUCSOL) LMAXSS, LMTOT, NNNBAS 279 IF (LMAXSS .LT. LMAXCU) THEN 280 WRITE (LUERR,'(//2A,2(/A,I5))') ' --- CC_SLVE ERROR,', 281 * ' insufficient number of intgrals on LUCSOL', 282 * ' l max from CC input :',LMAXCU, 283 * ' l max from LUCSOL file :',LMAXSS 284 CALL QUIT('CC_SLVE: lmax on LUCSOL is too small') 285 END IF 286 IF ((LMAXSS+1)**2 .NE. LMTOT) THEN 287 WRITE (LUERR,'(//2A,3(/A,I5))') ' --- CC_SLVE ERROR,', 288 * ' LUCSOL file info inconsistent', 289 * ' l_max :',LMAXSS, 290 * ' (l_max + 1) ** 2 :',(LMAXSS+1)**2, 291 * ' LMTOT :',LMTOT 292 CALL QUIT('CC_SLVE: LUCSOL info not internally consistent') 293 END IF 294 IF (NNNBAS .NE. NBAST) THEN 295 WRITE (LUERR,'(//2A,3(/A,I5))') ' --- CC_SLVE ERROR,', 296 * ' LUCSOL file info inconsistent with SIRIUS input', 297 * ' NBAST - LUCSOL :',NNNBAS, 298 * ' NBAST - SIRIUS :',NBAST 299 CALL QUIT('CC_SLVE: LUCSOL info not '// 300 & 'consistent with SIRIUS input.') 301 END IF 302 FIRST = .FALSE. 303 ELSE 304 READ (LUCSOL) 305 END IF 306 CALL READT(LUCSOL,NLMCU,ETLM(1,2)) 307C 308C------------------------------------- 309C Print out nuclear contributions. 310C------------------------------------- 311C 312 IF (IPRINT .GE. 10) THEN 313 WRITE(LUPRI,'(/A/)') 314 * ' l, m, Tn(lm) - the nuclear contributions :' 315 LM = 0 316 DO 220 L = 0,LMAXCU 317 DO 210 M = -L,L 318 LM = LM + 1 319 WRITE(LUPRI,'(2I5,F15.10)') L,M,ETLM(LM,2) 320 210 CONTINUE 321 WRITE(LUPRI,'()') 322 220 CONTINUE 323 END IF 324C 325C----------------------------------------- 326C Calculate electronic conctributions. 327C 1. Loop L 328C 2. Get symmetries of Tlm 329C 3. Loop M 330C----------------------------------------- 331C 332 LM = 0 333 DO 520 L = 0,LMAXCU 334 READ (LUCSOL) L1,(ISYTLM(M),M=1,2*L+1) 335 IF (L1 .NE. L) THEN 336 WRITE (LUERR,*) 337 & 'ERROR CC_SLVE: L from LUCSOL not as expected' 338 WRITE (LUERR,*) 'L from 520 loop:',L 339 WRITE (LUERR,*) 'L from LUCSOL :',L1 340 CALL QUIT('ERROR CC_SLVE: L from LUCSOL not as expected') 341 END IF 342C 343 DO 500 M = -L,L 344 LM = LM + 1 345 IF (IPRINT .GE. 15) THEN 346 WRITE(LUPRI,'(/A,2I5/A)') ' --- l, m :',L,M, 347 * ' ====================' 348 WRITE(LUPRI,'(A,I2)') ' Symmetry :',ABS(ISYTLM(L+M+1)) 349 END IF 350 IF (ISYTLM(L+M+1) .NE. 1) THEN 351Chj IF (ABS(ETLM(LM,2)) .GT. THRZER) THEN 352Chj numerical errors have seen to give 1.0D-12, 353Chj so we now use 1.0D-10 for test /Feb 2005 354 IF (ABS(ETLM(LM,2)) .GT. 1.0D-10) THEN 355 WRITE(LUPRI,*) 'ERROR CC_SLVE for l,m',L,M 356 WRITE(LUPRI,*) 'Symmetry :',ISYTLM(L+M+1) 357 WRITE(LUPRI,*) 'Tn(l,m) .ne. 0, but =',ETLM(LM,2) 358 CALL QUIT( 'ERROR CC_SLVE: Tn(l,m) not 0 as expected') 359 END IF 360 ETLM(LM,2) = ZERO 361C ... to fix round-off errors in Tn(l,m) calculation 362 IF (ISYTLM(L+M+1) .GT. 1) READ (LUCSOL) 363 GO TO 500 364 END IF 365C 366C-------------------------------- 367C Read T(l,m) in AO basis. 368C-------------------------------- 369C 370 KTLMAO = 1 371 KWRK1 = KTLMAO + N2BST(ISYMOP) 372 LWRK1 = LWORK - KWRK1 373 IF (LWRK1.LT.0) CALL QUIT( 'Too little work in CC_SLVE, 1') 374C 375 CALL CC_SLVINT(WORK(KTLMAO),WORK(KWRK1),LWRK1,ISYMOP) 376C 377 IF (IPRINT .GT. 50) THEN 378 CALL AROUND('cc_slve: Tlm_ao matrix: - cc storage') 379 CALL CC_PRFCKAO(WORK(KTLMAO),ISYMOP) 380 ENDIF 381C 382 IF (IPRINT .GT. 50) THEN 383 CALL AROUND('One electron density in cc_slve') 384 CALL CC_PRFCKAO(AODEN,ISYMOP) 385 ENDIF 386C 387C------------------------------------------------------ 388C Add electronic contribution TE(l,m) to T(l,m) 389C Contract CC density with integrals. 390C SLV98,OC 391C------------------------------------------------------ 392C 393 TELM = DDOT(N2BST(ISYMOP),WORK(KTLMAO),1,AODEN,1) 394C 395 IF (IPRINT .GE. 6) THEN 396 WRITE(LUPRI,'(A,2I5,/A,3F17.8)') 397 * ' --- l, m :',L,M, 398 * ' Te(lm), Tn(lm), T(lm) :', 399 * TELM,ETLM(LM,2),ETLM(LM,2)-TELM 400 END IF 401C 402 ETLM(LM,2) = ETLM(LM,2) - TELM 403C 404C To avoid numerical stability problems. 405C SLV98,OC Necessary??? 406C 407 IF (ABS(ETLM(LM,2)) .LE. THRZER) THEN 408 ETLM(LM,2) = ZERO 409 GO TO 500 410 END IF 411C 412 500 CONTINUE 413 520 CONTINUE 414C 415C--------------------------------- 416C Add up energy contributions. 417C Save <Tlm>'s. 418C--------------------------------- 419C 420 ESOLT = ZERO 421 LM = 0 422 DO 920 L = 0,LMAXCU 423 DO 900 M = -L,L 424 LM = LM + 1 425 CCTLM(LM) = ETLM(LM,2) 426 ETLM(LM,1) = GL(L,EPSTCU,RCAVCU)*ETLM(LM,2)*ETLM(LM,2) 427 IF (IPRINT.GT.5) THEN 428 WRITE(LUPRI,*) 'L,M,Energy cont.',L,M,ETLM(LM,1) 429 WRITE(LUPRI,*) 'ETLM2,GL',ETLM(LM,2),GL(L,EPSTCU,RCAVCU) 430 ENDIF 431 ESOLT = ESOLT + ETLM(LM,1) 432 900 CONTINUE 433 920 CONTINUE 434C 435 CALL CCSL_TLMPUT 436C 437 CALL GPCLOSE(LUCSOL,'KEEP') 438C 439 END 440C 441 DOUBLE PRECISION FUNCTION GL(L,EPS,A) 442C 443C----------------------------------------------------------------------------- 444C 445C Purpose: Calculate G_l(eps) factors for solvent calculations. 446C Cavity radius = A and dielectric constant = EPS 447C based on solfl 448C 449C Ove Christiansen, April 1998 450C 451C----------------------------------------------------------------------------- 452C 453#include "implicit.h" 454C 455C Parameters: FLFAC = - 1 / (8 pi epsilon_null) (S.I.) 456C = - 1 / 2 (a.u.) 457C 458 PARAMETER ( FLFAC = -0.5D0 , D1 = 1.0D0 ) 459C 460 RL = L 461 GL = FLFAC * A**(-(2*L+1))* (RL + D1) * (EPS - D1) 462 * / (RL + EPS*(RL + D1)) 463 END 464C 465 DOUBLE PRECISION FUNCTION GL2(L,EPSST,EPSOP,A) 466C 467C----------------------------------------------------------------------------- 468C 469C Purpose: Calculate G_l(epsst,epsop) factors for solvent calculations. 470C Ove Christiansen, April 1998 471C 472C----------------------------------------------------------------------------- 473C 474#include "implicit.h" 475C 476 EXTERNAL GL 477 DOUBLE PRECISION GL 478C 479 GL2 = GL(L,EPSST,A) - GL(L,EPSOP,A) 480C 481 END 482 SUBROUTINE CC_AOUP(AOUP,AOP,WORK,LWORK,ISYM) 483C 484C----------------------------------------------------------------------------- 485C Purpose: 486C Transform from packed to symmetry packed but unpacked form. 487C 488C Ove Christiansen, April 1998. 489C 490C----------------------------------------------------------------------------- 491C 492#include "implicit.h" 493#include "maxorb.h" 494#include "mxcent.h" 495#include "dummy.h" 496#include "iratdef.h" 497#include "priunit.h" 498#include "ccorb.h" 499#include "ccsdinp.h" 500#include "ccsdsym.h" 501C 502 DIMENSION AOUP(*),AOP(*),WORK(LWORK) 503C 504 IF (LWORK.LT.(NBAST**2)) CALL QUIT( 'too little work in CC_AOUP') 505C 506 CALL DSPTGE(NBAST,AOP(1),WORK(1)) 507C 508 IF (IPRINT.GT.50) THEN 509 CALL AROUND('CC_AOUP: Integrals') 510 CALL OUTPUT(WORK(1),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI) 511 END IF 512C 513 DO ISYMA = 1, NSYM 514 ISYMB = MULD2H(ISYM,ISYMA) 515 DO A = 1, NBAS(ISYMA) 516 DO B = 1, NBAS(ISYMB) 517 IOLD = (IBAS(ISYMA)+A-1)*NBAST + (IBAS(ISYMB)+B) 518 INEW = IAODIS(ISYMB,ISYMA) + NBAS(ISYMB)*(A-1) + B 519 AOUP(INEW) = WORK(IOLD) 520 END DO 521 END DO 522 END DO 523 END 524C /* Deck cc_slvint */ 525 SUBROUTINE CC_SLVINT(TLMAO,WORK,LWORK,ISYMT) 526C 527C----------------------------------------------------------------------------- 528C 529C Purpose: Readin solvent integrals in coupled cluster format. 530C 531C SLV98,OC 532C Ove Christiansen, April 1998. 533C 534C----------------------------------------------------------------------------- 535C 536#include "implicit.h" 537#include "maxorb.h" 538#include "mxcent.h" 539 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, IZERO = 0 , TWO = 2.0D0) 540#include "dummy.h" 541#include "iratdef.h" 542#include "priunit.h" 543#include "ccorb.h" 544#include "ccsdsym.h" 545#include "ccsdio.h" 546#include "ccsdinp.h" 547#include "ccfield.h" 548#include "exeinf.h" 549#include "ccfdgeo.h" 550#include "ccslvinf.h" 551#include "ccinftap.h" 552C 553 DIMENSION WORK(LWORK),TLMAO(*) 554C 555 IF (IPRINT.GT.10) THEN 556 WRITE(LUPRI,*) 'CC_SLVINT: Read in integrals' 557 WRITE(LUPRI,*) 'Input symmetry claimed', isymt 558 ENDIF 559C 560 KTLMAOP = 1 561 KWRK1 = KTLMAOP + NNBASX 562 LWRK1 = LWORK - KWRK1 563 IF (LWRK1.LT.2*NNBASX) CALL QUIT( 'Too little work in CC_SLVINT') 564C 565 CALL READT(LUCSOL,NNBASX,WORK(KTLMAOP)) 566C 567 IF (IPRINT .GE. 25) THEN 568 CALL AROUND('CC_SLVINT: Tlm_ao matrix: - packed ') 569 CALL OUTPAK(WORK(KTLMAOP),NBAST,1,LUPRI) 570 END IF 571C 572 CALL CC_AOUP(TLMAO,WORK(KTLMAOP),WORK(KWRK1),LWRK1,ISYMT) 573C 574 END 575C /* Deck ccsl_rhstg */ 576 SUBROUTINE CCSL_RHSTG(FOCK,WORK,LWORK) 577C 578C----------------------------------------------------------------------------- 579C 580C Purpose: Direct calculation of Coupled Cluster 581C solvent effects. 582C 583C SLV98,OC 584C Ove Christiansen, April 1998. 585C 586C----------------------------------------------------------------------------- 587C 588#include "implicit.h" 589#include "maxorb.h" 590#include "mxcent.h" 591 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, IZERO = 0 , TWO = 2.0D0) 592#include "dummy.h" 593#include "iratdef.h" 594#include "priunit.h" 595#include "ccorb.h" 596#include "ccsdsym.h" 597#include "ccsdio.h" 598#include "ccsdinp.h" 599#include "ccfield.h" 600#include "exeinf.h" 601#include "ccfdgeo.h" 602#include "ccslvinf.h" 603#include "ccinftap.h" 604C 605 DIMENSION WORK(LWORK),FOCK(*) 606C 607 LOGICAL FIRST 608 SAVE FIRST 609 DATA FIRST/.TRUE./ 610C 611 IF (IPRINT.GT.10) THEN 612 WRITE(LUPRI,*) 613 * 'CCSL_RHSTG: Solvent contribution to CC equations.' 614 ENDIF 615C 616C--------------------------------- 617C Read in integrals from file. 618C--------------------------------- 619C 620 LUCSOL = -1 621 CALL GPOPEN(LUCSOL,'AOSOLINT','OLD',' ','UNFORMATTED',IDUMMY, 622 & .FALSE.) 623 REWIND(LUCSOL) 624 CALL MOLLAB('SOLVRLM ',LUCSOL,LUERR) 625 IF (FIRST) THEN 626 READ (LUCSOL) LMAXSS, LMTOT, NNNBAS 627 IF (LMAXSS .LT. LMAXCU) THEN 628 WRITE (LUERR,'(//2A,2(/A,I5))') ' --- CCSL_RHSTG ERROR,', 629 * ' insufficient number of intgrals on LUCSOL', 630 * ' l max from CC input :',LMAXCU, 631 * ' l max from LUCSOL file :',LMAXSS 632 CALL QUIT('CCSL_RHSTG: lmax on LUCSOL is too small') 633 END IF 634 IF ((LMAXSS+1)**2 .NE. LMTOT) THEN 635 WRITE (LUERR,'(//2A,3(/A,I5))') ' --- CCSL_RHSTG ERROR,', 636 * ' LUCSOL file info inconsistent', 637 * ' l_max :',LMAXSS, 638 * ' (l_max + 1) ** 2 :',(LMAXSS+1)**2, 639 * ' LMTOT :',LMTOT 640 CALL QUIT('CCSL_RHSTG: LUCSOL info not internally consistent') 641 END IF 642 IF (NNNBAS .NE. NBAST) THEN 643 WRITE (LUERR,'(//2A,3(/A,I5))') ' --- CCSL_RHSTG ERROR,', 644 * ' LUCSOL file info inconsistent with SIRIUS input', 645 * ' NBAST - LUCSOL :',NNNBAS, 646 * ' NBAST - SIRIUS :',NBAST 647 CALL QUIT('CCSL_RHSTG: LUCSOL info not '// 648 & 'consistent with SIRIUS input.') 649 END IF 650 FIRST = .FALSE. 651 ELSE 652 READ (LUCSOL) 653 END IF 654 CALL READT(LUCSOL,NLMCU,WORK(1)) 655 IF (IPRINT.GT.15) THEN 656 WRITE(LUPRI,*) 'Common TLM' 657 LM = 0 658 DO L = 0,LMAXCU 659 DO M = -L,L 660 LM = LM + 1 661 WRITE(LUPRI,*) 'LM,TLM:',LM,CCTLM(LM) 662 ENDDO 663 ENDDO 664 ENDIF 665C 666 CALL CCSL_TLMGET 667C 668 IF (IPRINT.GT.15) THEN 669 WRITE(LUPRI,*) 'TLM from FIL' 670 LM = 0 671 DO L = 0,LMAXCU 672 DO M = -L,L 673 LM = LM + 1 674 WRITE(LUPRI,*) 'LM,TLM:',LM,CCTLM(LM) 675 ENDDO 676 ENDDO 677 ENDIF 678C 679 LM = 0 680 DO 520 L = 0,LMAXCU 681 READ (LUCSOL) L1,(ISYTLM(M),M=1,2*L+1) 682 IF (L1 .NE. L) THEN 683 WRITE (LUERR,*) 684 * 'ERROR CCSL_RHSTG: L from LUCSOL not as expected' 685 WRITE (LUERR,*) 'L from 520 loop:',L 686 WRITE (LUERR,*) 'L from LUCSOL :',L1 687 CALL QUIT('ERROR CCSL_RHSTG: L from LUCSOL not as expected') 688 END IF 689C 690 DO 500 M = -L,L 691 LM = LM + 1 692 IF (IPRINT .GE. 15) THEN 693 WRITE(LUPRI,'(/A,2I5/A)') ' --- l, m :',L,M, 694 * ' ====================' 695 WRITE(LUPRI,'(A,I2)') ' Symmetry :',ABS(ISYTLM(L+M+1)) 696 END IF 697 IF (ISYTLM(L+M+1) .NE. 1) THEN 698 READ (LUCSOL) 699 GO TO 500 700 ENDIF 701C 702C-------------------------------- 703C Read T(l,m) in AO basis. 704C-------------------------------- 705C 706 ISYMTLM = ISYTLM(L+M+1) 707 KTLMAO = 1 708 KWRK1 = KTLMAO + N2BST(ISYMTLM) 709 LWRK1 = LWORK - KWRK1 710 IF (LWRK1.LT.0) CALL QUIT( 'Too little work in CCSL_RHSTG, 1') 711C 712 CALL CC_SLVINT(WORK(KTLMAO),WORK(KWRK1),LWRK1,ISYMTLM) 713C 714 IF (IPRINT .GT. 50) THEN 715 CALL AROUND('ccsl_rhstg: Tlm_ao matrix: - cc storage') 716 CALL CC_PRFCKAO(WORK(KTLMAO),ISYMTLM) 717 ENDIF 718C 719C SLV98,OC Take care with sign of FACT! 720C 721 FACT = -2.0D0*GL(L,EPSTCU,RCAVCU)*CCTLM(LM) 722 IF (IPRINT.GT.10) THEN 723 WRITE(LUPRI,*) 'CCSL_RHSTG: GL ',GL(L,EPSTCU,RCAVCU) 724 WRITE(LUPRI,*) 'CCSL_RHSTG: TLM',CCTLM(LM) 725 WRITE(LUPRI,*) 'CCSL_RHSTG: FAC',FACT 726 ENDIF 727C 728C----------------------------------------------------- 729C Add contribution to effective AO fock matrix. 730C----------------------------------------------------- 731C 732 CALL DAXPY(N2BST(ISYMTLM),FACT,WORK(KTLMAO),1,FOCK,1) 733C 734 500 CONTINUE 735 520 CONTINUE 736C 737 CALL GPCLOSE(LUCSOL,'KEEP') 738 END 739 740C /* Deck ccsl_ltrb */ 741 SUBROUTINE CCSL_LTRB(RHO1,RHO2,CTR1,CTR2,ISYMTR,LR,WORK,LWORK) 742C 743C----------------------------------------------------------------------------- 744C 745C Purpose: Calculation of Coupled Cluster solvent T^gB contribution 746C to left and right Jacobian transformation. 747C <mu|exp(-T)T^g|CC> for LR = 0 748C F transformation for LR = F 749C P transformation for LR = P 750C 751C LR = 'L','R','0','F','P' 752C SLV98,OC 753C Ove Christiansen, April/mai 1998. 754C 755C----------------------------------------------------------------------------- 756C 757#include "implicit.h" 758#include "maxorb.h" 759#include "mxcent.h" 760 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, IZERO = 0 , TWO = 2.0D0) 761 PARAMETER (HALF = 0.5D0) 762#include "dummy.h" 763#include "iratdef.h" 764#include "priunit.h" 765#include "ccorb.h" 766#include "ccsdsym.h" 767#include "ccsdio.h" 768#include "ccsdinp.h" 769#include "ccfield.h" 770#include "exeinf.h" 771#include "ccfdgeo.h" 772#include "ccslvinf.h" 773#include "ccinftap.h" 774C 775 DIMENSION WORK(LWORK),RHO1(*),RHO2(*),CTR1(*),CTR2(*) 776C 777 CHARACTER*8 LABEL,LIST*(2),LR*(1) 778 CHARACTER*8 FILMME, FILMMX 779 LOGICAL LEXIST 780C 781 DOUBLE PRECISION GL 782 EXTERNAL GL 783C 784 IF (IPRINT.GT.10) THEN 785 WRITE(LUPRI,*)'CCSL_LTRB: Solvent contribution to CC L. transf.' 786 WRITE(LUPRI,*)'CCSL_LTRB: LWORK:', LWORK 787 WRITE(LUPRI,*)'CCSL_LTRB: LR:', LR 788 WRITE(LUPRI,*)'CCSL_LTRB: ISYMTR:', ISYMTR 789 ENDIF 790 IF ( DEBUG .OR.(IPRINT.GT.10)) THEN 791 RHO1N = DDOT(NT1AM(ISYMTR),RHO1,1,RHO1,1) 792 RHO2N = DDOT(NT2AM(ISYMTR),RHO2,1,RHO2,1) 793 WRITE(LUPRI,*) ' Norm af RHO1 in CCSL_LTGB on input:', RHO1N 794 WRITE(LUPRI,*) ' Norm af RHO2 in CCSL_LTGB on input:', RHO2N 795 IF (LR .NE. '0') THEN 796 RHO1N = DDOT(NT1AM(ISYMTR),CTR1,1,CTR1,1) 797 RHO2N = DDOT(NT2AM(ISYMTR),CTR2,1,CTR2,1) 798 WRITE(LUPRI,*) ' Norm af C1AM in CCSL_LTGB on input:', RHO1N 799 WRITE(LUPRI,*) ' Norm af C2AM in CCSL_LTGB on input:', RHO2N 800 END IF 801 ENDIF 802C 803C Note if CCSAV_LAM does not exist then 804C we have no contribution yet. 805C 806c INQUIRE(FILE='CCSAV_LAM',EXIST=LEXIST) 807c IF (.NOT.LEXIST) THEN 808c WRITE(LUPRI,*) ' CCSAV_LAM does not exits yet - no T^gB cont' 809c RETURN 810c ENDIF 811C 812C--------------------- 813C Init parameters. 814C--------------------- 815C 816 NTAMP1 = NT1AM(ISYMTR) 817 NTAMP2 = NT2AM(ISYMTR) 818 IF (CCS) NTAMP2 = 0 819 NTAMP = NTAMP1 + NTAMP2 820C 821 IF (DISCEX) THEN 822 LUMMET = -1 823 LUMMXI = -1 824 FILMME = 'CCSL_ETA' 825 FILMMX = 'CCSL__XI' 826 CALL WOPEN2(LUMMET, FILMME, 64, 0) 827 CALL WOPEN2(LUMMXI, FILMMX, 64, 0) 828 END IF 829C 830C--------------------------------- 831C Read in integrals from file. 832C--------------------------------- 833C 834 LUCSOL = -1 835 CALL GPOPEN(LUCSOL,'AOSOLINT','OLD',' ','UNFORMATTED',IDUMMY, 836 & .FALSE.) 837 REWIND(LUCSOL) 838 CALL MOLLAB('SOLVRLM ',LUCSOL,LUERR) 839 READ (LUCSOL) 840 CALL READT(LUCSOL,NLMCU,WORK(1)) 841C 842 KTGB = 1 843 KWRK1 = KTGB + N2BST(ISYMTR) 844 LWRK1 = LWORK - KWRK1 845 IF (LWRK1.LT.0) CALL QUIT( 'Too little work in CCSL_LTRB, 1') 846 KWRK2 = KWRK1 847 LWRK2 = LWRK1 848C 849 CALL DZERO(WORK(KTGB),N2BST(ISYMTR)) 850C 851 LM = 0 852 IADR1 = 1 853 IADR2 = 1 854C 855 DO 600 L = 0,LMAXCU 856 READ (LUCSOL) L1,(ISYTLM(M),M=1,2*L+1) 857 IF (L1 .NE. L) THEN 858 WRITE (LUERR,*) 859 * 'ERROR CCSL_LTRB: L from LUCSOL not as expected' 860 WRITE (LUERR,*) 'L from 600 loop:',L 861 WRITE (LUERR,*) 'L from LUCSOL :',L1 862 CALL QUIT('ERROR CCSL_LTRB: L from LUCSOL not as expected') 863 END IF 864C 865 DO 500 M = -L,L 866 LM = LM + 1 867 IF (IPRINT .GE. 15) THEN 868 WRITE(LUPRI,'(/A,2I5/A)') ' --- l, m :',L,M, 869 * ' ====================' 870 WRITE(LUPRI,'(A,I2)') ' Symmetry :',ABS(ISYTLM(L+M+1)) 871 END IF 872C 873 ISYMNY = ABS(ISYTLM(L+M+1)) 874 NTA1 = NT1AM(ISYMNY) 875 NTA2 = NT2AM(ISYMNY) 876 NTA = NTA1 + NTA2 877C 878C---------------------------------------------------------- 879C Symmetry should be equal to trial vector symmetry. 880C---------------------------------------------------------- 881C 882 IF (ISYTLM(L+M+1) .NE. ISYMTR) THEN 883 READ (LUCSOL) 884 IADR1 = IADR1 + NTA 885 IADR2 = IADR2 + NTA 886 GO TO 500 887 ENDIF 888C 889C------------------------------------------------------------------- 890C Calculate T^{gB}= -sum_lm 2G_l(ep)<B|T_lm|CC>T_lm operator. 891C 1. Readin integrals again. 892C 2. Calculate xi^T_lm 893C 3. Contract with B 894C 4. scale integrals with 2G_l(ep)<B|T_lm|CC> 895C 5. Add to T^{gB} integrals. 896C------------------------------------------------------------------- 897C 898C 899C-------------------------------- 900C 1. Readin integrals again. 901C Read T(l,m) in AO basis. 902C-------------------------------- 903C 904 KTLMAO = KWRK1 905 KWRK2 = KTLMAO + N2BST(ISYMTR) 906 LWRK2 = LWORK - KWRK2 907 IF (LWRK2.LT.0) CALL QUIT( 'Too little work in CCSL_LTRB, 2') 908C 909 CALL CC_SLVINT(WORK(KTLMAO),WORK(KWRK2),LWRK2,ISYMTR) 910C 911 IF (IPRINT .GT. 25) THEN 912 CALL AROUND('ccsl_ltrb: Tlm_ao matrix: - cc storage') 913 CALL CC_PRFCKAO(WORK(KTLMAO),ISYMTR) 914 ENDIF 915C 916C----------------------------------------------------------- 917C 2. Calculate xi^T_lm (for LR=R actually eta^T_lm ) 918C----------------------------------------------------------- 919C 920 KXI = KWRK2 921 KWRK3 = KXI + NTAMP 922 LWRK3 = LWORK - KWRK3 923 IF (LWRK3.LT.0) CALL QUIT( 'Too little work in CCSL_LTRB, 3') 924C 925 IF (.NOT. (DISCEX)) THEN 926 LABEL = 'GIVE INT' 927 IF ((LR.EQ.'L').OR.(LR.EQ.'0').OR.(LR.EQ.'P')) THEN 928 CALL CC_XKSI(WORK(KXI),LABEL,ISYMTR,0,WORK(KTLMAO), 929 * WORK(KWRK3),LWRK3) 930 ELSE IF ((LR.EQ.'R').OR.(LR.EQ.'F')) THEN 931 LIST = 'L0' 932 ILSTNR = 1 933 CALL CC_ETAC(ISYMTR,LABEL,WORK(KXI), 934 * LIST,ILSTNR,0,WORK(KTLMAO),WORK(KWRK3),LWRK3) 935 ENDIF 936 ELSE 937 IF ((LR.EQ.'L').OR.(LR.EQ.'P')) THEN 938 CALL GETWA2(LUMMXI,FILMMX,WORK(KXI),IADR2,NTAMP) 939 IADR2 = IADR2 + NTAMP 940 ELSE IF ((LR.EQ.'R').OR.(LR.EQ.'F')) THEN 941 LIST = 'L0' 942 ILSTNR = 1 943 CALL GETWA2(LUMMET,FILMME,WORK(KXI),IADR1,NTAMP) 944 IADR1 = IADR1 + NTAMP 945 ELSE IF (LR.EQ.'0') THEN 946 LABEL = 'GIVE INT' 947 CALL CC_XKSI(WORK(KXI),LABEL,ISYMTR,0,WORK(KTLMAO), 948 * WORK(KWRK3),LWRK3) 949 END IF 950 END IF 951C 952C--------------------------- 953C 3. Contract with B 954C--------------------------- 955C 956 IF (LR.NE.'0') THEN 957 KXI1 = KXI 958 KXI2 = KXI + NTAMP1 959C 960 BXILMD1= DDOT(NTAMP1,CTR1,1,WORK(KXI1),1) 961 BXILMD2= DDOT(NTAMP2,CTR2,1,WORK(KXI2),1) 962 BXILMD = BXILMD1 + BXILMD2 963C 964C---------------------------------------------------- 965C 4. Find 2G_l(ep)<B|T_lm|CC> factor 966C---------------------------------------------------- 967C 968 GLFAC = GL(L,EPOPCU,RCAVCU) 969C 970C SLV98,OC 971C 972 FACT = 2.0D0*GLFAC*BXILMD 973C 974 IF (IPRINT.GT.10) THEN 975 WRITE(LUPRI,*) 'CCSL_LTRB: GL ',GLFAC 976 WRITE(LUPRI,*) 'CCSL_LTRB: BTLM1 ',BXILMD1 977 WRITE(LUPRI,*) 'CCSL_LTRB: BTLM2 ',BXILMD2 978 WRITE(LUPRI,*) 'CCSL_LTRB: BTLM ',BXILMD 979 WRITE(LUPRI,*) 'CCSL_LTRB: FAC ',FACT 980 ENDIF 981C 982C-------------------------------------------------------------- 983C 5. Add to T^{gB} integrals.(for LR=R actually T^gC) 984C-------------------------------------------------------------- 985C 986 CALL DAXPY(N2BST(ISYMTR),FACT,WORK(KTLMAO),1, 987 * WORK(KTGB),1) 988 989 ELSE IF (LR.EQ.'0') THEN 990 KXI1 = KXI 991 KXI2 = KXI + NTAMP1 992 FACT = -2.0D0*GL(L,EPSTCU,RCAVCU)*CCTLM(LM) 993 CALL DAXPY(NT1AM(ISYMTR),FACT,WORK(KXI1),1, 994 * RHO1,1) 995 CALL DAXPY(NT2AM(ISYMTR),FACT,WORK(KXI2),1, 996 * RHO2,1) 997 ENDIF 998C 999 500 CONTINUE 1000 600 CONTINUE 1001C 1002 IF (IPRINT .GT. 50) THEN 1003 CALL AROUND('ccsl_ltrb: T^gB_ao matrix: ') 1004 CALL CC_PRFCKAO(WORK(KTGB),ISYMTR) 1005 ENDIF 1006C 1007 IF (DISCEX) THEN 1008 CALL WCLOSE2(LUMMET, FILMME, 'KEEP') 1009 CALL WCLOSE2(LUMMXI, FILMMX, 'KEEP') 1010 END IF 1011C 1012 IF (LR.NE.'0') THEN 1013C 1014C----------------------------------------------------- 1015C Calculate contribution from the T^gB operator. 1016C----------------------------------------------------- 1017C 1018 KETA = KWRK2 1019 KWRK4 = KETA + NTAMP 1020 LWRK4 = LWORK - KWRK4 1021 IF (LWRK4.LT.0) CALL QUIT( 'Too little work in CCSL_LTRB 4') 1022 KETA1 = KETA 1023 KETA2 = KETA + NTAMP1 1024C 1025 IF ((LR.EQ.'L').OR.(LR.EQ.'F')) THEN 1026 LIST = 'L0' 1027 LABEL = 'GIVE INT' 1028 CALL CC_ETAC(ISYMTR,LABEL,WORK(KETA), 1029 * LIST,1,0,WORK(KTGB),WORK(KWRK4),LWRK4) 1030 ELSE IF ((LR.EQ.'R').OR.(LR.EQ.'P')) THEN 1031 LABEL = 'GIVE INT' 1032 CALL CC_XKSI(WORK(KETA),LABEL,ISYMTR,0,WORK(KTGB), 1033 * WORK(KWRK4),LWRK4) 1034 IF (LR.EQ.'R') 1035 * CALL CCLR_DIASCL(WORK(KETA2),TWO,ISYMTR) 1036 ENDIF 1037C 1038 IF ( DEBUG .OR.(IPRINT.GT.10)) THEN 1039 RHO1N = DDOT(NT1AM(ISYMTR),WORK(KETA1),1,WORK(KETA1),1) 1040 RHO2N = DDOT(NT2AM(ISYMTR),WORK(KETA2),1,WORK(KETA2),1) 1041 WRITE(LUPRI,*) ' Norm af T^gB contribution to LHTR1:', RHO1N 1042 WRITE(LUPRI,*) ' Norm af T^gB contribution to LHTR2:', RHO2N 1043 ENDIF 1044C 1045 CALL DAXPY(NT1AM(ISYMTR),ONE,WORK(KETA1),1,RHO1,1) 1046 CALL DAXPY(NT2AM(ISYMTR),ONE,WORK(KETA2),1,RHO2,1) 1047 IF ( DEBUG .OR.(IPRINT.GT.10)) THEN 1048 RHO1N = DDOT(NT1AM(ISYMTR),RHO1,1,RHO1,1) 1049 RHO2N = DDOT(NT2AM(ISYMTR),RHO2,1,RHO2,1) 1050 WRITE(LUPRI,*) ' Norm af RHO1 in CCSL_LTGB:', RHO1N 1051 WRITE(LUPRI,*) ' Norm af RHO2 in CCSL_LTGB:', RHO2N 1052 ENDIF 1053C 1054 ENDIF 1055C 1056 CALL GPCLOSE(LUCSOL,'KEEP') 1057 END 1058C /* Deck ccsl_tlmget*/ 1059 SUBROUTINE CCSL_TLMGET 1060#include "implicit.h" 1061#include "maxorb.h" 1062#include "mxcent.h" 1063 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, IZERO = 0 , TWO = 2.0D0) 1064 PARAMETER (HALF = 0.5D0) 1065#include "dummy.h" 1066#include "iratdef.h" 1067#include "priunit.h" 1068#include "ccslvinf.h" 1069#include "ccinftap.h" 1070C 1071 LUCTLM = -1 1072 CALL GPOPEN(LUCTLM,'CC_TL','UNKNOWN',' ','FORMATTED',IDUMMY, 1073 & .FALSE.) 1074 LM = 0 1075 DO 820 L = 0,LMAXCU 1076 DO 800 M = -L,L 1077 LM = LM + 1 1078 READ(LUCTLM,'(I5,E25.15)',END=100,ERR=100) LM1,CCTLM(LM) 1079 IF (LM1.NE.LM) CALL QUIT( 'CCSL_TLMGET: Error on CC_TL') 1080 GOTO 800 1081 100 CCTLM(LM) = 0.0D0 1082C 1083 800 CONTINUE 1084 820 CONTINUE 1085C 1086 CALL GPCLOSE(LUCTLM,'KEEP') 1087 RETURN 1088C 1089 END 1090C 1091C /* Deck ccsl_tlmput*/ 1092 SUBROUTINE CCSL_TLMPUT 1093#include "implicit.h" 1094#include "maxorb.h" 1095#include "mxcent.h" 1096 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, IZERO = 0 , TWO = 2.0D0) 1097 PARAMETER (HALF = 0.5D0) 1098#include "dummy.h" 1099#include "iratdef.h" 1100#include "priunit.h" 1101#include "ccslvinf.h" 1102#include "ccinftap.h" 1103C 1104 LUCTLM = -1 1105 CALL GPOPEN(LUCTLM,'CC_TL','UNKNOWN',' ','FORMATTED',IDUMMY, 1106 & .FALSE.) 1107C 1108 LM = 0 1109 DO 920 L = 0,LMAXCU 1110 DO 900 M = -L,L 1111 LM = LM + 1 1112 WRITE(LUCTLM,'(I5,E25.15)') LM,CCTLM(LM) 1113 900 CONTINUE 1114 920 CONTINUE 1115C 1116 CALL GPCLOSE(LUCTLM,'KEEP') 1117C 1118 END 1119**************************************************************** 1120C /* Deck ccsl_pbtr */ 1121 SUBROUTINE CCSL_PBTR(RHO1,RHO2,ISYRES, 1122 * LISTL,IDLSTL,CTR1,ISYCTR, 1123 * LISTR,IDLSTR,BTR1,ISYBTR, 1124 * MODEL,WORK,LWORK) 1125C 1126C--------------------------------------------------------------- 1127C Calculates contributions from the T^gB , the ^CT^g 1128C and ^CT^gB operators. 1129C JK+OC, 03 1130C--------------------------------------------------------------- 1131C 1132#include "implicit.h" 1133#include "maxorb.h" 1134#include "mxcent.h" 1135 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, IZERO = 0 , TWO = 2.0D0) 1136 PARAMETER (HALF = 0.5D0) 1137#include "dummy.h" 1138#include "iratdef.h" 1139#include "priunit.h" 1140#include "ccorb.h" 1141#include "ccsdsym.h" 1142#include "ccsdio.h" 1143#include "ccsdinp.h" 1144#include "ccfield.h" 1145#include "exeinf.h" 1146#include "ccfdgeo.h" 1147#include "ccslvinf.h" 1148#include "ccinftap.h" 1149C 1150 DIMENSION WORK(LWORK),RHO1(*),RHO2(*),CTR1(*),BTR1(*) 1151C 1152 CHARACTER*(*) LISTL,LISTR,LIST*(2) 1153 CHARACTER*8 LABEL 1154 CHARACTER*10 MODEL 1155 LOGICAL LEXIST 1156 INTEGER IDLSTL, IDLSTR, ISYCTR, ISYBTR, ISYBC 1157C 1158 INTEGER KDUM 1159 PARAMETER( KDUM = +99 999 999 ) ! dummy address for work space 1160C 1161 IF (IPRINT.GT.10) THEN 1162 WRITE(LUPRI,*)'CCSL_PBTR: ISYRES =',ISYRES 1163 WRITE(LUPRI,*)'CCSL_PBTR: ISYCTR =',ISYCTR 1164 WRITE(LUPRI,*)'CCSL_PBTR: ISYBTR =',ISYBTR 1165 WRITE(LUPRI,*)'CCSL_PBTR: LISTL =',LISTL 1166 WRITE(LUPRI,*)'CCSL_PBTR: LISTR =',LISTR 1167 WRITE(LUPRI,*)'CCSL_PBTR: IDLSTL =',IDLSTL 1168 WRITE(LUPRI,*)'CCSL_PBTR: IDLSTR =',IDLSTR 1169 CALL FLSHFO(LUPRI) 1170 ENDIF 1171C 1172C--------------------- 1173C Init parameters. 1174C--------------------- 1175C 1176C For the B (right) trial vector 1177C 1178 NAMP1B = NT1AM(ISYBTR) 1179 NAMP2B = NT2AM(ISYBTR) 1180 NAMPB = NAMP1B + NAMP2B 1181C 1182C For the C (left) trial vector 1183C 1184 NAMP1C = NT1AM(ISYCTR) 1185 NAMP2C = NT2AM(ISYCTR) 1186 NAMPC = NAMP1C + NAMP2C 1187C 1188C For the F = C X B vector 1189C 1190 ISYBC = MULD2H(ISYBTR,ISYCTR) 1191 NAMP1F = NT1AM(ISYBC) 1192 NAMP2F = NT2AM(ISYBC) 1193 NAMPF = NAMP1F + NAMP2F 1194C 1195C----------------------------------------------------------- 1196C Calculate contribution from the effective operators. 1197C----------------------------------------------------------- 1198C 1199C First we concentrate on the effective operator T^gB 1200C ----------------------------------------------------- 1201C 1202 KTGB = 1 1203 KWRK1 = KTGB + N2BST(ISYBTR) 1204 LWRK1 = LWORK - KWRK1 1205 IF (LWRK1.LT.0) CALL QUIT( 'Too little work in CCSL_PBTR, 1') 1206C 1207 CALL DZERO(WORK(KTGB),N2BST(ISYBTR)) 1208C 1209 CALL CCSL_TGB(BTR1,ISYBTR,LISTR,IDLSTR,WORK(KTGB),'ET', 1210 * MODEL,WORK(KWRK1),LWRK1) 1211C 1212C Symmetry: 1213C --------- 1214C 1215 ISYBC = MULD2H(ISYBTR,ISYCTR) 1216C 1217 IF (ISYBC .NE. ISYRES) THEN 1218 CALL QUIT( 'Symmetry problem in CCSL_PBTR') 1219 END IF 1220C 1221 KETA = KWRK1 1222 KWRK2 = KETA + NAMPF 1223 LWRK2 = LWORK - KWRK2 1224 IF (LWRK2.LT.0) CALL QUIT( 'Too little work in CCSL_PBTR, 2') 1225C 1226C Note, LISTL .EQ. LE/L1 so the HF part of the following 1227C eta-transformation is skipped 1228C 1229 LABEL = 'GIVE INT' 1230 CALL CC_ETAC(ISYBTR,LABEL,WORK(KETA), 1231 * LISTL,IDLSTL,0,WORK(KTGB),WORK(KWRK2),LWRK2) 1232C 1233 KETA1 = KETA 1234 KETA2 = KETA1 + NAMP1F 1235C 1236 CALL DAXPY(NAMP1F,ONE,WORK(KETA1),1,RHO1,1) 1237 CALL DAXPY(NAMP2F,ONE,WORK(KETA2),1,RHO2,1) 1238C 1239C Next, the effective operator ^CT^gB is considered 1240C---------------------------------------------------------- 1241C 1242C Symmetry of the operator: 1243C ------------------------- 1244C 1245 ISYBC = MULD2H(ISYBTR,ISYCTR) 1246C 1247 KTGB = 1 1248 KWRK1 = KTGB + N2BST(ISYBC) 1249 LWRK1 = LWORK - KWRK1 1250 IF (LWRK1.LT.0) CALL QUIT( 'Too little work in CCSL_PBTR, 3') 1251C 1252 CALL DZERO(WORK(KTGB),N2BST(ISYBC)) 1253C 1254 CALL CCSL_CTGB(LISTL,IDLSTL,CTR1,ISYCTR, 1255 * LISTR,IDLSTR,BTR1,ISYBTR, 1256 * WORK(KTGB),MODEL,WORK(KWRK1),LWRK1) 1257C 1258 KETA = KWRK1 1259 KWRK2 = KETA + NAMPF 1260 LWRK2 = LWORK - KWRK2 1261 IF (LWRK2.LT.0) CALL QUIT( 'Too little work in CCSL_PBTR, 4') 1262C 1263 CALL DZERO(WORK(KETA),NAMPF) 1264C 1265 LABEL = 'GIVE INT' 1266 LIST = 'L0' 1267 IDLINO = 1 1268C 1269 CALL CC_ETAC(ISYBC,LABEL,WORK(KETA), 1270 * LIST,IDLINO,0,WORK(KTGB),WORK(KWRK2),LWRK2) 1271C 1272 KETA1 = KETA 1273 KETA2 = KETA1 + NAMP1F 1274C 1275 CALL DAXPY(NT1AM(ISYBC),ONE,WORK(KETA1),1,RHO1,1) 1276 CALL DAXPY(NT2AM(ISYBC),ONE,WORK(KETA2),1,RHO2,1) 1277C 1278C Finally, we calculate the third contribution which is 1279C a contribution from a T^gB operator but calculated 1280C as a xi vector element. 1281C----------------------------------------------------------- 1282C 1283 KTGB = 1 1284 KWRK1 = KTGB + N2BST(ISYCTR) 1285 LWRK1 = LWORK - KWRK1 1286 IF (LWRK1.LT.0) CALL QUIT( 'Too little work in CCSL_PBTR, 5') 1287C 1288 CALL DZERO(WORK(KTGB),N2BST(ISYCTR)) 1289C 1290 CALL CCSL_TGB(CTR1,ISYCTR,LISTL,IDLSTL,WORK(KTGB),'XI', 1291 * MODEL,WORK(KWRK1),LWRK1) 1292C 1293C Symmetry: 1294C --------- 1295C 1296 ISYBC = MULD2H(ISYBTR,ISYCTR) 1297C 1298 IF (ISYBC .NE. ISYRES) THEN 1299 CALL QUIT( 'Symmetry problem in CCSL_PBTR') 1300 END IF 1301 1302 LABEL = 'GIVE INT' 1303 LIST = 'L0' 1304 LISTNO = 1 1305C 1306C (NB: Result vector from the following transformation is 1307C given at the beginning of WORK(KWRK1).) 1308C 1309 CALL CCLR_FA(LABEL,ISYCTR,LISTR,IDLSTR, 1310 & LIST,LISTNO,WORK(KTGB),WORK(KWRK1),LWRK1) 1311C 1312C Jacob, Husk at DDOTog DAXPY virker paa denne her maade !! 1313C TEMP1 = DDOT(NT1AM(ISYBC),WORK(KWRK1),1,WORK(KWRK1),1) 1314C TEMP2 = DDOT(NT2AM(ISYBC),WORK(KWRK1+NT1AM(ISYBC)),1, 1315C * WORK(KWRK1+NT1AM(ISYBC)),1) 1316C 1317 CALL DAXPY(NT1AM(ISYBC),ONE,WORK(KWRK1),1,RHO1,1) 1318 CALL DAXPY(NT2AM(ISYBC),ONE, 1319 * WORK(KWRK1+NT1AM(ISYBC)),1,RHO2,1) 1320C 1321 END 1322**************************************************************** 1323C /* Deck ccsl_ctgb */ 1324 SUBROUTINE CCSL_CTGB(LISTL,IDLSTL,CTR1,ISYCTR, 1325 * LISTR,IDLSTR,BTR1,ISYBTR, 1326 * TGB,MODEL,WORK,LWORK) 1327C 1328C JK+OC, 03 1329C----------------------------------------------------------------------------- 1330C 1331C Construcs effective operator equal to 1332C ^CT^gB = -2 Sum_lm g_l * Sum_sigma Sum_my 1333C t^C_my <bar my|[T_lm,Tau_sigma]|CC> t^B_sigma T_lm 1334C 1335C----------------------------------------------------------------------------- 1336C 1337#include "implicit.h" 1338#include "maxorb.h" 1339#include "mxcent.h" 1340 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, IZERO = 0 , TWO = 2.0D0) 1341 PARAMETER (HALF = 0.5D0) 1342#include "dummy.h" 1343#include "iratdef.h" 1344#include "priunit.h" 1345#include "ccorb.h" 1346#include "ccsdsym.h" 1347#include "ccsdio.h" 1348#include "ccsdinp.h" 1349#include "ccfield.h" 1350#include "exeinf.h" 1351#include "ccfdgeo.h" 1352#include "ccslvinf.h" 1353#include "ccinftap.h" 1354C 1355 DIMENSION WORK(LWORK),BTR1(*),CTR1(*) 1356 DIMENSION TGB(*) 1357C 1358 CHARACTER*8 LABEL, LIST*(2) 1359 CHARACTER*(*) LISTL, LISTR 1360 CHARACTER*10 MODEL 1361 INTEGER IDLSTL, IDLSTR, ISYCTR, ISYBTR 1362 LOGICAL LEXIST 1363C 1364 DOUBLE PRECISION GL 1365 EXTERNAL GL 1366C 1367 IF (IPRINT.GT.10) THEN 1368 WRITE(LUPRI,*)'CCSL_PBTR: ISYCTR =',ISYCTR 1369 WRITE(LUPRI,*)'CCSL_PBTR: ISYBTR =',ISYBTR 1370 WRITE(LUPRI,*)'CCSL_PBTR: LISTL =',LISTL 1371 WRITE(LUPRI,*)'CCSL_PBTR: LISTR =',LISTR 1372 WRITE(LUPRI,*)'CCSL_PBTR: IDLSTL =',IDLSTL 1373 WRITE(LUPRI,*)'CCSL_PBTR: IDLSTR =',IDLSTR 1374 CALL FLSHFO(LUPRI) 1375 ENDIF 1376C 1377C--------------------- 1378C Init parameters. 1379C--------------------- 1380C 1381C For the B (right) trial vector 1382C 1383 NAMP1B = NT1AM(ISYBTR) 1384 NAMP2B = NT2AM(ISYBTR) 1385 NAMPB = NAMP1B + NAMP2B 1386C 1387C For the C (left) trial vector 1388C 1389 NAMP1C = NT1AM(ISYCTR) 1390 NAMP2C = NT2AM(ISYCTR) 1391 NAMPC = NAMP1C + NAMP2C 1392C 1393C For the F = B x C vector 1394C 1395 ISYBC = MULD2H(ISYBTR,ISYCTR) 1396 NAMP1F = NT1AM(ISYBC) 1397 NAMP2F = NT2AM(ISYBC) 1398 NAMPF = NAMP1F + NAMP2F 1399C 1400C----------------------------------------------- 1401C Readin response vector (B, right) from file 1402C----------------------------------------------- 1403C 1404 KT2AMPA = 1 1405 KWRK1 = KT2AMPA + NT2AM(ISYBTR) 1406 LWRK1 = LWORK - KWRK1 1407C 1408 IF (LWRK1 .LT. 0) THEN 1409 CALL QUIT('Insuff. work in CCSL_TGB') 1410 END IF 1411C 1412 IOPT = 2 1413 CALL CC_RDRSP(LISTR,IDLSTR,ISYBTR,IOPT,MODEL, 1414 * WORK(KWRK1),WORK(KT2AMPA)) 1415C 1416C First we perform the contraction of of the eta' 1417C element from the left hand side. The ' means that 1418C the HF contribution is skipped. 1419C 1420C--------------------------------- 1421C Read in integrals from file. 1422C--------------------------------- 1423C 1424 LUCSOL = -1 1425 CALL GPOPEN(LUCSOL,'AOSOLINT','OLD',' ','UNFORMATTED',IDUMMY, 1426 & .FALSE.) 1427 REWIND(LUCSOL) 1428 CALL MOLLAB('SOLVRLM ',LUCSOL,LUERR) 1429 READ (LUCSOL) 1430 CALL READT(LUCSOL,NLMCU,WORK(KWRK1)) 1431C 1432 LM = 0 1433 DO 600 L = 0,LMAXCU 1434 READ (LUCSOL) L1,(ISYTLM(M),M=1,2*L+1) 1435 IF (L1 .NE. L) THEN 1436 WRITE (LUERR,*) 1437 * 'ERROR CCSL_CTGB: L from LUCSOL not as expected' 1438 WRITE (LUERR,*) 'L from 600 loop:',L 1439 WRITE (LUERR,*) 'L from LUCSOL :',L1 1440 CALL QUIT('ERROR CCSL_CTGB: L from LUCSOL not as expected') 1441 END IF 1442C 1443 DO 500 M = -L,L 1444 LM = LM + 1 1445 IF (IPRINT .GE. 15) THEN 1446 WRITE(LUPRI,'(/A,2I5/A)') ' --- l, m :',L,M, 1447 * ' ====================' 1448 WRITE(LUPRI,'(A,I2)') ' Symmetry :',ABS(ISYTLM(L+M+1)) 1449 END IF 1450C 1451C---------------------------------------------- 1452C Symmetry should be equal to symmetry 1453C of the final contracted vector (ISYBC). 1454C---------------------------------------------- 1455C 1456 ISYBC = MULD2H(ISYBTR,ISYCTR) 1457C 1458 IF (ISYTLM(L+M+1) .NE. ISYBC) THEN 1459 READ (LUCSOL) 1460 GO TO 500 1461 ENDIF 1462C 1463C-------------------------------- 1464C Read T(l,m) in AO basis. 1465C-------------------------------- 1466C 1467 KTLMAO = KWRK1 1468 KWRK2 = KTLMAO + N2BST(ISYBC) 1469 LWRK2 = LWORK - KWRK2 1470 IF (LWRK2.LT.0) CALL QUIT( 'Too little work in CCSL_CTGB, 2') 1471C 1472 CALL CC_SLVINT(WORK(KTLMAO),WORK(KWRK2),LWRK2,ISYBC) 1473C 1474 IF (IPRINT .GT. 25) THEN 1475 CALL AROUND('CCSL_CTGB : Tlm_ao matrix: - cc storage') 1476 CALL CC_PRFCKAO(WORK(KTLMAO),ISYBC) 1477 ENDIF 1478C 1479C--------------------------------------------------------------------- 1480C Calculate eta^T_lm 1481C (LIST = 'LE' ensures that the HF part of eta^T_lm is skipped. 1482C--------------------------------------------------------------------- 1483C 1484 KXI = KWRK2 1485 KWRK3 = KXI + NAMPB 1486 LWRK3 = LWORK - KWRK3 1487 IF (LWRK3.LT.0) CALL QUIT( 'Too little work in CCSL_CTGB, 3') 1488C 1489 LABEL = 'GIVE INT' 1490 CALL CC_ETAC(ISYBC,LABEL,WORK(KXI), 1491 * LISTL,IDLSTL,0,WORK(KTLMAO),WORK(KWRK3),LWRK3) 1492C 1493C--------------------------------------------------------------- 1494C Contract with trial vector B (from the right hand side) 1495C--------------------------------------------------------------- 1496C 1497 KXI1 = KXI 1498 KXI2 = KXI + NT1AM(ISYBTR) 1499C 1500 BXILMD1= DDOT(NT1AM(ISYBTR),BTR1,1,WORK(KXI1),1) 1501 BXILMD2= DDOT(NT2AM(ISYBTR),WORK(KT2AMPA),1,WORK(KXI2),1) 1502 BXILMD = BXILMD1 + BXILMD2 1503C 1504C-------------------------------------------------------------- 1505C Find factor from contracted vector and G_lm factor 1506C-------------------------------------------------------------- 1507C 1508 GLFAC = GL(L,EPOPCU,RCAVCU) 1509 FACT = 2.0D0*GLFAC*BXILMD 1510C 1511C-------------------------------------------------------------- 1512C Put factor on integrals and (for a given l,m) the 1513C effective operator has now been constructed. 1514C-------------------------------------------------------------- 1515C 1516 CALL DAXPY(N2BST(ISYBC),FACT,WORK(KTLMAO),1, 1517 * TGB,1) 1518C--------------------------------------------------------------- 1519 1520 500 CONTINUE 1521 600 CONTINUE 1522C 1523 IF (IPRINT .GT. 50) THEN 1524 CALL AROUND('ccsl_ctgb : T^gB_ao matrix: ') 1525 CALL CC_PRFCKAO(TGB,ISYBC) 1526 ENDIF 1527C 1528 CALL GPCLOSE(LUCSOL,'KEEP') 1529 END 1530******************************************************************* 1531C /* Deck ccsl_tgb */ 1532 SUBROUTINE CCSL_TGB(CTR1,ISYMTR,LISTIN,IDLIST,TGB,LR, 1533 * MODEL,WORK,LWORK) 1534C 1535C JK+OC, 03 1536C----------------------------------------------------------------------------- 1537C 1538C IF (LR.EQ.'ET') then the following effective operator is constructed 1539C T^gB = -2 Sum_lm g_l * Sum_sigma 1540C <Lambda|[T_lm,Tau_sigma]|CC> t^B_sigma T_lm 1541C 1542C IF (LR.EQ.'XI') then the following effective operator is constructed 1543C ^BT^g = -2 Sum_lm g_l * Sum_sigma 1544C t^B_sigma <bar sigma|T_lm|CC> T_lm 1545C 1546C----------------------------------------------------------------------------- 1547C 1548#include "implicit.h" 1549#include "maxorb.h" 1550#include "mxcent.h" 1551 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, IZERO = 0 , TWO = 2.0D0) 1552 PARAMETER (HALF = 0.5D0) 1553#include "dummy.h" 1554#include "iratdef.h" 1555#include "priunit.h" 1556#include "ccorb.h" 1557#include "ccsdsym.h" 1558#include "ccsdio.h" 1559#include "ccsdinp.h" 1560#include "ccfield.h" 1561#include "exeinf.h" 1562#include "ccfdgeo.h" 1563#include "ccslvinf.h" 1564#include "ccinftap.h" 1565C 1566 DIMENSION WORK(LWORK),CTR1(*) 1567 DIMENSION TGB(*) 1568C 1569 INTEGER ISYMTR, IDLIST 1570 INTEGER KDUM 1571 PARAMETER( KDUM = +99 999 999 ) ! dummy address for work space 1572C 1573 CHARACTER*(*) LISTIN 1574 CHARACTER*10 MODEL 1575C 1576 CHARACTER*8 LABEL,LR*(2),LIST 1577 CHARACTER*(8) FILMME, FILMMX 1578 LOGICAL LEXIST 1579C 1580 DOUBLE PRECISION GL 1581 EXTERNAL GL 1582C 1583 IF (IPRINT.GT.10) THEN 1584 WRITE(LUPRI,*)'CCSL_TGB : ISYMTR:', ISYMTR 1585 WRITE(LUPRI,*)'CCSL_TGB : LISTIN:', LISTIN 1586 ENDIF 1587C 1588C--------------------- 1589C Init parameters. 1590C--------------------- 1591C 1592 NTAMP1 = NT1AM(ISYMTR) 1593 NTAMP2 = NT2AM(ISYMTR) 1594 IF (CCS) NTAMP2 = 0 1595 NTAMP = NTAMP1 + NTAMP2 1596C 1597C 1598 IF (DISCEX) THEN 1599 LUMMET = -1 1600 LUMMXI = -1 1601 FILMME = 'CCSL_ETA' 1602 FILMMX = 'CCSL__XI' 1603 CALL WOPEN2(LUMMET, FILMME, 64, 0) 1604 CALL WOPEN2(LUMMXI, FILMMX, 64, 0) 1605 END IF 1606C 1607C-------------------------------------- 1608C Readin trial vector from file 1609C-------------------------------------- 1610C 1611 KT2AMPA = 1 1612 KWRK1 = KT2AMPA + NT2AM(ISYMTR) 1613 LWRK1 = LWORK - KWRK1 1614C 1615 IF (LWRK1 .LT. 0) THEN 1616 CALL QUIT('Insuff. work in CCSL_TGB') 1617 END IF 1618C 1619 IOPT = 2 1620 CALL CC_RDRSP(LISTIN,IDLIST,ISYMTR,IOPT,MODEL, 1621 * WORK(KDUM),WORK(KT2AMPA)) 1622C 1623C--------------------------------- 1624C Read in integrals from file. 1625C--------------------------------- 1626C 1627 LUCSOL = -1 1628 CALL GPOPEN(LUCSOL,'AOSOLINT','OLD',' ','UNFORMATTED',IDUMMY, 1629 & .FALSE.) 1630 REWIND(LUCSOL) 1631 CALL MOLLAB('SOLVRLM ',LUCSOL,LUERR) 1632 READ (LUCSOL) 1633 CALL READT(LUCSOL,NLMCU,WORK(KWRK1)) 1634C 1635 TEMP = 0.00 1636 LM = 0 1637 IADR1 = 1 1638 IADR2 = 1 1639C 1640 DO 600 L = 0,LMAXCU 1641 READ (LUCSOL) L1,(ISYTLM(M),M=1,2*L+1) 1642 IF (L1 .NE. L) THEN 1643 WRITE (LUERR,*) 1644 * 'ERROR CCSL_TGB: L from LUCSOL not as expected' 1645 WRITE (LUERR,*) 'L from 600 loop:',L 1646 WRITE (LUERR,*) 'L from LUCSOL :',L1 1647 CALL QUIT('ERROR CCSL_TGB: L from LUCSOL not as expected') 1648 END IF 1649C 1650 DO 500 M = -L,L 1651 LM = LM + 1 1652 IF (IPRINT .GE. 15) THEN 1653 WRITE(LUPRI,'(/A,2I5/A)') ' --- l, m :',L,M, 1654 * ' ====================' 1655 WRITE(LUPRI,'(A,I2)') ' Symmetry :',ABS(ISYTLM(L+M+1)) 1656 END IF 1657C 1658 ISYMNY = ABS(ISYTLM(L+M+1)) 1659 NTA1 = NT1AM(ISYMNY) 1660 NTA2 = NT2AM(ISYMNY) 1661 NTA = NTA1 + NTA2 1662C 1663C---------------------------------------------------------- 1664C Symmetry should be equal to trial vector symmetry. 1665C---------------------------------------------------------- 1666C 1667 IF (ISYTLM(L+M+1) .NE. ISYMTR) THEN 1668 READ (LUCSOL) 1669 IADR1 = IADR1 + NTA 1670 IADR2 = IADR2 + NTA 1671 GO TO 500 1672 ENDIF 1673C 1674C-------------------------------- 1675C Read T(l,m) in AO basis. 1676C-------------------------------- 1677C 1678 KTLMAO = KWRK1 1679 KWRK2 = KTLMAO + N2BST(ISYMTR) 1680 LWRK2 = LWORK - KWRK2 1681 IF (LWRK2.LT.0) CALL QUIT( 'Too little work in CCSL_CTGB, 2') 1682C 1683 CALL CC_SLVINT(WORK(KTLMAO),WORK(KWRK2),LWRK2,ISYMTR) 1684C 1685 IF (IPRINT .GT. 25) THEN 1686 CALL AROUND('CCSL_TGB : Tlm_ao matrix: - cc storage') 1687 CALL CC_PRFCKAO(WORK(KTLMAO),ISYMTR) 1688 ENDIF 1689C 1690C--------------------------------------------------------------------- 1691C (LR.EQ.'ET') Calculate eta^T_lm 1692C (LR.EQ.'XI') Calculate xi^T_lm 1693C--------------------------------------------------------------------- 1694C 1695 KXI = KWRK2 1696 KWRK3 = KXI + NTAMP 1697 LWRK3 = LWORK - KWRK3 1698 IF (LWRK3.LT.0) CALL QUIT( 'Too little work in CCSL_TGB, 3') 1699C 1700 IF (.NOT. (DISCEX)) THEN 1701 IF (LR.EQ.'ET') THEN 1702 LABEL = 'GIVE INT' 1703 LIST = 'L0' 1704 IDLINO = 1 1705 CALL CC_ETAC(ISYMTR,LABEL,WORK(KXI), 1706 * LIST,IDLINO,0,WORK(KTLMAO),WORK(KWRK3),LWRK3) 1707C 1708 ELSE IF (LR.EQ.'XI') THEN 1709 LABEL = 'GIVE INT' 1710 CALL CC_XKSI(WORK(KXI),LABEL,ISYMTR,0,WORK(KTLMAO), 1711 * WORK(KWRK3),LWRK3) 1712 END IF 1713 ELSE 1714 IF (LR.EQ.'ET') THEN 1715 CALL GETWA2(LUMMET,FILMME,WORK(KXI),IADR1,NTAMP) 1716 IADR1 = IADR1 + NTAMP 1717 ELSE IF (LR.EQ.'XI') THEN 1718 CALL GETWA2(LUMMXI,FILMMX,WORK(KXI),IADR2,NTAMP) 1719 IADR2 = IADR2 + NTAMP 1720 END IF 1721 END IF 1722C 1723C-------------------------------------------------------------- 1724C Contract with trial vector: 1725C 1726C (If LR.EQ.'ET' then from the right hand side) 1727C (If LR.EQ.'XI' then from the left hand side) 1728C 1729C-------------------------------------------------------------- 1730C 1731 KXI1 = KXI 1732 KXI2 = KXI + NTAMP1 1733C 1734C 1735 CXILMD1= DDOT(NTAMP1,CTR1,1,WORK(KXI1),1) 1736 CXILMD2= DDOT(NTAMP2,WORK(KT2AMPA),1,WORK(KXI2),1) 1737 CXILMD = CXILMD1 + CXILMD2 1738C 1739C-------------------------------------------------------------- 1740C Find factor from contracted vector and G_lm factor 1741C-------------------------------------------------------------- 1742C 1743 GLFAC = GL(L,EPOPCU,RCAVCU) 1744 FACT = 2.0D0*GLFAC*CXILMD 1745C 1746C-------------------------------------------------------------- 1747C Put factor on integrals and (for a given l,m) the 1748C effective operator has now been constructed. 1749C-------------------------------------------------------------- 1750C 1751 CALL DAXPY(N2BST(ISYMTR),FACT,WORK(KTLMAO),1, 1752 * TGB,1) 1753C--------------------------------------------------------------- 1754 1755 500 CONTINUE 1756 600 CONTINUE 1757C 1758 IF (IPRINT .GT. 50) THEN 1759 CALL AROUND('ccsl_tgb : T^gB_ao matrix: ') 1760 CALL CC_PRFCKAO(TGB,ISYMTR) 1761 ENDIF 1762C 1763 CALL GPCLOSE(LUCSOL,'KEEP') 1764C 1765 IF (DISCEX) THEN 1766 CALL WCLOSE2(LUMMET, FILMME, 'KEEP') 1767 CALL WCLOSE2(LUMMXI, FILMMX, 'KEEP') 1768 END IF 1769C 1770 END 1771******************************************************************* 1772**************************************************************** 1773C /* Deck ccsl_gtr */ 1774 SUBROUTINE CCSL_GTR(RHO1,RHO2,ISYRES, 1775 * LISTA,IDLSTA,ISYMTA, 1776 * LISTB,IDLSTB,ISYMTB, 1777 * LISTC,IDLSTC,ISYMTC, 1778 * MODEL,WORK,LWORK) 1779C 1780C----------------------------------------------------------------- 1781C JK+OC, feb.03 1782C Calculates the contribution to the CC/dielectric G matrx 1783C transformations. 1784C 1785C G_my,ny,sigma = 1/2 P(my,ny,sigma) 1786C <Lambda|[[T^g,mu,tau_ny],tau_sigma]|CC > 1787C 1788C The permutation creates 3 different contributions. 1789C 1790C------------------------------------------------------------------ 1791#include "implicit.h" 1792#include "maxorb.h" 1793#include "mxcent.h" 1794 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, IZERO = 0 , TWO = 2.0D0) 1795 PARAMETER (HALF = 0.5D0) 1796#include "dummy.h" 1797#include "iratdef.h" 1798#include "priunit.h" 1799#include "ccorb.h" 1800#include "ccsdsym.h" 1801#include "ccsdio.h" 1802#include "ccsdinp.h" 1803#include "ccfield.h" 1804#include "exeinf.h" 1805#include "ccfdgeo.h" 1806#include "ccslvinf.h" 1807#include "ccinftap.h" 1808C 1809 DIMENSION WORK(LWORK),RHO1(*),RHO2(*) 1810C 1811 CHARACTER*(3) LISTA, LISTB, LISTC 1812 CHARACTER*8 LABEL 1813 CHARACTER*10 MODEL 1814 LOGICAL LEXIST, LSAME 1815 INTEGER ISYCB,ISYRES,ISYMTA,ISYMTB,ISYMTC 1816 INTEGER IDLSTA,IDLSTB,IDLSTC 1817C 1818 INTEGER KDUM 1819 PARAMETER( KDUM = +99 999 999 ) ! dummy address for work space 1820C 1821 IF (IPRINT.GT.10) THEN 1822 WRITE(LUPRI,*)'CCSL_GTR: MODEL =',MODEL(1:4) 1823 WRITE(LUPRI,*)'CCSL_GTR: ISYRES =',ISYRES 1824 WRITE(LUPRI,*)'CCSL_GTR: ISYMTA =',ISYMTA 1825 WRITE(LUPRI,*)'CCSL_GTR: ISYMTB =',ISYMTB 1826 WRITE(LUPRI,*)'CCSL_GTR: ISYMTC =',ISYMTC 1827 WRITE(LUPRI,*)'CCSL_GTR: LISTA =',LISTA 1828 WRITE(LUPRI,*)'CCSL_GTR: LISTB =',LISTB 1829 WRITE(LUPRI,*)'CCSL_GTR: LISTC =',LISTC 1830 WRITE(LUPRI,*)'CCSL_GTR: IDLSTA =',IDLSTA 1831 WRITE(LUPRI,*)'CCSL_GTR: IDLSTB =',IDLSTB 1832 WRITE(LUPRI,*)'CCSL_GTR: IDLSTC =',IDLSTC 1833 CALL FLSHFO(LUPRI) 1834 ENDIF 1835C 1836C Symmetry: 1837C --------- 1838C 1839 ISYCB = MULD2H(ISYMTB,ISYMTC) 1840C 1841 IF (ISYCB .NE. ISYRES) THEN 1842 CALL QUIT( 'Symmetry problem in CCSL_GTR 1') 1843 END IF 1844C 1845C 1846C First we concentrate on the effective operator T^gB 1847C ----------------------------------------------------- 1848C 1849 KT1AMPA = 1 1850 KTGB = KT1AMPA + NT1AM(ISYMTB) 1851 KWRK1 = KTGB + N2BST(ISYMTB) 1852 LWRK1 = LWORK - KWRK1 1853C 1854 IF (LWRK1 .LT. 0) THEN 1855 CALL QUIT('Insuff. work in CCSL_GTR 1') 1856 END IF 1857C 1858 CALL DZERO(WORK(KT1AMPA),NT1AM(ISYMTB)) 1859 CALL DZERO(WORK(KTGB),N2BST(ISYMTB)) 1860C 1861C Read one electron excitation part of response vector 1862C 1863 IOPT = 1 1864 CALL CC_RDRSP(LISTB,IDLSTB,ISYMTB,IOPT,MODEL, 1865 * WORK(KT1AMPA),WORK(KDUM)) 1866C 1867 CALL CCSL_TGB(WORK(KT1AMPA),ISYMTB,LISTB,IDLSTB,WORK(KTGB),'ET', 1868 * MODEL,WORK(KWRK1),LWRK1) 1869C 1870 LABEL = 'GIVE INT' 1871 CALL CCLR_FA(LABEL,ISYMTB,LISTC,IDLSTC, 1872 & LISTA,IDLSTA,WORK(KTGB),WORK(KWRK1),LWRK1) 1873C 1874C 1875 LSAME = (LISTC.EQ.LISTB .AND. IDLSTC.EQ.IDLSTB) 1876 IF (LSAME) THEN 1877 FACTSLV = TWO 1878 ELSE 1879 FACTSLV = ONE 1880 END IF 1881C 1882 CALL DAXPY(NT1AM(ISYCB),FACTSLV,WORK(KWRK1),1,RHO1,1) 1883 CALL DAXPY(NT2AM(ISYCB),FACTSLV, 1884 * WORK(KWRK1+NT1AM(ISYCB)),1,RHO2,1) 1885C 1886 IF (.NOT. (LSAME)) THEN 1887C 1888C second we consider the effective operator T^gC 1889C ----------------------------------------------------- 1890C 1891 KT1AMPA = 1 1892 KTGC = KT1AMPA + NT1AM(ISYMTC) 1893 KWRK1 = KTGC + N2BST(ISYMTC) 1894 LWRK1 = LWORK - KWRK1 1895C 1896 IF (LWRK1 .LT. 0) CALL QUIT('Insuff. work in CCSL_GTR 2') 1897C 1898 CALL DZERO(WORK(KT1AMPA),NT1AM(ISYMTC)) 1899 CALL DZERO(WORK(KTGC),N2BST(ISYMTC)) 1900C 1901C Read one electron excitation part of response vector 1902C 1903 IOPT = 1 1904 CALL CC_RDRSP(LISTC,IDLSTC,ISYMTC,IOPT,MODEL, 1905 * WORK(KT1AMPA),WORK(KDUM)) 1906C 1907 CALL CCSL_TGB(WORK(KT1AMPA),ISYMTC,LISTC,IDLSTC,WORK(KTGC),'ET', 1908 * MODEL,WORK(KWRK1),LWRK1) 1909C 1910 LABEL = 'GIVE INT' 1911 CALL CCLR_FA(LABEL,ISYMTC,LISTB,IDLSTB, 1912 & LISTA,IDLSTA,WORK(KTGC),WORK(KWRK1),LWRK1) 1913C 1914C 1915 CALL DAXPY(NT1AM(ISYCB),ONE,WORK(KWRK1),1,RHO1,1) 1916 CALL DAXPY(NT2AM(ISYCB),ONE, 1917 * WORK(KWRK1+NT1AM(ISYCB)),1,RHO2,1) 1918C 1919 END IF 1920C 1921C finally, we consider the effective operator T^g sigma 1922C ----------------------------------------------------- 1923C 1924 KTGCB = 1 1925 KWRK1 = KTGCB + N2BST(ISYCB) 1926 LWRK1 = LWORK - KWRK1 1927 IF (LWRK1.LT.0) CALL QUIT( 'Too little work in CCSL_GTR, 3') 1928C 1929 CALL DZERO(WORK(KTGCB),N2BST(ISYCB)) 1930C 1931 CALL CCSL_TGCB(ISYRES,LISTA,LISTB,LISTC, 1932 * IDLSTA,IDLSTB,IDLSTC, 1933 * ISYMTA,ISYMTB,ISYMTC,WORK(KTGCB), 1934 * MODEL,WORK(KWRK1),LWRK1) 1935C 1936 NAMPF = NT1AM(ISYCB) + NT2AM(ISYCB) 1937C 1938 KETA = KWRK1 1939 KWRK2 = KETA + NAMPF 1940 LWRK2 = LWORK - KWRK2 1941 IF (LWRK2.LT.0) CALL QUIT( 'Too little work in CCSL_GTR, 4') 1942C 1943 CALL DZERO(WORK(KETA),NAMPF) 1944C 1945 LABEL = 'GIVE INT' 1946 CALL CC_ETAC(ISYCB,LABEL,WORK(KETA), 1947 * LISTA,IDLSTA,0,WORK(KTGCB),WORK(KWRK2),LWRK2) 1948C 1949 KETA1 = KETA 1950 KETA2 = KETA1 + NT1AM(ISYCB) 1951C 1952 CALL DAXPY(NT1AM(ISYCB),ONE,WORK(KETA1),1,RHO1,1) 1953 CALL DAXPY(NT2AM(ISYCB),ONE,WORK(KETA2),1,RHO2,1) 1954C 1955 END 1956************************************************************** 1957C /* Deck ccsl_tgcb */ 1958 SUBROUTINE CCSL_TGCB(ISYRES,LISTA,LISTB,LISTC, 1959 * IDLSTA,IDLSTB,IDLSTC, 1960 * ISYMTA,ISYMTB,ISYMTC,TGCB, 1961 * MODEL,WORK,LWORK) 1962C 1963C Constructs effective operator equal to 1964C T^gCB = -2 sum_lm g_l <Lambda|[[T_lm,C],B]|CC>T_lm 1965C 1966C JK+OC, 03 1967C----------------------------------------------------------------------------- 1968C 1969#include "implicit.h" 1970#include "maxorb.h" 1971#include "mxcent.h" 1972 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, IZERO = 0 , TWO = 2.0D0) 1973 PARAMETER (HALF = 0.5D0) 1974#include "dummy.h" 1975#include "iratdef.h" 1976#include "priunit.h" 1977#include "ccorb.h" 1978#include "ccsdsym.h" 1979#include "ccsdio.h" 1980#include "ccsdinp.h" 1981#include "ccfield.h" 1982#include "exeinf.h" 1983#include "ccfdgeo.h" 1984#include "ccslvinf.h" 1985#include "ccinftap.h" 1986C 1987 DIMENSION WORK(LWORK) 1988 DIMENSION TGCB(*) 1989C 1990 INTEGER KDUM 1991 INTEGER IDLSTA,IDLSTB,IDLSTC,ISYMTC,ISYMTA,ISYMTB,ISYRES 1992 PARAMETER( KDUM = +99 999 999 ) ! dummy address for work space 1993C 1994 CHARACTER*(*) LISTA,LISTB,LISTC 1995 CHARACTER*10 MODEL 1996C 1997 CHARACTER*8 LABEL 1998 LOGICAL LEXIST 1999C 2000 DOUBLE PRECISION GL 2001 EXTERNAL GL 2002C 2003 IF (IPRINT.GT.10) THEN 2004 WRITE(LUPRI,*)'CCSL_TGCB: ISYRES =',ISYRES 2005 WRITE(LUPRI,*)'CCSL_TGCB: ISYTMA =',ISYMTA 2006 WRITE(LUPRI,*)'CCSL_TGCB: ISYTMB =',ISYMTB 2007 WRITE(LUPRI,*)'CCSL_TGCB: ISYTMC =',ISYMTC 2008 WRITE(LUPRI,*)'CCSL_TGCB: LISTA =',LISTA 2009 WRITE(LUPRI,*)'CCSL_TGCB: LISTB =',LISTB 2010 WRITE(LUPRI,*)'CCSL_TGCB: LISTC =',LISTC 2011 WRITE(LUPRI,*)'CCSL_TGCB: IDLSTA =',IDLSTA 2012 WRITE(LUPRI,*)'CCSL_TGCB: IDLSTB =',IDLSTB 2013 WRITE(LUPRI,*)'CCSL_TGCB: IDLSTC =',IDLSTC 2014 CALL FLSHFO(LUPRI) 2015 ENDIF 2016C-------------------------------------- 2017C Readin trial vector (B) from file 2018C-------------------------------------- 2019C 2020 KT1AMB = 1 2021 KT2AMB = KT1AMB + NT1AM(ISYMTB) 2022 KWRK1 = KT2AMB + NT2AM(ISYMTB) 2023 LWRK1 = LWORK - KWRK1 2024C 2025 IF (LWRK1 .LT. 0) THEN 2026 CALL QUIT('Insuff. work in CCSL_TGCB 1') 2027 END IF 2028C 2029 IOPT = 3 2030 CALL CC_RDRSP(LISTB,IDLSTB,ISYMTB,IOPT,MODEL, 2031 * WORK(KT1AMB),WORK(KT2AMB)) 2032C 2033C------------------ 2034C Symmetry 1 2035C------------------ 2036C 2037 ISYCB = MULD2H(ISYMTB,ISYMTC) 2038C 2039 IF (ISYCB .NE. ISYRES) THEN 2040 CALL QUIT( 'Symmetry problem in CCSL_TGCB') 2041 END IF 2042C 2043C--------------------------------- 2044C Read in integrals from file. 2045C--------------------------------- 2046C 2047 LUCSOL = -1 2048 CALL GPOPEN(LUCSOL,'AOSOLINT','OLD',' ','UNFORMATTED',IDUMMY, 2049 & .FALSE.) 2050 REWIND(LUCSOL) 2051 CALL MOLLAB('SOLVRLM ',LUCSOL,LUERR) 2052 READ (LUCSOL) 2053 CALL READT(LUCSOL,NLMCU,WORK(KWRK1)) 2054C 2055 TEMP = 0.00 2056 LM = 0 2057 DO 600 L = 0,LMAXCU 2058 READ (LUCSOL) L1,(ISYTLM(M),M=1,2*L+1) 2059 IF (L1 .NE. L) THEN 2060 WRITE (LUERR,*) 2061 * 'ERROR CCSL_TGCB: L from LUCSOL not as expected' 2062 WRITE (LUERR,*) 'L from 600 loop:',L 2063 WRITE (LUERR,*) 'L from LUCSOL :',L1 2064 CALL QUIT('ERROR CCSL_TGCB: L from LUCSOL not as expected') 2065 END IF 2066C 2067 DO 500 M = -L,L 2068 LM = LM + 1 2069 IF (IPRINT .GE. 15) THEN 2070 WRITE(LUPRI,'(/A,2I5/A)') ' --- l, m :',L,M, 2071 * ' ====================' 2072 WRITE(LUPRI,'(A,I2)') ' Symmetry :',ABS(ISYTLM(L+M+1)) 2073 END IF 2074C 2075C------------------ 2076C Symmetry 2 2077C------------------ 2078C 2079 IF (ISYTLM(L+M+1) .NE. ISYCB) THEN 2080 READ (LUCSOL) 2081 GO TO 500 2082 ENDIF 2083C 2084C-------------------------------- 2085C Read T(l,m) in AO basis. 2086C-------------------------------- 2087C 2088 KTLMAO = KWRK1 2089 KWRK2 = KTLMAO + N2BST(ISYCB) 2090 LWRK2 = LWORK - KWRK2 2091 IF (LWRK2.LT.0) CALL QUIT( 'Too little work in CCSL_TGCB, 2') 2092C 2093 CALL CC_SLVINT(WORK(KTLMAO),WORK(KWRK2),LWRK2,ISYCB) 2094C 2095 IF (IPRINT .GT. 25) THEN 2096 CALL AROUND('CCSL_TGCB : Tlm_ao matrix: - cc storage') 2097 CALL CC_PRFCKAO(WORK(KTLMAO),ISYCB) 2098 ENDIF 2099C 2100C------------------------------------------------------------------- 2101C 2102 LABEL = 'GIVE INT' 2103 CALL CCLR_FA(LABEL,ISYCB,LISTC,IDLSTC, 2104 & LISTA,IDLSTA,WORK(KTLMAO),WORK(KWRK2),LWRK2) 2105C 2106C 2107 CBDOT1 = DDOT(NT1AM(ISYMTB),WORK(KWRK2),1,WORK(KT1AMB),1) 2108 CBDOT2 = DDOT(NT2AM(ISYMTB),WORK(KWRK2+NT1AM(ISYMTB)),1, 2109 * WORK(KT2AMB),1) 2110C 2111 CBDOT = CBDOT1 + CBDOT2 2112C 2113C------------------------- 2114C Find G_lm factor 2115C------------------------- 2116C 2117 GLFAC = GL(L,EPOPCU,RCAVCU) 2118 FACT = 2.0D0*GLFAC*CBDOT 2119C 2120C-------------------------------------------------------------- 2121C Put factor on integrals and (for a given l,m) the 2122C effective operator has now been constructed. 2123C-------------------------------------------------------------- 2124C 2125 CALL DAXPY(N2BST(ISYCB),FACT,WORK(KTLMAO),1, 2126 * TGCB,1) 2127C--------------------------------------------------------------- 2128 2129 500 CONTINUE 2130 600 CONTINUE 2131C 2132 IF (IPRINT .GT. 50) THEN 2133 CALL AROUND('ccsl_tgb : T^gCB_ao matrix: ') 2134 CALL CC_PRFCKAO(TGCB,ISYCB) 2135 ENDIF 2136C 2137 CALL GPCLOSE(LUCSOL,'KEEP') 2138C 2139 END 2140********************************************************************** 2141C /* Deck ccsl_btr */ 2142 SUBROUTINE CCSL_BTR(RHO1,RHO2,ISYMAB, 2143 * LISTA,IDLSTA,ISYMA, 2144 * LISTB,IDLSTB,ISYMB, 2145 * MODEL,WORK,LWORK) 2146C 2147C JK+OC, 03 2148C--------------------------------------------------------------- 2149C 2150#include "implicit.h" 2151#include "maxorb.h" 2152#include "mxcent.h" 2153 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, IZERO = 0 , TWO = 2.0D0) 2154 PARAMETER (HALF = 0.5D0) 2155#include "dummy.h" 2156#include "iratdef.h" 2157#include "priunit.h" 2158#include "ccorb.h" 2159#include "ccsdsym.h" 2160#include "ccsdio.h" 2161#include "ccsdinp.h" 2162#include "ccfield.h" 2163#include "exeinf.h" 2164#include "ccfdgeo.h" 2165#include "ccslvinf.h" 2166#include "ccinftap.h" 2167C 2168 DIMENSION WORK(LWORK),RHO1(*),RHO2(*) 2169C 2170 CHARACTER*(*) LISTA,LISTB 2171 CHARACTER*(3) LISTC 2172 CHARACTER*8 LABEL 2173 CHARACTER*10 MODEL 2174 LOGICAL LEXIST, LSAME 2175 INTEGER IDLSTA, IDLSTB, ISYMAB, ISYMA, ISYMB 2176 PARAMETER( KDUM = +99 999 999 ) ! dummy address for work space 2177C 2178 IF (IPRINT.GT.10) THEN 2179 WRITE(LUPRI,*)'CCSL_BTR: ISYMAB =',ISYMAB 2180 WRITE(LUPRI,*)'CCSL_BTR: ISYMA =',ISYMA 2181 WRITE(LUPRI,*)'CCSL_BTR: ISYMB =',ISYMB 2182 WRITE(LUPRI,*)'CCSL_BTR: LISTA =',LISTA 2183 WRITE(LUPRI,*)'CCSL_BTR: LISTB =',LISTB 2184 WRITE(LUPRI,*)'CCSL_BTR: IDLSTA =',IDLSTA 2185 WRITE(LUPRI,*)'CCSL_BTR: IDLSTB =',IDLSTB 2186 CALL FLSHFO(LUPRI) 2187 ENDIF 2188C 2189C First we concentrate on the effective operator T^gA 2190C ----------------------------------------------------- 2191C 2192 KT1AMPA = 1 2193 KTGA = KT1AMPA + NT1AM(ISYMA) 2194 KWRK1 = KTGA + N2BST(ISYMA) 2195 LWRK1 = LWORK - KWRK1 2196C 2197 IF (LWRK1 .LT. 0) THEN 2198 CALL QUIT('Insuff. work in CCSL_BTR 1') 2199 END IF 2200C 2201 CALL DZERO(WORK(KT1AMPA),NT1AM(ISYMA)) 2202 CALL DZERO(WORK(KTGA),N2BST(ISYMA)) 2203C 2204C Read one electron excitation part of response vector 2205C 2206 IOPT = 1 2207 CALL CC_RDRSP(LISTA,IDLSTA,ISYMA,IOPT,MODEL, 2208 * WORK(KT1AMPA),WORK(KDUM)) 2209C 2210 CALL CCSL_TGB(WORK(KT1AMPA),ISYMA,LISTA,IDLSTA,WORK(KTGA),'ET', 2211 * MODEL,WORK(KWRK1),LWRK1) 2212C 2213 LABEL = 'GIVE INT' 2214 CALL CCCR_AA(LABEL,ISYMA,LISTB,IDLSTB, 2215 & WORK(KTGA),WORK(KWRK1),LWRK1) 2216 CALL CCLR_DIASCL(WORK(KWRK1+NT1AM(ISYMAB)),2.0d0,ISYMAB) 2217C 2218C 2219 LSAME = (LISTA.EQ.LISTB .AND. IDLSTA.EQ.IDLSTB) 2220 IF (LSAME) THEN 2221 FACTSLV = TWO 2222 ELSE 2223 FACTSLV = ONE 2224 END IF 2225C 2226 CALL DAXPY(NT1AM(ISYMAB),FACTSLV,WORK(KWRK1),1,RHO1,1) 2227 CALL DAXPY(NT2AM(ISYMAB),FACTSLV, 2228 * WORK(KWRK1+NT1AM(ISYMAB)),1,RHO2,1) 2229C 2230C 2231 IF (.NOT.(LSAME)) THEN 2232C 2233C Next we concentrate on the effective operator T^gB 2234C ----------------------------------------------------- 2235C 2236 KT1AMPB = 1 2237 KTGB = KT1AMPB + NT1AM(ISYMB) 2238 KWRK1 = KTGB + N2BST(ISYMB) 2239 LWRK1 = LWORK - KWRK1 2240C 2241 IF (LWRK1 .LT. 0) THEN 2242 CALL QUIT('Insuff. work in CCSL_BTR 2') 2243 END IF 2244C 2245 CALL DZERO(WORK(KT1AMPB),NT1AM(ISYMB)) 2246 CALL DZERO(WORK(KTGB),N2BST(ISYMB)) 2247C 2248C Read one electron excitation part of response vector 2249C 2250 IOPT = 1 2251 CALL CC_RDRSP(LISTB,IDLSTB,ISYMB,IOPT,MODEL, 2252 * WORK(KT1AMPB),WORK(KDUM)) 2253 2254 CALL CCSL_TGB(WORK(KT1AMPB),ISYMB,LISTB,IDLSTB,WORK(KTGB),'ET', 2255 * MODEL,WORK(KWRK1),LWRK1) 2256C 2257 LABEL = 'GIVE INT' 2258 CALL CCCR_AA(LABEL,ISYMB,LISTA,IDLSTA, 2259 & WORK(KTGB),WORK(KWRK1),LWRK1) 2260 CALL CCLR_DIASCL(WORK(KWRK1+NT1AM(ISYMAB)),2.0d0,ISYMAB) 2261C 2262C 2263 CALL DAXPY(NT1AM(ISYMAB),ONE,WORK(KWRK1),1,RHO1,1) 2264 CALL DAXPY(NT2AM(ISYMAB),ONE, 2265 * WORK(KWRK1+NT1AM(ISYMAB)),1,RHO2,1) 2266C 2267 END IF 2268C 2269C Finally, we concentrate on the effective operator T^gAB 2270C ------------------------------------------------------- 2271C 2272 KTGAB = 1 2273 KWRK1 = KTGAB + N2BST(ISYMAB) 2274 LWRK1 = LWORK - KWRK1 2275 IF (LWRK1.LT.0) CALL QUIT( 'Too little work in CCSL_BTR, 3') 2276C 2277 CALL DZERO(WORK(KTGAB),N2BST(ISYMAB)) 2278C 2279 LISTC = 'L0' 2280 IDLSTC = 1 2281 ISYMC = ILSTSYM(LISTC,IDLSTC) 2282C 2283 CALL CCSL_TGCB(ISYMAB,LISTC,LISTA,LISTB, 2284 * IDLSTC,IDLSTA,IDLSTB, 2285 * ISYMC,ISYMA,ISYMB,WORK(KTGAB), 2286 * MODEL,WORK(KWRK1),LWRK1) 2287C 2288 NAMPF = NT1AM(ISYMAB) + NT2AM(ISYMAB) 2289C 2290 KXI = KWRK1 2291 KWRK2 = KXI + NAMPF 2292 LWRK2 = LWORK - KWRK2 2293 IF (LWRK2.LT.0) CALL QUIT( 'Too little work in CCSL_BTR, 4') 2294C 2295 CALL DZERO(WORK(KXI),NAMPF) 2296C 2297 LABEL = 'GIVE INT' 2298 CALL CC_XKSI(WORK(KXI),LABEL,ISYMAB,0,WORK(KTGAB), 2299 * WORK(KWRK2),LWRK2) 2300C 2301 KXI1 = KXI 2302 KXI2 = KXI1 + NT1AM(ISYMAB) 2303C 2304 CALL CCLR_DIASCL(WORK(KXI2),2.0d0,ISYMAB) 2305C 2306 CALL DAXPY(NT1AM(ISYMAB),ONE,WORK(KXI1),1,RHO1,1) 2307 CALL DAXPY(NT2AM(ISYMAB),ONE,WORK(KXI2),1,RHO2,1) 2308 2309 END 2310******************************************************************** 2311C /* Deck ccsl_xieta */ 2312 SUBROUTINE CCSL_XIETA(WORK,LWORK) 2313C 2314C----------------------------------------------------------------- 2315C JK, april .03 2316C 2317C Calculates the CCSL ETA and XI vectors and write them to files. 2318C 2319C----------------------------------------------------------------------------- 2320C 2321#include "implicit.h" 2322#include "maxorb.h" 2323#include "mxcent.h" 2324 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, IZERO = 0 , TWO = 2.0D0) 2325 PARAMETER (HALF = 0.5D0) 2326#include "dummy.h" 2327#include "iratdef.h" 2328#include "priunit.h" 2329#include "ccorb.h" 2330#include "ccsdsym.h" 2331#include "ccsdio.h" 2332#include "ccsdinp.h" 2333#include "ccfield.h" 2334#include "exeinf.h" 2335#include "ccfdgeo.h" 2336#include "ccslvinf.h" 2337#include "ccinftap.h" 2338C 2339 DIMENSION WORK(LWORK) 2340C 2341 INTEGER LUMMET, LUMMXI, ISYMTR 2342 INTEGER KDUM 2343 PARAMETER( KDUM = +99 999 999 ) ! dummy address for work space 2344C 2345 CHARACTER*(8) FILMME, FILMMX, LABEL 2346 CHARACTER*2 LIST 2347C 2348 LUMMET = -1 2349 LUMMXI = -1 2350 FILMME = 'CCSL_ETA' 2351 FILMMX = 'CCSL__XI' 2352C 2353 CALL WOPEN2(LUMMET, FILMME, 64, 0) 2354 CALL WOPEN2(LUMMXI, FILMMX, 64, 0) 2355 2356C--------------------------------- 2357C Read in integrals from file. 2358C--------------------------------- 2359C 2360 LUCSOL = -1 2361 CALL GPOPEN(LUCSOL,'AOSOLINT','OLD',' ','UNFORMATTED',IDUMMY, 2362 & .FALSE.) 2363 REWIND(LUCSOL) 2364 CALL MOLLAB('SOLVRLM ',LUCSOL,LUERR) 2365 READ (LUCSOL) 2366 CALL READT(LUCSOL,NLMCU,WORK(1)) 2367C 2368 TEMP = 0.00 2369 LM = 0 2370 IADR1 = 1 2371 IADR2 = 1 2372C 2373 DO 600 L = 0,LMAXCU 2374 READ (LUCSOL) L1,(ISYTLM(M),M=1,2*L+1) 2375 IF (L1 .NE. L) THEN 2376 WRITE (LUERR,*) 2377 * 'ERROR CCSL_XIETA: L from LUCSOL not as expected' 2378 WRITE (LUERR,*) 'L from 600 loop:',L 2379 WRITE (LUERR,*) 'L from LUCSOL :',L1 2380 CALL QUIT('ERROR CCSL_XIETA: L from LUCSOL not as expected') 2381 END IF 2382C 2383 DO 500 M = -L,L 2384 LM = LM + 1 2385C IF (IPRINT .GE. 15) THEN 2386 WRITE(LUPRI,'(/A,2I5/A)') ' --- l, m :',L,M, 2387 * ' ====================' 2388 WRITE(LUPRI,'(A,I2)') ' Symmetry :',ABS(ISYTLM(L+M+1)) 2389C END IF 2390C 2391C 2392C Symmetry of the T_lm operator: 2393C ------------------------------ 2394C 2395 ISYMTR = ABS(ISYTLM(L+M+1)) 2396C 2397C Number of amplitudes in the given symmetry: 2398C ------------------------------------------- 2399C 2400 NTAMP1 = NT1AM(ISYMTR) 2401 NTAMP2 = NT2AM(ISYMTR) 2402 IF (CCS) NTAMP2 = 0 2403 NTAMP = NTAMP1 + NTAMP2 2404C 2405 KTLMAO = 1 2406 KWRK1 = KTLMAO + N2BST(ISYMTR) 2407 LWRK1 = LWORK - KWRK1 2408 IF (LWRK1.LT.0) CALL QUIT( 'Too little work in CCSL_XIETA, 1') 2409 CALL DZERO(WORK(KTLMAO),N2BST(ISYMTR)) 2410C 2411 CALL CC_SLVINT(WORK(KTLMAO),WORK(KWRK1),LWRK1,ISYMTR) 2412C 2413 IF (IPRINT .GT. 25) THEN 2414 CALL AROUND('CCSL_XIETA : Tlm_ao matrix: - cc storage') 2415 CALL CC_PRFCKAO(WORK(KTLMAO),ISYMTR) 2416 ENDIF 2417C 2418C Allocate space for the result vector: 2419C ------------------------------------- 2420C 2421 KXI = KWRK1 2422 KWRK2 = KXI + NTAMP 2423 LWRK2 = LWORK - KWRK2 2424 IF (LWRK2.LT.0) CALL QUIT( 'Too little work in CCSL_XIETA, 2') 2425 CALL DZERO(WORK(KXI),NTAMP) 2426C 2427C Perform the calculation 2428C ----------------------- 2429C 2430 LABEL = 'GIVE INT' 2431 LIST = 'L0' 2432 IDLINO = 1 2433C 2434 CALL CC_ETAC(ISYMTR,LABEL,WORK(KXI), 2435 * LIST,IDLINO,0,WORK(KTLMAO),WORK(KWRK2),LWRK2) 2436 CALL PUTWA2(LUMMET,FILMME,WORK(KXI),IADR1,NTAMP) 2437 IADR1 = IADR1 + NTAMP 2438C 2439 CALL CC_XKSI(WORK(KXI),LABEL,ISYMTR,0,WORK(KTLMAO), 2440 * WORK(KWRK2),LWRK2) 2441 CALL PUTWA2(LUMMXI,FILMMX,WORK(KXI),IADR2,NTAMP) 2442 IADR2 = IADR2 + NTAMP 2443C 2444 500 CONTINUE 2445 600 CONTINUE 2446C 2447 CALL GPCLOSE(LUCSOL,'KEEP') 2448C 2449 CALL WCLOSE2(LUMMET, FILMME, 'KEEP') 2450 CALL WCLOSE2(LUMMXI, FILMMX, 'KEEP') 2451C 2452 END 2453C*********************************************************************** 2454