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_qm3 */ 20 SUBROUTINE CC_QM3(AODEN,CONVERGED,WORK,LWORK) 21C********************************************************************* 22C 23C Purpose: Calculation of CC/MM wave-function and energy. 24C CCS(CIS/HF)(nci), MP2(nci), CC2(nci), CCSD, CC3(nci), MCC2(nci) 25C QM3, JK+AO+OC+KMI,01 26C 27C********************************************************************* 28C 29#include "implicit.h" 30#include "maxorb.h" 31#include "mxcent.h" 32#include "dummy.h" 33#include "iratdef.h" 34#include "priunit.h" 35#include "ccorb.h" 36#include "ccsdsym.h" 37#include "ccsdio.h" 38#include "ccsdinp.h" 39#include "ccfield.h" 40#include "exeinf.h" 41#include "ccfdgeo.h" 42#include "ccslvinf.h" 43#include "ccinftap.h" 44#include "qm3.h" 45C 46 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, IZERO = 0 , TWO = 2.0D0) 47C 48 LOGICAL CONVERGED 49 DIMENSION WORK(LWORK),AODEN(*),FFJ(3) 50 CHARACTER MODEL*10 51 CHARACTER MODELPRI*4 52 CHARACTER LABEL*8 53C 54#if defined (SYS_CRAY) 55 REAL DEMMMM 56#else 57 DOUBLE PRECISION DEMMMM 58#endif 59C 60C-------------------------------------------------- 61C Header for the CCMM output in each iteration: 62C-------------------------------------------------- 63C 64 WRITE (LUPRI,'(1X,A,/)') ' ' 65 WRITE (LUPRI,'(1X,A)') 66 *'*****************************************************'// 67 *'************' 68 WRITE (LUPRI,'(1X,A)') 69 *'**** Output from Coupled Cluster/Molecular Mechanics program'// 70 *' ****' 71 WRITE (LUPRI,'(1X,A)') 72 *'*****************************************************'// 73 *'************' 74C 75C----------------------------------------------- 76C Testing if the CC model is CC2 or CCSD: 77C----------------------------------------------- 78C 79 MODEL = 'CCSD' 80C 81 IF (CCSD) THEN 82 CALL AROUND('Coupled Cluster model is: CCSD') 83 MODEL = 'CCSD' 84 MODELPRI = 'CCSD' 85 ENDIF 86C 87 IF (CC2.AND.(.NOT.MCC2)) THEN 88 CALL AROUND('Coupled Cluster model is: CC2') 89 MODEL = 'CC2' 90 MODELPRI = ' CC2' 91 ENDIF 92C 93 IF (.NOT. (CCSD.OR.CC2)) THEN 94 CALL QUIT('CCMM only implemented for the ' // 95 & 'CC2 and CCSD parameterizations.') 96 ENDIF 97C 98 IF (IQM3PR .GT .10) 99 &WRITE(LUPRI,*) 'CC_MM-1: Workspace:',LWORK 100C 101C----------------------------------------------- 102C Calculate the E(QM/MM) part of the energy: 103C----------------------------------------------- 104C 105 ECMCON = ZERO 106 107 IF (.NOT. NYQMMM) THEN 108 CALL CC_QM3E(ECMCON,EPOL,EEL,EELEL,EREP,AODEN, 109 & WORK,LWORK) 110 ELSE IF (NYQMMM) THEN 111 CALL CC_QMMME(ECMCON,EPOL,EEL,EELEL,EREP,AODEN, 112 & WORK,LWORK) 113 END IF 114C 115C---------------------------- 116C Calculate new E(QM/MM): 117C---------------------------- 118C 119 IF (.NOT. (OLDTG) .AND. .NOT. (NYQMMM)) THEN 120 ECMCU = (ECCGRS - EELEL) + ECMCON 121 ECCGRS1 = (ECCGRS - EELEL) 122 ELSE 123 ECMCU = ECCGRS + ECMCON 124 ECCGRS1 = ECCGRS 125 END IF 126C 127 IF (ABS(ECMCU-ECCPR).LT.CVGESOL) LSLECVG = .TRUE. 128C 129 WRITE(LUPRI,'(//1X,A,I3,A,F20.10)') 130 * 'E(QM/MM) contribution in iteration',ICCSLIT,': ', 131 * ECMCON 132C 133 WRITE(LUPRI,'(1X,A,F20.10)') 134 * 'CC energy in the current CCMM iteration: ',ECMCU 135C 136 WRITE(LUPRI,'(1X,A,F20.10)') 137 * 'CC energy in the previous CCMM iteration: ',ECCPR 138C 139 WRITE(LUPRI,'(1X,A,E20.6//)') 140 * 'Change in Total energy in this CCMM it.: ', 141 * ECMCU - ECCPR 142C 143 ECCPR = ECMCU 144C 145 CONVERGED = .FALSE. 146 IF (LSLECVG.AND.LSLTCVG.AND.LSLLCVG) THEN 147 CONVERGED = .TRUE. 148 IF (LOEC3) THEN 149C 150 EPOL = ZERO 151C 152 WRITE(LUPRI,'(1X,A)') 153 & 'Perturbative correction to polarization energy is calculated' 154 & ,'(Model SPC_EC3)' 155C 156 CALL EPERT2(EPOL,WORK,LWORK) 157 END IF 158C 159C Here we see if any static fields are to be added 160C 161 FFJ(1) = 0.0D0 162 FFJ(2) = 0.0D0 163 FFJ(3) = 0.0D0 164C 165 IF (RELMOM) THEN 166 DO 330 IF =1, NFIELD 167 IF (LFIELD(IF) .EQ. 'XDIPLEN') FFJ(1) = FFJ(1) + EFIELD(IF) 168 IF (LFIELD(IF) .EQ. 'YDIPLEN') FFJ(2) = FFJ(2) + EFIELD(IF) 169 IF (LFIELD(IF) .EQ. 'ZDIPLEN') FFJ(3) = FFJ(3) + EFIELD(IF) 170 330 CONTINUE 171 END IF 172C 173 CALL QM3_EMMMM(DEMMMM,TEMPOL,FFJ) 174C 175 IF (.NOT.(RELMOM)) PEDIP1 = 0.0D0 176C 177 EMMPOL = TEMPOL 178 IF (.NOT.(RELMOM)) THEN 179 EMM_MM = EMMELC + EMMVDW + EMMPOL 180 ELSE 181 EMM_MM = EMMELC + EMMVDW + EMMPOL + PEDIP1 182 END IF 183C 184C--------------------------------------------------- 185C Final output for the QM/MM results with the 186C converged wave function: 187C--------------------------------------------------- 188C 189C First the MM/MM energy with the QM induced 190C dipole moments: 191C------------------------------------------------- 192C 193 IF (.NOT. NYQMMM) THEN 194 WRITE(LUPRI,'(//1X,72A)') ('-',J=1,72) 195 WRITE(LUPRI,'(1X,A56,A16)') 196 * '| --- The MM/MM classical interaction energy ', 197 * '--- |' 198 WRITE(LUPRI,'(1X,72A)') ('-',J=1,72) 199 WRITE(LUPRI,'(1X,A1,A50,A7,F12.8,A2)') '|', 200 * ' Eelec = Sum_n,s[ (Q_n*Q_s)/|R_n - R_s| ] ', 201 * '| ',EMMELC,' |' 202 IF (RELMOM) THEN 203 WRITE(LUPRI,'(1X,72A)') ('-',J=1,72) 204 WRITE(LUPRI,'(1X,A1,A50,A7,F12.8,A2)') '|', 205 * ' Epol = - 1/2*Sum_a[ MYind_a*(E^s_a + E^ext)] ', 206 * '| ',EMMPOL,' |' 207 WRITE(LUPRI,'(1X,72A)') ('-',J=1,72) 208 WRITE(LUPRI,'(1X,A1,A50,A7,F12.8,A2)') '|', 209 * ' E_field = - Mu_perm* E^ext ', 210 * '| ',PEDIP1,' |' 211 ELSE 212 WRITE(LUPRI,'(1X,72A)') ('-',J=1,72) 213 WRITE(LUPRI,'(1X,A1,A50,A7,F12.8,A2)') '|', 214 * ' Epol = - 1/2*Sum_a[ MYind_a*E^site_a ] ', 215 * '| ',EMMPOL,' |' 216 END IF 217 WRITE(LUPRI,'(1X,72A)') ('-',J=1,72) 218 WRITE(LUPRI,'(1X,A1,A50,A7,F12.8,A2)') '|', 219 * ' Evdw = Sum_a[ A_ma/|R_ma|^12 - B_ma/|R_ma|^6 ] ', 220 * '| ',EMMVDW,' |' 221 WRITE(LUPRI,'(1X,72A)') ('-',J=1,72) 222 WRITE(LUPRI,'(1X,A1,50A,A1,A20)') 223 * '|',('*',J=1,50),'|', 224 * '*******************|' 225 WRITE(LUPRI,'(1X,72A)') ('-',J=1,72) 226 IF (.NOT.(RELMOM)) THEN 227 WRITE(LUPRI,'(1X,A1,A50,A7,F12.8,A2)') '|', 228 * ' E(MM/MM) = Eelec + Epol + Evdw ', 229 * '| ',EMM_MM,' |' 230 WRITE(LUPRI,'(1X,72A,//)') ('-',J=1,72) 231 ELSE 232 WRITE(LUPRI,'(1X,A1,A50,A7,F12.8,A2)') '|', 233 * ' E(MM/MM) = Eelec + Epol + E_field + Evdw ', 234 * '| ',EMM_MM,' |' 235 WRITE(LUPRI,'(1X,72A,//)') ('-',J=1,72) 236 END IF 237 END IF 238C 239C-------------------------------------------- 240C Second the CC/MM interaction energy terms: 241C-------------------------------------------- 242C 243 WRITE(LUPRI,'(1X,A1,70A,A1)') '+',('=',J=1,70),'+' 244C 245 IF (MODEL .EQ. 'CC2')THEN 246 WRITE(LUPRI,'(1X,A1,A15,A18,A3,A19,A15,A1)') 247 * '|','--------------- ', 248 * 'Final output from ',MODEL,'/MM energy program ', 249 * '---------------','|' 250 ELSE IF (MODEL .EQ. 'CCSD') THEN 251 WRITE(LUPRI,'(1X,A1,A15,A18,A4,A19,A14,A1)') 252 * '|','-------------- ', 253 * 'Final output from ',MODEL,'/MM energy program ', 254 * '--------------','|' 255 END IF 256C 257 WRITE(LUPRI,'(1X,A1,70A,A1)') '+',('=',J=1,70),'+' 258C 259 WRITE(LUPRI,'(1X,A2,A15,A3,A15,A3,A15,A3,A14,A2)') 260 * '| ',' Eelec ',' | ',' Epol ', 261 * ' | ',' Evdw ',' | ',' E(QM/MM) ', 262 * ' |' 263 WRITE(LUPRI,'(1X,A1,70A,A1)') '+',('-',J=1,70),'+' 264 WRITE(LUPRI,'(1X,A2,F15.10,A3,F15.10,A3,F15.10, 265 * A3,F14.10,A2)') 266 * '| ',EEL,' | ',EPOL,' | ',ECLVDW,' | ',ECMCON, 267 * ' |' 268 WRITE(LUPRI,'(1X,A1,70A,A1)') '+',('-',J=1,70),'+' 269C 270 WRITE(LUPRI,'(1X,A1,70A,A1)') '+',('=',J=1,70),'+' 271C 272 WRITE(LUPRI,'(1X,A1,70A,A1)') '+',('-',J=1,70),'+' 273 WRITE(LUPRI,'(1X,A2,A15,A3,A15,A3,A15,A3,A14,A2)') 274 * '| ',' <L|H(vac)|CC> ',' | ','<H(qm)+H(qmmm)>', 275 * ' | ','Delta E(mm/mm) ',' | ','E_rep ', 276 * ' |' 277 WRITE(LUPRI,'(1X,A1,70A,A1)') '+',('-',J=1,70),'+' 278 WRITE(LUPRI,'(1X,A2,F15.10,A3,F15.10,A3,F15.10, 279 * A3,F14.10,A2)') 280 * '| ',ECCGRS1,' | ',ECMCU,' | ',DEMMMM,' | ',EREP, 281 * ' |' 282 WRITE(LUPRI,'(1X,A1,70A,A1//)') '+',('=',J=1,70),'+' 283C 284 IF (RELMOM) THEN 285 WRITE(LUPRI,*)' ' 286 WRITE(LUPRI,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' 287 WRITE(LUPRI,*)' ' 288 WRITE(LUPRI,*)'Pind_a is determined without external field' 289 WRITE(LUPRI,*)'Delta E(mm/mm) is WRONG !!! ' 290 WRITE(LUPRI,*)' ' 291 WRITE(LUPRI,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' 292 WRITE(LUPRI,*)' ' 293 END IF 294C------------------------------------------------------ 295C Third information concerning the convergence is 296C written: 297C------------------------------------------------------ 298C 299 IF (LOITER) THEN 300 WRITE(LUPRI,'(1X,A,I3,A)') 301 * 'Maximum inner iterations for t set to ', 302 * MXTINIT,' in each outer iteration' 303 WRITE(LUPRI,'(1X,A,I3,A/)') 304 * 'Maximum inner iterations for t-bar set to ', 305 * MXLINIT,' in each outer iteration' 306 END IF 307C 308 WRITE(LUPRI,'(/,1X,A,I3,A)') 309 * ' CCMM equations are converged in ',ICCSLIT, 310 * ' outer iterations' 311 WRITE(LUPRI,'(1X,A,I3,A)') 312 * ' CCMM equations are converged in ',NSLVINIT, 313 * ' inner iterations' 314C 315 WRITE(LURES,'(12X,A4,A,F20.10)') 316 * MODELPRI,' Total energy: ',ECMCU 317C 318 WRITE(LURES,'(12X,A4,A,F20.10,/)') 319 * MODELPRI,' E(QM/MM) : ',ECMCON 320C 321 ECCGRS1 = ECMCU 322 LABEL = MODELPRI//'/MM ' 323 CALL WRIPRO(ECCGRS1,MODEL,0, 324 * LABEL,LABEL,LABEL,LABEL, 325 * ZERO,ZERO,ZERO,1,0,0,0) 326C 327 LABEL = 'E_QM/MM ' 328 CALL WRIPRO(ECMCON,MODEL,0, 329 * LABEL,LABEL,LABEL,LABEL, 330 * ZERO,ZERO,ZERO,1,0,0,0) 331 332C------------------------------------------------------ 333C Print information concerning the reaction field 334C------------------------------------------------------ 335 336 CALL RFQMMM() 337 338C------------------------------------------------------ 339 340 ELSE 341C 342 ICCSLIT = ICCSLIT + 1 343C 344 IF (ICCSLIT.GT.MXCCSLIT) THEN 345 WRITE(LUPRI,*) 'Maximum number of CCMM iterations', 346 * MXCCSLIT,' is reached' 347 CALL QUIT( 'Maximum number of CCMM iterations reached') 348 ENDIF 349C 350 ENDIF 351C 352C-------------------------------------- 353C Ending header for each iteration: 354C-------------------------------------- 355C 356 WRITE (LUPRI,'(1X,A)') 357 *'*****************************************************'// 358 *'************' 359 WRITE (LUPRI,'(1X,A)') 360 *'******* End of Coupled Cluster/Molecular Mechanics program'// 361 *' ******' 362 WRITE (LUPRI,'(1X,A)') 363 *'*****************************************************'// 364 *'************' 365C 366 END 367C 368C************************************************************** 369C /* Deck cc_qm3e */ 370 SUBROUTINE CC_QM3E(ECCMM,EPOL,EEL,EELEL,EREP,AODEN,WORK,LWORK) 371C************************************************************** 372C 373C QM3: JK+AO+OC+KMI,01 374C Purpose: Only calculation of E(QM/MM) 375C QMMM: sneskov, 09 376C Purpose: Updates the induced dipoles besides calculating E(QM/MM) 377C 378C In new QMMM implementation (NYQMMM) the strategy is as follows 379C 1) Calculate F^sta=F^el+F^mul+F^nuc 380C 2) Calculate induced dipole by transforming with relay 381C matrix(MMMAT) or by iterative means (MMITER) 382C 3) Put induced moments to file 383C 4) Calculate Epol=-1/2sum_a\mu^ind_aF^sta_a 384C 5) Calculate Eelesta=-sum_s<N_s> 385C 386C Also the LGFILE keyword is reset such that a new G matrix is 387C written to file when TGT is called next time. 388C 389C************************************************************** 390C 391#include "implicit.h" 392#include "priunit.h" 393#include "dummy.h" 394#include "mxcent.h" 395#include "qmmm.h" 396#include "qm3.h" 397#include "iratdef.h" 398#include "maxorb.h" 399#include "orgcom.h" 400#include "nuclei.h" 401#include "ccinftap.h" 402#include "ccorb.h" 403#include "ccsdsym.h" 404#include "ccsdio.h" 405#include "ccsdinp.h" 406#include "ccfield.h" 407#include "exeinf.h" 408#include "ccfdgeo.h" 409#include "thrzer.h" 410 411C 412 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, IZERO = 0 , TWO = 2.0D0) 413 PARAMETER (HALF = 1.0D0/2.0D0) 414 PARAMETER (D3I = 1.0D0/3.0D0, D6I = 1.0D0/6.0D0) 415C 416 DIMENSION WORK(LWORK),AODEN(*) 417 CHARACTER LABEL*8 418 419 DIMENSION CCNS(MXQM3),ALPHAMM(MXQM3),FFJ(3) 420 INTEGER KS 421C 422 CALL QENTER('CC_QM3E') 423C------------------------------------------------------------------ 424C Init parameters 425C------------------------------------------------------------------ 426 427 ISYMOP = 1 428C 429C------------------------------------------------------------------ 430C Dynamical allocation for CCMM 431C------------------------------------------------------------------ 432C 433 434C------------------------------------------------------------------ 435C Dynamical allocation. If all models are SPC (LOSPC = .TRUE.) 436C the electric field from the electrons are not needed. Then 437C only space for Ns in AO basis is allocated. 438C------------------------------------------------------------------ 439C 440 441 IF ( .NOT. (LOSPC) ) THEN 442C 443 KRAAOx = 1 444 KRAAOy = KRAAOx + N2BST(ISYMOP) 445 KRAAOz = KRAAOy + N2BST(ISYMOP) 446 KOMMSx = KRAAOz + N2BST(ISYMOP) 447 KOMMSy = KOMMSx + NCOMS 448 KOMMSz = KOMMSy + NCOMS 449 KNSAO = KOMMSz + NCOMS 450 KRAx = KNSAO + N2BST(ISYMOP) 451 KRAy = KRAx + NCOMS 452 KRAz = KRAy + NCOMS 453 KWRK1 = KRAz + NCOMS 454 LWRK1 = LWORK - KWRK1 455C 456 CALL DZERO(WORK(KRAAOx),3*N2BST(ISYMOP)) 457 CALL DZERO(WORK(KOMMSx),3*NCOMS) 458 CALL DZERO(WORK(KNSAO),N2BST(ISYMOP)) 459 CALL DZERO(WORK(KRAx),3*NCOMS) 460C 461 ELSE 462C 463 KNSAO = 1 464 KWRK1 = KNSAO + N2BST(ISYMOP) 465 LWRK1 = LWORK - KWRK1 466C 467 CALL DZERO(WORK(KNSAO),N2BST(ISYMOP)) 468C 469 END IF 470C 471 IF (LWRK1.LT.0) CALL QUIT( 'Too little work in CC_QM3RAINT') 472C 473C---------------------------- 474C Read Rr(a) in AO basis: 475C---------------------------- 476C 477C First we see if any static fields are to be added 478C 479 FFJ(1) = 0.0D0 480 FFJ(2) = 0.0D0 481 FFJ(3) = 0.0D0 482C 483 IF (RELMOM) THEN 484 DO 330 IF =1, NFIELD 485 IF (LFIELD(IF).EQ. 'XDIPLEN') FFJ(1)=FFJ(1)+EFIELD(IF) 486 IF (LFIELD(IF).EQ. 'YDIPLEN') FFJ(2)=FFJ(2)+EFIELD(IF) 487 IF (LFIELD(IF).EQ. 'ZDIPLEN') FFJ(3)=FFJ(3)+EFIELD(IF) 488 330 CONTINUE 489 END IF 490C 491 IF (.NOT. LOSPC) THEN 492 CALL CC_QM3RAINT(WORK(KRAAOx),WORK(KRAAOy),WORK(KRAAOz), 493 & AODEN,WORK(KRAx),WORK(KRAy),WORK(KRAz), 494 & WORK(KWRK1),LWRK1,ISYMOP) 495 END IF 496C 497C------------------------------ 498C Read N(s) in AO basis: 499C------------------------------ 500C 501 CALL CC_QM3NSINT(WORK(KNSAO),AODEN,CCNS, 502 & WORK(KWRK1),LWRK1,ISYMOP) 503C 504C------------------------------------- 505C Add up the energy contributions: 506C 1) Epol: 507C------------------------------------- 508C 509 ECCMM = ZERO 510C 511 IF (.NOT. LOSPC) THEN 512 CALL CC_GET31('OBARFILE',NCOMS, 513 * WORK(KOMMSx),WORK(KOMMSy),WORK(KOMMSz)) 514C 515 KS = 0 516 517 DO 110 I = 1, ISYTP 518 IF (MDLWRD(I)(1:5) .EQ. 'SPC_E') THEN 519 DO 120 J = NSYSBG(I), NSYSED(I) 520 DO 121 K = 1, NUALIS(I) 521 KS = KS+1 522 RAT = ZERO 523 RAT = 0.5D0 * ALPIMM(I,K)*WORK(KRAx+KS-1) * 524 & (WORK(KRAx + KS -1) + WORK(KOMMSx + KS - 1)) 525 & + 0.5D0 * ALPIMM(I,K) * WORK(KRAy + KS -1) * 526 & (WORK(KRAy + KS -1) + WORK(KOMMSy + KS - 1)) 527 & + 0.5D0 * ALPIMM(I,K) * WORK(KRAz + KS -1) * 528 & (WORK(KRAz + KS -1) + WORK(KOMMSz + KS - 1)) 529 ECCMM = ECCMM - RAT 530 121 CONTINUE 531 120 CONTINUE 532 END IF 533 110 CONTINUE 534C 535 TEMP1 = ECCMM 536 CALL QM3_OTILDE(OTILDE,FFJ) 537 ELSE 538 OTILDE = 0.0D0 539 TEMP1 = 0.0D0 540 END IF 541C 542 ECCMM = ECCMM + OTILDE 543 EPOL = ECCMM 544C 545C-------- 546C 2) Eel: 547C-------- 548C 549 ETEMP = 0.0D0 550 L = 0 551C 552 DO 130 I = 1, ISYTP 553 IF (MDLWRD(I)(1:3) .EQ. 'SPC') THEN 554 DO 140 J = NSYSBG(I), NSYSED(I) 555 DO 150 K = 1, NSISY(I) 556 L = L + 1 557 ETEMP = ETEMP - CCNS(L) 558 150 CONTINUE 559 140 CONTINUE 560 END IF 561 130 CONTINUE 562C 563 EELEL = ETEMP 564 565 ECCMM = ECCMM + ENUQM3 + EELEL 566C 567 EEL = ENUQM3 + EELEL 568C 569C------------- 570C 3) E(VDW): 571C------------- 572C 573 ECCMM = ECCMM + ECLVDW 574C 575C----------------------------------------------- 576C 4) E_repulsion 577C 578 CALL CCMM_REP2(EREP,AODEN,WORK(KWRK1),LWRK1) 579 ECCMM = ECCMM + EREP 580C 581C----------------------------------------------- 582C Testing the energy contributions one by one: 583C----------------------------------------------- 584C 585 IF (IQM3PR .GT. 5) THEN 586 WRITE(LUPRI,'(//,A)') 587 WRITE(LUPRI,'(/A)') 588 * ' +=====================================' 589 * //'======================================' 590 * //'==================================+' 591 WRITE(LUPRI,'(1X,A)') 592 * '|--------- ' 593 * //'The different energy contributions outlined' 594 * //' ---------|' 595 WRITE(LUPRI,'(A)') 596 * ' +=====================================' 597 * //'======================================' 598 * //'==================================+' 599 WRITE(LUPRI,'(1X,A)') 600 * '| Evdw | E-nuclear |' 601 * //' -Sum_s <N_s> |' 602 * //'| -alpha*Sum_a[<Rra> | O-tilde |' 603 WRITE(LUPRI,'(1X,A)') 604 * '| contribution | contribution |' 605 * //' |' 606 * //'*{<Rra>+OmmS}] | contribution |' 607 WRITE(LUPRI,'(A)') 608 * ' +-------------------------------------' 609 * //'--------------------------------------' 610 * //'----------------------------------+' 611 WRITE(LUPRI,'(1X,A,F20.16,A,F20.16,A,F20.16,A, 612 & F20.16,A,F20.16,A)') 613 * '|',ECLVDW,' |',ENUQM3,' |',EELEL,' |',TEMP1,' |', 614 * OTILDE,' |' 615 WRITE(LUPRI,'(A)') 616 * ' +=====================================' 617 * //'======================================' 618 * //'==================================+' 619 WRITE(LUPRI,'(//,A)') 620 END IF 621 622 CALL QEXIT('CC_QM3E') 623C 624 END 625C 626C************************************************************** 627C /* Deck cc_qm3raint */ 628 SUBROUTINE CC_QM3RAINT(RAAOx,RAAOy,RAAOz,AODEN, 629 & RAx,RAy,RAz,WORK,LWORK, 630 & ISYMT) 631C************************************************************** 632C 633C Readin mm-gradient integrals in coupled cluster format 634C and evaluate the expectation value of Rra. Calculates 635C the induced dipole moment at the center of mass a and 636C evaluates the O(mmss) factor at the center of mass a 637C 638C************************************************************** 639C 640#include "implicit.h" 641#include "maxorb.h" 642#include "mxcent.h" 643#include "nuclei.h" 644#include "dummy.h" 645#include "iratdef.h" 646#include "priunit.h" 647#include "ccorb.h" 648#include "ccsdsym.h" 649#include "ccsdio.h" 650#include "ccsdinp.h" 651#include "ccfield.h" 652#include "exeinf.h" 653#include "ccfdgeo.h" 654#include "qm3.h" 655#include "ccslvinf.h" 656#include "ccinftap.h" 657C 658 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, IZERO = 0 , TWO = 2.0D0) 659C 660 DIMENSION WORK(LWORK) 661 DIMENSION RAAOx(*),RAAOy(*),RAAOz(*) 662 DIMENSION AODEN(*) 663 DIMENSION RAx(*) 664 DIMENSION RAy(*) 665 DIMENSION RAz(*) 666 DIMENSION EELEC(3,MXQM3) 667 DIMENSION FFJ(3) 668 CHARACTER*8 LABEL 669 INTEGER KL 670 LOGICAL LOINDM 671C 672C-------------------------------------------------------- 673C Initializing the electronic electric field to zero: 674C-------------------------------------------------------- 675C 676 DO 879 I = 1, MXQM3 677 DO 880 J = 1, 3 678 EELEC(J,I) = 0.0D0 679 880 CONTINUE 680 879 CONTINUE 681C 682 IF (IQM3PR .GT. 10) THEN 683 WRITE(LUPRI,*) 'CC_QM3RAINT: Read in integrals' 684 WRITE(LUPRI,*) 'Input symmetry claimed', isymt 685 ENDIF 686C 687 LUCCEF = -1 688 CALL GPOPEN(LUCCEF,'ELFDMM','OLD',' ', 689 & 'UNFORMATTED',IDUMMY,.FALSE.) 690 REWIND(LUCCEF) 691C 692 KL = 0 693C 694 DO 950 I = 1, ISYTP 695 IF (MDLWRD(I)(1:5) .EQ. 'SPC_E') THEN 696 DO 951 J = NSYSBG(I), NSYSED(I) 697 DO 952 K = 1, NUALIS(I) 698C 699 KL = KL +1 700C 701 LABEL = 'READ INT' 702 CALL CCMMINT(LABEL,LUCCEF,RAAOx,DUMMY,IRREP, 703 * ISYM,IERR,WORK,LWORK) 704 RAx(KL) = DDOT(N2BST(ISYMOP),RAAOx,1,AODEN,1) 705C 706 LABEL = 'READ INT' 707 CALL CCMMINT(LABEL,LUCCEF,RAAOy,DUMMY,IRREP, 708 * ISYM,IERR,WORK,LWORK) 709 RAy(KL) = DDOT(N2BST(ISYMOP),RAAOy,1,AODEN,1) 710C 711 LABEL = 'READ INT' 712 CALL CCMMINT(LABEL,LUCCEF,RAAOz,DUMMY,IRREP, 713 * ISYM,IERR,WORK,LWORK) 714 RAz(KL) = DDOT(N2BST(ISYMOP),RAAOz,1,AODEN,1) 715C 716 IF (IQM3PR .GT. 10) THEN 717 WRITE(LUPRI,*)'RAx(KL) =',RAx(KL) 718 WRITE(LUPRI,*)'RAy(KL) =',RAy(KL) 719 WRITE(LUPRI,*)'RAz(KL) =',RAz(KL) 720 END IF 721C 722 952 CONTINUE 723 951 CONTINUE 724 END IF 725 950 CONTINUE 726C 727 IF (KL .NE. NCOMS) THEN 728 CALL QUIT('Error in no. of center of masses in CC_QM3RAINT') 729 END IF 730C 731 DO 534 LM = 1,NCOMS 732 EELEC(1,LM) = RAx(LM) 733 EELEC(2,LM) = RAy(LM) 734 EELEC(3,LM) = RAz(LM) 735 534 CONTINUE 736C 737 CALL GPCLOSE(LUCCEF,'KEEP') 738C 739C If RELMOM is true we want to include the external field(s) in 740C the determination of the induced dipole moments 741C 742 FFJ(1) = 0.0D0 743 FFJ(2) = 0.0D0 744 FFJ(3) = 0.0D0 745C 746 IF (RELMOM) THEN 747 DO 330 IF =1, NFIELD 748 IF (LFIELD(IF) .EQ. 'XDIPLEN') FFJ(1) = FFJ(1) + EFIELD(IF) 749 IF (LFIELD(IF) .EQ. 'YDIPLEN') FFJ(2) = FFJ(2) + EFIELD(IF) 750 IF (LFIELD(IF) .EQ. 'ZDIPLEN') FFJ(3) = FFJ(3) + EFIELD(IF) 751 330 CONTINUE 752 END IF 753C 754 IF (FIXMOM) THEN 755 WRITE(LUPRI,'(/A)')'FIXMOM: NO ITER. DETERM. OF MYIND' 756 ELSE IF (LOSPC) THEN 757 WRITE(LUPRI,'(/A)')'All MM models are SPC, INDMOM not called' 758 ELSE 759 LOINDM = .FALSE. 760 CALL INDMOM(EELEC,LOINDM,FFJ) 761 END IF 762C 763 CALL QM3_OBAR(FFJ) 764 CALL CC_PUT31('CC_RA',NCOMS,RAx,RAy,RAz) 765C 766 END 767C 768C************************************************************** 769C /* Deck cc_qm3nsint */ 770 SUBROUTINE CC_QM3NSINT(NSAO,AODEN,CCNS,WORK,LWORK,ISYMT) 771C************************************************************** 772C 773C Readin mm-potentiale energy integrals in coupled cluster 774C format and calculate the expectation value of Ns. 775C 776C************************************************************** 777C 778#include "implicit.h" 779#include "maxorb.h" 780#include "mxcent.h" 781#include "nuclei.h" 782#include "dummy.h" 783#include "iratdef.h" 784#include "priunit.h" 785#include "ccorb.h" 786#include "ccsdsym.h" 787#include "ccsdio.h" 788#include "ccsdinp.h" 789#include "ccfield.h" 790#include "exeinf.h" 791#include "ccfdgeo.h" 792#include "ccinftap.h" 793#include "qm3.h" 794C 795 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, IZERO = 0 , TWO = 2.0D0) 796C 797 DIMENSION WORK(LWORK),NSAO(*) 798 DIMENSION AODEN(*) 799 DIMENSION CCNS(*) 800C 801 CHARACTER*8 LABEL 802C 803 LUCCPO = -1 804 CALL GPOPEN(LUCCPO,'POTMM','UNKNOWN',' ', 805 & 'UNFORMATTED',IDUMMY,.FALSE.) 806 REWIND(LUCCPO) 807C 808 LM = 0 809C 810 DO 940 I= 1, ISYTP 811 IF (MDLWRD(I)(1:3) .EQ. 'SPC') THEN 812 DO 941 J = NSYSBG(I), NSYSED(I) 813 DO 942 K = 1, NSISY(I) 814 LM = LM + 1 815 LABEL = 'READ INT' 816 CALL CCMMINT(LABEL,LUCCPO,NSAO,DUMMY,IRREP, 817 & ISYM,IERR,WORK,LWORK) 818 CCNS(LM) = DDOT(N2BST(ISYMOP),NSAO,1,AODEN,1) 819 942 CONTINUE 820 941 CONTINUE 821 END IF 822 940 CONTINUE 823C 824 CALL GPCLOSE(LUCCPO,'KEEP') 825C 826C------------------- 827C Print section: 828C------------------- 829C 830 IF (IQM3PR .GT. 5) THEN 831 WRITE(LUPRI,'(//,A)') 832 * ' +======================+' 833 WRITE(LUPRI,'(A)') 834 * ' |Site| <N_s> |' 835 WRITE(LUPRI,'(1X,A)') 836 * '+======================+' 837C 838 LS = 0 839C 840 DO 215 I = 1, ISYTP 841 IF (MDLWRD(I)(1:3) .EQ. 'SPC') THEN 842 DO 216 J = NSYSBG(I), NSYSED(I) 843 DO 217 K = 1, NSISY(I) 844C 845 LS = LS + 1 846C 847 WRITE(LUPRI,'(1X,A,I3,A,F16.10,A)') 848 * '| ', LS,'|', CCNS(LS),' |' 849 WRITE(LUPRI,'(1X,A)') 850 * '+----------------------+' 851 217 CONTINUE 852 216 CONTINUE 853 END IF 854 215 CONTINUE 855 WRITE(LUPRI,'(//,A)') 856 END IF 857C 858 END 859C 860C************************************************************** 861C /* Deck ccmm_rhstg */ 862 SUBROUTINE CCMM_RHSTG(FOCK,WORK,LWORK) 863C************************************************************** 864C 865C Direct calculation of CC/MM solvent effects. If 866C OLDTG = .TRUE. the Ns part is included in the T^g 867C operator. However, if OLDTG = .FALSE. the contrinution is 868C allready in the h1 "vacuum" operator. 869C 870C************************************************************** 871C 872#include "implicit.h" 873#include "maxorb.h" 874#include "mxcent.h" 875#include "dummy.h" 876#include "iratdef.h" 877#include "priunit.h" 878#include "ccorb.h" 879#include "ccsdsym.h" 880#include "ccsdio.h" 881#include "ccsdinp.h" 882#include "ccfield.h" 883#include "exeinf.h" 884#include "ccfdgeo.h" 885#include "qm3.h" 886#include "ccinftap.h" 887C 888 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, IZERO = 0 , TWO = 2.0D0) 889C 890 DIMENSION WORK(LWORK),FOCK(*) 891 DIMENSION ALPHAMM(MXQM3) 892C 893 CHARACTER*8 LABEL 894 CHARACTER*9 FILMMR 895C 896 LOGICAL FIRST 897 SAVE FIRST 898C 899C------------------------------------------------------------------ 900C If a system is model SPC_ECX (X=1,2,3,4) the polarization of 901C the system does not contribute to the optimization of the 902C wave function. Also, if OLDTG = .FALSE. the partial charges 903C are included in h1. So if all models are SPC -> RETURN from 904C CCMM_RHSTG. 905C------------------------------------------------------------------ 906C 907 IF ( (.NOT. (OLDTG)) .AND. (LOSPC) ) RETURN 908C 909C---------------------------- 910C Dynamical allocation: 911C---------------------------- 912C 913 IF ( (.NOT. (OLDTG)) .AND. (.NOT. (LOSPC)) ) THEN 914C 915C-------------------------- 916C No space needed for 917C Ns in AO basis: 918C-------------------------- 919C 920 KTAO = 1 921 KRAAOx = KTAO + NNBASX 922 KRAAOy = KRAAOx + NNBASX 923 KRAAOz = KRAAOy + NNBASX 924 KENSAx = KRAAOz + NNBASX 925 KENSAy = KENSAx + NCOMS 926 KENSAz = KENSAy + NCOMS 927 KRAx = KENSAz + NCOMS 928 KRAy = KRAx + NCOMS 929 KRAz = KRAy + NCOMS 930 KWRK1 = KRAz + NCOMS 931 LWRK1 = LWORK - KWRK1 932C 933 CALL DZERO(WORK(KTAO),4*NNBASX) 934 CALL DZERO(WORK(KENSAx),3*NCOMS) 935 CALL DZERO(WORK(KRAx),3*NCOMS) 936C 937 ELSE IF ( (OLDTG) .AND. (.NOT. (LOSPC)) ) THEN 938C 939C-------------------------- 940C Space needed for 941C Ns in AO basis: 942C-------------------------- 943C 944 KTAO = 1 945 KRAAOx = KTAO + NNBASX 946 KRAAOy = KRAAOx + NNBASX 947 KRAAOz = KRAAOy + NNBASX 948 KENSAx = KRAAOz + NNBASX 949 KENSAy = KENSAx + NCOMS 950 KENSAz = KENSAy + NCOMS 951 KNSAO = KENSAz + NCOMS 952 KRAx = KNSAO + NNBASX 953 KRAy = KRAx + NCOMS 954 KRAz = KRAy + NCOMS 955 KWRK1 = KRAz + NCOMS 956 LWRK1 = LWORK - KWRK1 957C 958 CALL DZERO(WORK(KTAO),4*NNBASX) 959 CALL DZERO(WORK(KENSAx),3*NCOMS) 960 CALL DZERO(WORK(KNSAO),NNBASX) 961 CALL DZERO(WORK(KRAx),3*NCOMS) 962C 963 ELSE IF ( (OLDTG) .AND. (LOSPC) ) THEN 964C 965C-------------------------- 966C No space needed for 967C alpha contributers: 968C-------------------------- 969C 970 KTAO = 1 971 KNSAO = KTAO + NNBASX 972 KWRK1 = KNSAO + NNBASX 973 LWRK1 = LWORK - KWRK1 974C 975 CALL DZERO(WORK(KTAO),2*NNBASX) 976C 977 END IF 978C 979 IF (LWRK1.LT.0) CALL QUIT( 'Too little work in CCMM_RHSTG' ) 980C 981C------------------------------------------------- 982C If OLDTG and LOSPC are .TRUE. -> GOTO 888 983C which is after the alpha part in CCMM_RHSTG. 984C------------------------------------------------- 985C 986 IF ( (OLDTG) .AND. (LOSPC) ) GOTO 888 987C 988 CALL CC_GET31('CC_RA',NCOMS, 989 * WORK(KRAx),WORK(KRAy),WORK(KRAz)) 990C 991C------------------- 992C Print section: 993C------------------- 994C 995 IF (IQM3PR .GT. 10) THEN 996 WRITE(LUPRI,'(//,A)') 997 * ' +============================================' 998 * //'==============+' 999 WRITE(LUPRI,'(A)') 1000 * ' | COM| <Rra>_x | <Rra>_y |' 1001 * //' <Rra>_z |' 1002 WRITE(LUPRI,'(1X,A)') 1003 * '+============================================' 1004 * //'==============+' 1005C 1006 NUM = 0 1007C 1008 DO 215 I = 1, ISYTP 1009 IF (MDLWRD(I)(1:5) .EQ. 'SPC_E') THEN 1010 DO 216 J = NSYSBG(I), NSYSED(I) 1011 DO 217 K=1,NUALIS(I) 1012C 1013 NUM = NUM + 1 1014C 1015 WRITE(LUPRI,'(1X,A,I3,A,F16.10,A, 1016 & F16.10,A,F16.10,A)') 1017 * '| ', J,'|', WORK(KRAx + NUM - 1),' |', 1018 * WORK(KRAy + NUM - 1),' |', 1019 * WORK(KRAz + NUM - 1),' |' 1020 WRITE(LUPRI,'(1X,A)') 1021 * '+---------------------------------------------' 1022 * //'-------------+' 1023 217 CONTINUE 1024 216 CONTINUE 1025 END IF 1026 215 CONTINUE 1027 WRITE(LUPRI,'(//,A)') 1028 END IF 1029C 1030 CALL CC_GET31('ENSAFILE',NCOMS, 1031 * WORK(KENSAx),WORK(KENSAy),WORK(KENSAz)) 1032 1033C 1034 IF (IQM3PR .GT. 10) THEN 1035 WRITE(LUPRI,'(//,A)') 1036 * ' +============================================' 1037 * //'==============+' 1038 WRITE(LUPRI,'(A)') 1039 * ' | COM| ENSA_x | ENSA_y |' 1040 * //' ENSA_z |' 1041 WRITE(LUPRI,'(1X,A)') 1042 * '+============================================' 1043 * //'==============+' 1044C 1045 NUM = 0 1046C 1047 DO 415 I = 1, ISYTP 1048 IF (MDLWRD(I)(1:5) .EQ. 'SPC_E') THEN 1049C 1050 NUM = NUM + 1 1051C 1052 DO 416 J = NSYSBG(I), NSYSED(I) 1053 DO 417 K=1,NUALIS(I) 1054 WRITE(LUPRI,'(1X,A,I3,A,F16.10,A, 1055 & F16.10,A,F16.10,A)') 1056 * '| ', J,'|', WORK(KENSAx + NUM - 1),' |', 1057 * WORK(KENSAy + NUM - 1),' |', 1058 * WORK(KENSAz + NUM - 1),' |' 1059 WRITE(LUPRI,'(1X,A)') 1060 * '+---------------------------------------------' 1061 * //'-------------+' 1062 417 CONTINUE 1063 416 CONTINUE 1064 END IF 1065 415 CONTINUE 1066 WRITE(LUPRI,'(//,A)') 1067 END IF 1068C 1069C--------------------------------------- 1070C Read Rr(s) in ao basis and add to tao: 1071C--------------------------------------- 1072C 1073 LUCCEF = -1 1074 CALL GPOPEN(LUCCEF,'ELFDMM','OLD',' ', 1075 & 'UNFORMATTED',IDUMMY,.FALSE.) 1076 REWIND(LUCCEF) 1077C 1078 LM = 0 1079C 1080 DO 520 I = 1, ISYTP 1081 IF (MDLWRD(I)(1:5) .EQ. 'SPC_E') THEN 1082 DO 521 J = NSYSBG(I), NSYSED(I) 1083 DO 522 K = 1, NUALIS(I) 1084 LM = LM + 1 1085C 1086 CALL READT(LUCCEF,NNBASX,WORK(KRAAOx)) 1087C 1088 IF (IQM3PR .GE. 15) THEN 1089 WRITE (LUPRI,'(/A)') ' Rra_ao x matrix:' 1090 CALL OUTPAK(WORK(KRAAOx),NBAST,1,LUPRI) 1091 END IF 1092C 1093 CALL READT(LUCCEF,NNBASX,WORK(KRAAOy)) 1094C 1095 IF (IQM3PR .GE. 15) THEN 1096 WRITE (LUPRI,'(/A)') ' Rra_ao y matrix:' 1097 CALL OUTPAK(WORK(KRAAOy),NBAST,1,LUPRI) 1098 END IF 1099C 1100 CALL READT(LUCCEF,NNBASX,WORK(KRAAOz)) 1101C 1102 IF (IQM3PR .GE. 15) THEN 1103 WRITE (LUPRI,'(/A)') ' Rra_ao z matrix:' 1104 CALL OUTPAK(WORK(KRAAOz),NBAST,1,LUPRI) 1105 END IF 1106C 1107 IF (MDLWRD(I) .EQ. 'SPC_E01') THEN 1108 FACx = - ALPIMM(I,K) * 1109 & (WORK(KRAx + LM - 1) 1110 & + 0.5D0 * WORK(KENSAx + LM - 1)) 1111 FACy = - ALPIMM(I,K) * 1112 & (WORK(KRAy + LM - 1) 1113 & + 0.5D0 * WORK(KENSAy + LM - 1)) 1114 FACz = - ALPIMM(I,K) * 1115 & (WORK(KRAz + LM - 1) 1116 & + 0.5D0 * WORK(KENSAz + LM - 1)) 1117C 1118 CALL DAXPY(NNBASX,FACx,WORK(KRAAOx), 1119 * 1,WORK(KTAO),1) 1120C 1121 CALL DAXPY(NNBASX,FACy,WORK(KRAAOy), 1122 * 1,WORK(KTAO),1) 1123C 1124 CALL DAXPY(NNBASX,FACz,WORK(KRAAOz), 1125 * 1,WORK(KTAO),1) 1126 END IF 1127 522 CONTINUE 1128 521 CONTINUE 1129 END IF 1130 520 CONTINUE 1131C 1132 IF (IQM3PR.GT.14) THEN 1133 WRITE (LUPRI,*) ' NORM of TAO matrix /alpha contr.: ', 1134 * DDOT(NNBASX,WORK(KTAO),1,WORK(KTAO),1) 1135 WRITE (LUPRI,'(/A)') ' TAO matrix: ' 1136 CALL OUTPAK(WORK(KTAO),NBAST,1,LUPRI) 1137 ENDIF 1138C 1139 CALL GPCLOSE(LUCCEF,'KEEP') 1140C 1141C-------------------------------------- 1142C End of alpha part in CCMM_RHSTG: 1143C-------------------------------------- 1144C 1145 888 CONTINUE 1146C 1147C----------------------------------------------------------- 1148C If OLDTG we want to include the partial MM charges in 1149C the T^g operator in the opt. of the wave function. 1150C----------------------------------------------------------- 1151C 1152 IF (OLDTG) THEN 1153C 1154 LUCCPO = -1 1155 CALL GPOPEN(LUCCPO,'POTMM','UNKNOWN',' ', 1156 & 'UNFORMATTED',IDUMMY,.FALSE.) 1157 REWIND (LUCCPO) 1158C 1159 FAC1 = -1.0D0 1160C 1161 L = 0 1162C 1163 DO 525 I = 1, ISYTP 1164 IF (MDLWRD(I)(1:3) .EQ. 'SPC') THEN 1165 DO 526 J = NSYSBG(I), NSYSED(I) 1166 DO 527 K = 1,NSISY(I) 1167C 1168 L = L +1 1169C 1170 CALL READT(LUCCPO,NNBASX,WORK(KNSAO)) 1171C 1172 IF (IQM3PR .GE. 4) THEN 1173 WRITE (LUPRI,'(/A,I3,A)') 1174 * ' N(',L,')_ao matrix: ' 1175 CALL OUTPAK(WORK(KNSAO),NBAST,1,LUPRI) 1176 END IF 1177C 1178 CALL DAXPY(NNBASX,FAC1,WORK(KNSAO), 1179 * 1,WORK(KTAO),1) 1180 527 CONTINUE 1181 526 CONTINUE 1182 END IF 1183 525 CONTINUE 1184C 1185 IF (IQM3PR.GT.14) THEN 1186 WRITE (LUPRI,*) ' NORM of TAO matrix /nsao: ', 1187 * DDOT(NNBASX,WORK(KTAO),1,WORK(KTAO),1) 1188 WRITE (LUPRI,'(/A)') ' TAO matrix: ' 1189 CALL OUTPAK(WORK(KTAO),NBAST,1,LUPRI) 1190 ENDIF 1191C 1192 CALL GPCLOSE(LUCCPO,'KEEP') 1193C 1194 END IF 1195C 1196C-------------------------------------------------- 1197C Add contribution to effective AO fock matrix: 1198C-------------------------------------------------- 1199C 1200 FACT = 1.0D0 1201C 1202C-------------------------- 1203C Dynamical allocation: 1204C-------------------------- 1205C 1206 KINT = KWRK1 1207 KWRK2 = KINT + N2BST(ISYMOP) 1208 LWRK2 = LWORK - KWRK2 1209C 1210 IF (LWRK2.LT.0) CALL QUIT( 'Too little work in CCMM_RHSTG' ) 1211C 1212 CALL DZERO(WORK(KINT),N2BST(ISYMOP)) 1213C 1214 IF (IQM3PR .GT. 14) THEN 1215 WRITE (LUPRI,*) 'NORM of TAO matrix. CCMMINT: ', 1216 * DDOT(NNBASX,WORK(KTAO),1,WORK(KTAO),1) 1217 WRITE (LUPRI,'(/A)') ' TAO matrix: ' 1218 CALL OUTPAK(WORK(KTAO),NBAST,1,LUPRI) 1219 ENDIF 1220C 1221 LABEL= 'GIVE INT' 1222C 1223 CALL CCMMINT(LABEL,LUCCEF,WORK(KINT),WORK(KTAO), 1224 & IRREP,ISYM,IERR,WORK(KWRK2),LWRK2) 1225C 1226 IF (IQM3PR .GT.14) THEN 1227 CALL AROUND( 'CCMM cont. to AO matrix' ) 1228 CALL CC_PRFCKAO(WORK(KINT),ISYMOP) 1229 ENDIF 1230C 1231 IF (IQM3PR .GT.14) THEN 1232 CALL AROUND( 'Usual Fock AO matrix' ) 1233 CALL CC_PRFCKAO(FOCK,ISYMOP) 1234 ENDIF 1235C 1236 CALL DAXPY(N2BST(ISYMOP),FACT,WORK(KINT),1,FOCK,1) 1237C 1238 IF (IQM3PR .GT.14) THEN 1239 CALL AROUND( 'Modified Fock AO matrix before rep.') 1240 CALL CC_PRFCKAO(FOCK,ISYMOP) 1241 ENDIF 1242C 1243 IF (LOSHAW) THEN 1244 CALL DZERO(WORK(KINT),N2BST(ISYMOP)) 1245C 1246 FILMMR = 'MMREPOP_0' 1247 LUMMRE = -1 1248 IADRFIL = 1 1249 NDATA = N2BST(ISYMOP) 1250C 1251 CALL WOPEN2(LUMMRE,FILMMR,64,0) 1252C 1253 CALL GETWA2(LUMMRE,FILMMR,WORK(KINT),IADRFIL,NDATA) 1254C 1255 IF (REPTST) THEN 1256 NBASLO = MAX(NBAST,1) 1257C 1258 KINT1 = KWRK2 1259 KINT2 = KINT1 + N2BST(ISYMOP) 1260 KWRK3 = KINT2 + N2BST(ISYMOP) 1261 LWRK3 = LWORK - KWRK3 1262C 1263 IF (LWRK3 .LT. 0) THEN 1264 CALL QUIT( 'Too litle work in CCMM_RHSTG Rep. 1') 1265 END IF 1266C 1267 CALL DZERO(WORK(KINT1),2*N2BST(ISYMOP)) 1268 CALL DCOPY(N2BST(ISYMOP),WORK(KINT),1,WORK(KINT1),1) 1269C 1270 DO 668 KR=1,NREPMT 1271 CALL DGEMM('N','N',NBAST,NBAST,NBAST, 1272 * ONE,WORK(KINT),NBASLO, 1273 * WORK(KINT1),NBASLO, 1274 * ZERO,WORK(KINT2),NBASLO) 1275C 1276 CALL DCOPY(N2BST(ISYMOP),WORK(KINT2),1,WORK(KINT),1) 1277 668 CONTINUE 1278 END IF 1279C 1280 TAL1 = DDOT(N2BST(ISYMOP),WORK(KINT),1,WORK(KINT),1) 1281 TAL2 = SQRT(TAL1) 1282C 1283 CALL DAXPY(N2BST(ISYMOP),FACT,WORK(KINT),1,FOCK,1) 1284C 1285 IF (IQM3PR .GT.14) THEN 1286 CALL AROUND( 'Modified Fock AO matrix after rep.') 1287 CALL CC_PRFCKAO(FOCK,ISYMOP) 1288 ENDIF 1289C 1290 CALL WCLOSE2(LUMMRE,FILMMR, 'KEEP') 1291C 1292 END IF !LOSHAW 1293C 1294 END 1295C 1296C************************************************************** 1297C /* Deck ccmm_ltrb */ 1298 SUBROUTINE CCMM_LTRB(RHO1,RHO2,CTR1,CTR2, 1299 & ISYMTR,LR,WORK,LWORK) 1300C************************************************************** 1301C 1302C Calculation of CCMM T^gB contribution to left and right 1303C Jacobian transformation: 1304C 1305C <mu|exp(-T)T^g|CC> for LR = 0 1306C F transformation for LR = F 1307C P transformation for LR = P 1308C 1309C LR = 'L','R','0','F','P' 1310C 1311C************************************************************** 1312C 1313#include "implicit.h" 1314#include "maxorb.h" 1315#include "mxcent.h" 1316#include "dummy.h" 1317#include "iratdef.h" 1318#include "priunit.h" 1319#include "ccorb.h" 1320#include "ccsdsym.h" 1321#include "ccsdio.h" 1322#include "ccsdinp.h" 1323#include "ccfield.h" 1324#include "exeinf.h" 1325#include "ccfdgeo.h" 1326#include "qm3.h" 1327#include "ccinftap.h" 1328#include "ccslvinf.h" 1329C 1330 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, IZERO = 0 , TWO = 2.0D0) 1331 PARAMETER (HALF = 0.5D0) 1332C 1333 DIMENSION WORK(LWORK),RHO1(*),RHO2(*),CTR1(*),CTR2(*) 1334 DIMENSION ALPHAMM(MXQM3) 1335C 1336 CHARACTER*8 LABEL,LABEL1,LIST*(2),LR*(1) 1337 CHARACTER*(8) FILMME, FILMMX 1338 CHARACTER*9 FILMMR 1339 LOGICAL LEXIST 1340C 1341C---------------------------------------------- 1342C If all models are SPC and OLDTG = .FALSE. 1343C -> RETURN from CCMM_LTRB: 1344C---------------------------------------------- 1345C 1346 IF ( (.NOT. (OLDTG)) .AND. (LOSPC) ) RETURN 1347C 1348C----------------------------------------------------- 1349C Also, if all models are SPC and OLDTG are .TRUE. 1350C only enter CCMM_LTRB if LR .NE.'0': 1351C----------------------------------------------------- 1352C 1353 IF ( (LOSPC) .AND. (OLDTG) .AND. (LR .NE.'0') ) RETURN 1354C 1355 IF (IQM3PR .GT. 10) THEN 1356 WRITE(LUPRI,*)'CCMM_LTRB: CC/MM contribution to CC L. transf.' 1357 WRITE(LUPRI,*)'CCMM_LTRB: LWORK:', LWORK 1358 WRITE(LUPRI,*)'CCMM_LTRB: LR:', LR 1359 ENDIF 1360C 1361 IF ( DEBUG .OR.(IQM3PR .GT. 10)) THEN 1362 RHO1N = DDOT(NT1AM(ISYMTR),RHO1,1,RHO1,1) 1363 RHO2N = DDOT(NT2AM(ISYMTR),RHO2,1,RHO2,1) 1364 WRITE(LUPRI,*) ' Norm of RHO1 in CCMM_LTGB on input:', RHO1N 1365 WRITE(LUPRI,*) ' Norm of RHO2 in CCMM_LTGB on input:', RHO2N 1366 RHO1N = DDOT(NT1AM(ISYMTR),CTR1,1,CTR1,1) 1367 RHO2N = DDOT(NT2AM(ISYMTR),CTR2,1,CTR2,1) 1368 WRITE(LUPRI,*) ' Norm af C1AM in CCMM_LTGB on input:', RHO1N 1369 WRITE(LUPRI,*) ' Norm af C2AM in CCMM_LTGB on input:', RHO2N 1370 ENDIF 1371C 1372C--------------------- 1373C Init parameters: 1374C--------------------- 1375C 1376 NTAMP1 = NT1AM(ISYMTR) 1377 NTAMP2 = NT2AM(ISYMTR) 1378 NTAMP = NTAMP1 + NTAMP2 1379C 1380 IF (CCS) NT2AM(ISYMTR) = IZERO 1381C------------------------------ 1382C 1383 IF (DISCEX) THEN 1384 LUMMET = -1 1385 LUMMXI = -1 1386 FILMME = 'CCMM_ETA' 1387 FILMMX = 'CCMM__XI' 1388 CALL WOPEN2(LUMMET, FILMME, 64, 0) 1389 CALL WOPEN2(LUMMXI, FILMMX, 64, 0) 1390 END IF 1391C 1392C----------------------------------------- 1393C 1394 KTGB = 1 1395 KWRK1 = KTGB + N2BST(ISYMTR) 1396 LWRK1 = LWORK - KWRK1 1397C 1398 IF (LWRK1 .LT. 0) CALL QUIT( ' Too litle work in CCMM_LTRB 1' ) 1399C 1400 CALL DZERO(WORK(KTGB),N2BST(ISYMTR)) 1401C 1402C-------------------------------------------------------- 1403C If OLDTG, LOSPC = .TRUE. and (LR .EQ.'0') only the 1404C contribution due to the partial charges in the T^g 1405C operator -> GOTO 888 (after the alpha part). 1406C-------------------------------------------------------- 1407C 1408 IF ( (LOSPC) .AND. (OLDTG) .AND. (LR .EQ.'0') ) GOTO 888 1409C 1410C-------------------------------------------------------- 1411C Calculate T^{gB}= - alpha sum_a <B|Rra|CC>Rra 1412C and contributions from this term. For LR = 0 1413C calculates <mu|T^g|CC>. 1414C-------------------------------------------------------- 1415C 1416C-------------------------- 1417C Dynamical allocation: 1418C-------------------------- 1419C 1420 KRAAOx = KWRK1 1421 KRAAOy = KRAAOx + N2BST(ISYMTR) 1422 KRAAOz = KRAAOy + N2BST(ISYMTR) 1423 KENSAx = KRAAOz + N2BST(ISYMTR) 1424 KENSAy = KENSAx + NCOMS 1425 KENSAz = KENSAy + NCOMS 1426 KRAx = KENSAz + NCOMS 1427 KRAy = KRAx + NCOMS 1428 KRAz = KRAy + NCOMS 1429 KWRK2 = KRAz + NCOMS 1430 LWRK2 = LWORK - KWRK2 1431C 1432 IF (LWRK2.LT.0) CALL QUIT( 'Too little work in CCMM_LTRB, 2') 1433C 1434 KXIx = KWRK2 1435 KXIy = KXIx + NTAMP 1436 KXIz = KXIy + NTAMP 1437 KWRK3 = KXIz + NTAMP 1438 LWRK3 = LWORK - KWRK3 1439C 1440 IF (LWRK3.LT.0) CALL QUIT( 'Too little work in CCMM_LTRB, 3') 1441C 1442 CALL DZERO(WORK(KRAAOx),3*N2BST(ISYMTR)) 1443 CALL DZERO(WORK(KENSAx),3*NCOMS) 1444 CALL DZERO(WORK(KRAx),3*NCOMS) 1445 CALL DZERO(WORK(KXIx),3*NTAMP) 1446C 1447 CALL CC_GET31('ENSAFILE',NCOMS, 1448 * WORK(KENSAx),WORK(KENSAy),WORK(KENSAz)) 1449 CALL CC_GET31('CC_RA',NCOMS, 1450 * WORK(KRAx),WORK(KRAy),WORK(KRAz)) 1451C 1452 LUCCEF = -1 1453 CALL GPOPEN(LUCCEF,'ELFDMM','OLD',' ', 1454 & 'UNFORMATTED',IDUMMY,.FALSE.) 1455 REWIND(LUCCEF) 1456C 1457C---------------------------- 1458C Readin integrals again, 1459C Read Rra in AO basis: 1460C---------------------------- 1461C 1462 LM = 0 1463 IADR1 = 1 1464 IADR2 = 1 1465 LEN = NTAMP 1466C 1467 DO 700 I = 1, ISYTP 1468 IF (MDLWRD(I)(1:5) .EQ. 'SPC_E') THEN 1469 DO 701 J = NSYSBG(I), NSYSED(I) 1470 DO 702 K = 1, NUALIS(I) 1471C 1472 LM = LM + 1 1473 LABEL = 'READ INT' 1474C 1475 CALL CCMMINT(LABEL,LUCCEF,WORK(KRAAOx),DUMMY,IRREP, 1476 & ISYM,IERR,WORK(KWRK3),LWRK3) 1477 CALL CCMMINT(LABEL,LUCCEF,WORK(KRAAOy),DUMMY,IRREP, 1478 & ISYM,IERR,WORK(KWRK3),LWRK3) 1479 CALL CCMMINT(LABEL,LUCCEF,WORK(KRAAOz),DUMMY,IRREP, 1480 & ISYM,IERR,WORK(KWRK3),LWRK3) 1481C 1482C------------------------------------------------------------ 1483C Calculate xi^Rra (for LR=R actually eta^Rra ) 1484C 1485C Only if MDLWRD(I) = 'SPC_E01' are we going to calculate 1486C contributions from the T^g operator to the response 1487C equations. 1488C------------------------------------------------------------ 1489C 1490 IF (MDLWRD(I) .EQ. 'SPC_E01') THEN 1491 IF (.NOT. (DISCEX)) THEN 1492C 1493 LABEL = 'GIVE INT' 1494 IF ( (LR.EQ.'L').OR.(LR.EQ.'0') 1495 & .OR.(LR.EQ.'P') ) THEN 1496 CALL CC_XKSI(WORK(KXIx), 1497 * LABEL,ISYMTR,0,WORK(KRAAOx), 1498 * WORK(KWRK3),LWRK3) 1499 ELSE IF ((LR.EQ.'R').OR.(LR.EQ.'F')) THEN 1500 LIST = 'L0' 1501 ILSTNR = 1 1502 CALL CC_ETAC(ISYMTR,LABEL,WORK(KXIx), 1503 * LIST,ILSTNR,0,WORK(KRAAOx), 1504 * WORK(KWRK3),LWRK3) 1505 END IF 1506C 1507 IF ( (LR.EQ.'L').OR.(LR.EQ.'0') 1508 & .OR.(LR.EQ.'P') ) THEN 1509 CALL CC_XKSI(WORK(KXIy),LABEL,ISYMTR, 1510 * 0,WORK(KRAAOy), 1511 * WORK(KWRK3),LWRK3) 1512 ELSE IF ((LR.EQ.'R').OR.(LR.EQ.'F')) THEN 1513 LIST = 'L0' 1514 ILSTNR = 1 1515 CALL CC_ETAC(ISYMTR,LABEL,WORK(KXIy), 1516 * LIST,ILSTNR,0,WORK(KRAAOy), 1517 * WORK(KWRK3),LWRK3) 1518 END IF 1519C 1520 IF ( (LR.EQ.'L').OR.(LR.EQ.'0') 1521 & .OR.(LR.EQ.'P') ) THEN 1522 CALL CC_XKSI(WORK(KXIz),LABEL,ISYMTR, 1523 * 0,WORK(KRAAOz), 1524 * WORK(KWRK3),LWRK3) 1525 ELSE IF ((LR.EQ.'R').OR.(LR.EQ.'F')) THEN 1526 LIST = 'L0' 1527 ILSTNR = 1 1528 CALL CC_ETAC(ISYMTR,LABEL,WORK(KXIz), 1529 * LIST,ILSTNR,0,WORK(KRAAOz),WORK(KWRK3),LWRK3) 1530 END IF 1531C 1532 ELSE 1533C 1534 IF ( (LR.EQ.'L').OR.(LR.EQ.'P') ) THEN 1535 CALL GETWA2(LUMMXI,FILMMX,WORK(KXIx),IADR2,LEN) 1536 IADR2 = IADR2 + LEN 1537 ELSE IF ( (LR.EQ.'R').OR.(LR.EQ.'F') ) THEN 1538 CALL GETWA2(LUMMET,FILMME,WORK(KXIx),IADR1,LEN) 1539 IADR1 = IADR1 + LEN 1540 ELSE IF (LR.EQ.'0') THEN 1541 LABEL = 'GIVE INT' 1542 CALL CC_XKSI(WORK(KXIx),LABEL,ISYMTR,0,WORK(KRAAOx), 1543 * WORK(KWRK3),LWRK3) 1544 END IF 1545C 1546 IF ( (LR.EQ.'L').OR.(LR.EQ.'P') ) THEN 1547 CALL GETWA2(LUMMXI,FILMMX,WORK(KXIy),IADR2,LEN) 1548 IADR2 = IADR2 + LEN 1549 ELSE IF ( (LR.EQ.'R').OR.(LR.EQ.'F') ) THEN 1550 CALL GETWA2(LUMMET,FILMME,WORK(KXIy),IADR1,LEN) 1551 IADR1 = IADR1 + LEN 1552 ELSE IF (LR.EQ.'0') THEN 1553 LABEL = 'GIVE INT' 1554 CALL CC_XKSI(WORK(KXIy),LABEL,ISYMTR,0,WORK(KRAAOy), 1555 * WORK(KWRK3),LWRK3) 1556 END IF 1557C 1558 IF ( (LR.EQ.'L').OR.(LR.EQ.'P') ) THEN 1559 CALL GETWA2(LUMMXI,FILMMX,WORK(KXIz),IADR2,LEN) 1560 IADR2 = IADR2 + LEN 1561 ELSE IF ((LR.EQ.'R').OR.(LR.EQ.'F')) THEN 1562 CALL GETWA2(LUMMET,FILMME,WORK(KXIz),IADR1,LEN) 1563 IADR1 = IADR1 + LEN 1564 ELSE IF (LR.EQ.'0') THEN 1565 LABEL = 'GIVE INT' 1566 CALL CC_XKSI(WORK(KXIz),LABEL,ISYMTR,0,WORK(KRAAOz), 1567 * WORK(KWRK3),LWRK3) 1568 END IF 1569 END IF 1570C 1571C------------------------------- 1572C Contract with B: 1573C------------------------------- 1574C 1575 IF (LR.NE.'0') THEN 1576C 1577 KXI1x = KXIx 1578 KXI2x = KXIx + NTAMP1 1579 BXILMD1x= DDOT(NTAMP1,CTR1,1,WORK(KXI1x),1) 1580 BXILMD2x= DDOT(NTAMP2,CTR2,1,WORK(KXI2x),1) 1581 BXILMDx = BXILMD1x + BXILMD2x 1582C 1583 KXI1y = KXIy 1584 KXI2y = KXIy + NTAMP1 1585 BXILMD1y = DDOT(NTAMP1,CTR1,1,WORK(KXI1y),1) 1586 BXILMD2y = DDOT(NTAMP2,CTR2,1,WORK(KXI2y),1) 1587 BXILMDy = BXILMD1y + BXILMD2y 1588C 1589 KXI1z = KXIz 1590 KXI2z = KXIz + NTAMP1 1591 BXILMD1z = DDOT(NTAMP1,CTR1,1,WORK(KXI1z),1) 1592 BXILMD2z = DDOT(NTAMP2,CTR2,1,WORK(KXI2z),1) 1593 BXILMDz = BXILMD1z + BXILMD2z 1594C 1595 1596 FACx = - ALPIMM(I,K) * (BXILMDx) 1597 FACy = - ALPIMM(I,K) * (BXILMDy) 1598 FACz = - ALPIMM(I,K) * (BXILMDz) 1599C 1600 IF (IQM3PR .GT. 5) THEN 1601 WRITE(LUPRI,*)'---------------------------' // 1602 * '-------------------' 1603 WRITE(LUPRI,*)'The factors calculated at cent.' // 1604 * ' of mass no.',LM 1605 WRITE(LUPRI,*)'---------------------------' // 1606 * '-------------------' 1607 WRITE(LUPRI,*)'FACx=',FACx 1608 WRITE(LUPRI,*)'FACy=',FACy 1609 WRITE(LUPRI,*)'FACz=',FACz 1610 WRITE(LUPRI,*)'---------------------------' // 1611 * '-------------------' 1612 END IF 1613C 1614C------------------------------------------ 1615C Add to T^{qB} integrals. 1616C (for LR=R actually T^gC): 1617C------------------------------------------ 1618C 1619 CALL DAXPY(N2BST(ISYMTR),FACx,WORK(KRAAOx),1, 1620 * WORK(KTGB),1) 1621C 1622 CALL DAXPY(N2BST(ISYMTR),FACy,WORK(KRAAOy),1, 1623 * WORK(KTGB),1) 1624C 1625 CALL DAXPY(N2BST(ISYMTR),FACz,WORK(KRAAOz),1, 1626 * WORK(KTGB),1) 1627C 1628 END IF 1629 IF (LR.EQ.'0') THEN 1630C 1631 KXI1x = KXIx 1632 KXI2x = KXIx + NTAMP1 1633C 1634 FACx = - ALPIMM(I,K) * (WORK(KRAx + LM - 1) 1635 & + 0.5D0 * WORK(KENSAx + LM - 1)) 1636C 1637 CALL DAXPY(NT1AM(ISYMTR),FACx,WORK(KXI1x),1, 1638 * RHO1,1) 1639 CALL DAXPY(NT2AM(ISYMTR),FACx,WORK(KXI2x),1, 1640 * RHO2,1) 1641C 1642 KXI1y= KXIy 1643 KXI2y = KXIy + NTAMP1 1644C 1645 FACy = - ALPIMM(I,K) * (WORK(KRAy + LM - 1) 1646 & + 0.5D0 * WORK(KENSAy + LM - 1)) 1647C 1648 CALL DAXPY(NT1AM(ISYMTR),FACy,WORK(KXI1y),1, 1649 * RHO1,1) 1650 CALL DAXPY(NT2AM(ISYMTR),FACy,WORK(KXI2y),1, 1651 * RHO2,1) 1652C 1653 KXI1z = KXIz 1654 KXI2Z = KXIz + NTAMP1 1655C 1656 FACz = - ALPIMM(I,K) * (WORK(KRAz + LM - 1) 1657 & + 0.5D0 * WORK(KENSAz + LM - 1)) 1658C 1659 CALL DAXPY(NT1AM(ISYMTR),FACz,WORK(KXI1z),1, 1660 * RHO1,1) 1661 CALL DAXPY(NT2AM(ISYMTR),FACz,WORK(KXI2z),1, 1662 * RHO2,1) 1663C 1664 END IF 1665 END IF 1666 702 CONTINUE 1667 701 CONTINUE 1668 END IF 1669 700 CONTINUE 1670C 1671 CALL GPCLOSE(LUCCEF,'KEEP') 1672C 1673 IF (DISCEX) THEN 1674 CALL WCLOSE2(LUMMET, FILMME, 'KEEP') 1675 CALL WCLOSE2(LUMMXI, FILMMX, 'KEEP') 1676 END IF 1677C 1678C----------------------- 1679C End of alpha part: 1680C----------------------- 1681C 1682 888 CONTINUE 1683C 1684C------------------------------ 1685C If OLDTG = .TRUE. 1686C N_s part -> T^g operator: 1687C------------------------------ 1688C 1689 IF (OLDTG) THEN 1690 IF (LR.EQ.'0') THEN 1691C 1692 KNSAO = KWRK1 1693 KWRK2 = KNSAO + N2BST(ISYMTR) 1694 LWRK2 = LWORK - KWRK2 1695C 1696 CALL DZERO(WORK(KNSAO),N2BST(ISYMTR)) 1697C 1698 IF (LWRK2 .LT. 0) 1699 & CALL QUIT( ' Too litle work in CCMM_LTRB 4b' ) 1700C 1701 KXI = KWRK2 1702 KWRK3 = KXI + NTAMP 1703 LWRK3 = LWORK - KWRK3 1704C 1705 CALL DZERO(WORK(KXI),NTAMP) 1706C 1707 IF (LWRK3.LT.0) 1708 & CALL QUIT( 'Too little work in CCMM_LTRB, 5') 1709C 1710 LUCCPO = -1 1711 CALL GPOPEN(LUCCPO,'POTMM','OLD',' ', 1712 & 'UNFORMATTED',IDUMMY,.FALSE.) 1713 REWIND (LUCCPO) 1714C 1715 FAC1 = -1.0D0 1716 LM = 0 1717C 1718 DO 800 I = 1, ISYTP 1719 IF (MDLWRD(I)(1:3) .EQ. 'SPC') THEN 1720 DO 801 J = NSYSBG(I), NSYSED(I) 1721 DO 802 K = 1,NSISY(I) 1722C 1723 LM = LM +1 1724 LABEL = 'READ INT' 1725C 1726 CALL CCMMINT(LABEL,LUCCPO,WORK(KNSAO),DUMMY,IRREP, 1727 & ISYM,IERR,WORK(KWRK3),LWRK3) 1728C 1729C--------------------------------- 1730C Calculate xi^Ns: 1731C--------------------------------- 1732C 1733 LABEL = 'GIVE INT' 1734C 1735 CALL CC_XKSI(WORK(KXI),LABEL,ISYMTR,0, 1736 * WORK(KNSAO),WORK(KWRK3),LWRK3) 1737C 1738 KXI1 = KXI 1739 KXI2 = KXI + NTAMP1 1740C 1741 CALL DAXPY(NT1AM(ISYMTR),FAC1,WORK(KXI1),1, 1742 * RHO1,1) 1743 CALL DAXPY(NT2AM(ISYMTR),FAC1,WORK(KXI2),1, 1744 * RHO2,1) 1745 802 CONTINUE 1746 801 CONTINUE 1747 END IF 1748 800 CONTINUE 1749 IF ( DEBUG .OR.(IQM3PR .GT. 10)) THEN 1750 RHO1N = DDOT(NT1AM(ISYMTR),WORK(KXI1),1,WORK(KXI1),1) 1751 RHO2N = DDOT(NT2AM(ISYMTR),WORK(KXI2),1,WORK(KXI2),1) 1752 WRITE(LUPRI,*) ' Norm af elstat contribution to LHTR1:', RHO1N 1753 WRITE(LUPRI,*) ' Norm af elstat contribution to LHTR2:', RHO2N 1754 END IF 1755C 1756 CALL GPCLOSE(LUCCPO,'KEEP') 1757C 1758 END IF 1759 END IF 1760C 1761C------------------------------------------------------- 1762C 1763C Introduce repulsion component into T^g 1764C 1765C------------------------------------------------------- 1766C 1767 IF (LOSHAW) THEN 1768 IF (LR.EQ.'0') THEN 1769C 1770 KINT = KWRK1 1771 KWRK2 = KINT + N2BST(ISYMOP) 1772 LWRK2 = LWORK - KWRK2 1773C 1774 IF (LWRK2 .LT. 0) THEN 1775 CALL QUIT( 'Too litle work in CCMM_LTRB Rep. 1') 1776 END IF 1777C 1778 CALL DZERO(WORK(KINT),N2BST(ISYMOP)) 1779C 1780 FILMMR = 'MMREPOP_0' 1781 LUMMRE = -1 1782 IADRFIL = 1 1783 NDATA = N2BST(ISYMOP) 1784C 1785 CALL WOPEN2(LUMMRE,FILMMR,64,0) 1786 CALL GETWA2(LUMMRE,FILMMR,WORK(KINT),IADRFIL,NDATA) 1787C 1788 IF (REPTST) THEN 1789 NBASLO = MAX(NBAST,1) 1790C 1791 KINT1 = KWRK2 1792 KINT2 = KINT1 + N2BST(ISYMOP) 1793 KWRK3 = KINT2 + N2BST(ISYMOP) 1794 LWRK3 = LWORK - KWRK3 1795C 1796 IF (LWRK3.LT.0) THEN 1797 CALL QUIT( 'Too little work in CCMM_LTRB Rep. 1a') 1798 END IF 1799C 1800 CALL DZERO(WORK(KINT1),2*N2BST(ISYMOP)) 1801 CALL DCOPY(N2BST(ISYMOP),WORK(KINT),1,WORK(KINT1),1) 1802C 1803 DO 668 KR=1,NREPMT 1804 CALL DGEMM('N','N',NBAST,NBAST,NBAST, 1805 * ONE,WORK(KINT),NBASLO, 1806 * WORK(KINT1),NBASLO, 1807 * ZERO,WORK(KINT2),NBASLO) 1808C 1809 CALL DCOPY(N2BST(ISYMOP),WORK(KINT2),1,WORK(KINT),1) 1810 668 CONTINUE 1811 END IF 1812C 1813 IF (IQM3PR .GT.14) THEN 1814 TAL1 = DDOT(N2BST(ISYMOP),WORK(KINT),1,WORK(KINT),1) 1815 TAL2 = SQRT(TAL1) 1816 WRITE(lupri,*)'Norm of rep. operator:',TAL2 1817 END IF 1818C 1819 KXI = KWRK2 1820 KWRK3 = KXI + NTAMP 1821 LWRK3 = LWORK - KWRK3 1822C 1823 IF (LWRK3.LT.0) THEN 1824 CALL QUIT( 'Too little work in CCMM_LTRB Rep. 2') 1825 END IF 1826C 1827 CALL DZERO(WORK(KXI),NTAMP) 1828C 1829 LABEL = 'GIVE INT' 1830 CALL CC_XKSI(WORK(KXI),LABEL,ISYMTR,0, 1831 * WORK(KINT),WORK(KWRK3),LWRK3) 1832C 1833 KXI1 = KXI 1834 KXI2 = KXI + NTAMP1 1835C 1836 TAL3 = DDOT(NT1AM(ISYMTR),WORK(KXI1),1,WORK(KXI1),1) 1837 TAL4 = DDOT(NT2AM(ISYMTR),WORK(KXI2),1,WORK(KXI2),1) 1838C 1839 FAC2 = 1.0D0 1840 CALL DAXPY(NT1AM(ISYMTR),FAC2,WORK(KXI1),1, 1841 * RHO1,1) 1842 CALL DAXPY(NT2AM(ISYMTR),FAC2,WORK(KXI2),1, 1843 * RHO2,1) 1844C 1845 CALL WCLOSE2(LUMMRE,FILMMR, 'KEEP') 1846C 1847C test test test 1848 END IF 1849C test slut 1850 END IF 1851C 1852C--------------------------------------------------- 1853C Calculate contribution from the T^gB operator: 1854C--------------------------------------------------- 1855C 1856 IF (LR.NE.'0') THEN 1857C 1858 KETA = KWRK1 1859 KWRK2 = KETA + NTAMP 1860 LWRK2 = LWORK - KWRK2 1861C 1862 IF (LWRK2 .LT. 0) 1863 & CALL QUIT( ' Too litle work in CCMM_LTRB 6' ) 1864C 1865 CALL DZERO(WORK(KETA),NTAMP) 1866C 1867 KETA1 = KETA 1868 KETA2 = KETA + NTAMP1 1869C 1870 IF ((LR.EQ.'L').OR.(LR.EQ.'F')) THEN 1871C 1872 LIST = 'L0' 1873 LABEL = 'GIVE INT' 1874C 1875 CALL CC_ETAC(ISYMTR,LABEL,WORK(KETA), 1876 * LIST,1,0,WORK(KTGB),WORK(KWRK2),LWRK2) 1877 ELSE IF ((LR.EQ.'R').OR.(LR.EQ.'P')) THEN 1878C 1879 LABEL = 'GIVE INT' 1880C 1881 CALL CC_XKSI(WORK(KETA),LABEL,ISYMTR,0,WORK(KTGB), 1882 * WORK(KWRK2),LWRK2) 1883 IF (LR.EQ.'R') 1884 * CALL CCLR_DIASCL(WORK(KETA2),TWO,ISYMTR) 1885 END IF 1886C 1887 IF ( DEBUG .OR.(IQM3PR .GT. 10)) THEN 1888 RHO1N = DDOT(NT1AM(ISYMTR),WORK(KETA1),1,WORK(KETA1),1) 1889 RHO2N = DDOT(NT2AM(ISYMTR),WORK(KETA2),1,WORK(KETA2),1) 1890 WRITE(LUPRI,*) ' Norm af T^gB contribution to LHTR1:', RHO1N 1891 WRITE(LUPRI,*) ' Norm af T^gB contribution to LHTR2:', RHO2N 1892 END IF 1893C 1894 CALL DAXPY(NT1AM(ISYMTR),ONE,WORK(KETA1),1,RHO1,1) 1895 CALL DAXPY(NT2AM(ISYMTR),ONE,WORK(KETA2),1,RHO2,1) 1896C 1897 IF ( DEBUG .OR.(IQM3PR .GT. 10)) THEN 1898 RHO1N = DDOT(NT1AM(ISYMTR),RHO1,1,RHO1,1) 1899 RHO2N = DDOT(NT2AM(ISYMTR),RHO2,1,RHO2,1) 1900 WRITE(LUPRI,*) ' Norm af RHO1 in CCMM_LTGB:', RHO1N 1901 WRITE(LUPRI,*) ' Norm af RHO2 in CCMM_LTGB:', RHO2N 1902 END IF 1903 END IF 1904C 1905 END 1906C 1907C******************************************************************** 1908C /* Deck ccmmint */ 1909 SUBROUTINE CCMMINT(LABEL,LU,RAAOx,Xint,IRREP, 1910 * ISYM,IERR,WORK,LWORK) 1911C******************************************************************** 1912C 1913C Purpose: read property one-electron AO integrals from file or copy 1914C from array and transform to CC format. 1915C 1916C input: LU -- Logical unit number for file to be read. 1917C 1918C output: property AO integrals in coupled cluster 1919C storage scheme (dimension of ouput vector 1920C is N2BST(IRREP) 1921C IRREP -- irrep symmetry of the operator 1922C ISYM -- +1 for a symmetric operator 1923C -1 for an antisymmetric operator 1924C 0 if unknown 1925C 1926C******************************************************************** 1927C 1928#if defined (IMPLICIT_NONE) 1929 IMPLICIT NONE 1930#else 1931# include "implicit.h" 1932#endif 1933#include "priunit.h" 1934#include "ccsdinp.h" 1935#include "ccsdsym.h" 1936#include "maxorb.h" 1937#include "ccorb.h" 1938#include "inftap.h" 1939#include "ccisao.h" 1940#include "dummy.h" 1941#include "qm3.h" 1942C 1943C local parameters: 1944C 1945 LOGICAL LOCDBG 1946 PARAMETER (LOCDBG = .FALSE.) 1947C 1948#if defined (SYS_CRAY) 1949 REAL CKMXPR 1950#else 1951 DOUBLE PRECISION CKMXPR 1952#endif 1953 PARAMETER (CKMXPR = 1.0d-12) 1954C 1955C input: 1956C 1957 CHARACTER*8 LABEL 1958 CHARACTER*8 RRAO 1959 INTEGER IRREP, ISYM, IERR, LWORK, LU 1960C 1961#if defined (SYS_CRAY) 1962 REAL WORK(LWORK) 1963 REAL RAAOx(*),Xint(*) 1964#else 1965 DOUBLE PRECISION WORK(LWORK) 1966 DOUBLE PRECISION RAAOx(*),Xint(*) 1967#endif 1968C 1969 LOGICAL LOPENED 1970 INTEGER IDX, IDXI, IDXJ 1971 INTEGER KEND0, LEND0, KINTEG, ISYMA, ISYMB, IOLD, INEW 1972C 1973C functions: 1974C 1975 INTEGER IDAMAX 1976C 1977C-------------------------- 1978C Dynamical allocation: 1979C-------------------------- 1980C 1981 KINTEG = 1 1982 KEND0 = KINTEG + N2BASX 1983 LEND0 = LWORK - KEND0 1984C 1985 IF (LEND0 .LT. 0) CALL QUIT('Insufficient memory in CCRRA.') 1986C 1987 CALL DZERO(WORK(KINTEG),N2BASX) 1988C 1989 ISYM = +1 1990C 1991 IF (LABEL.EQ.'GIVE INT') THEN 1992 IF (IQM3PR .GT. 15) WRITE (LUPRI,*) ' Starting to copy!' 1993 CALL DCOPY(NNBASX,XINT,1,WORK(KEND0),1) 1994 ELSE 1995 CALL READT(LU,NNBASX,WORK(KEND0)) 1996 END IF 1997C 1998 CALL DSPTGE(NBAST,WORK(KEND0),WORK(KINTEG)) 1999C 2000 IF (IQM3PR .GT. 99 .OR. LOCDBG) THEN 2001 CALL AROUND('AOQM3INT: -integrals') 2002 CALL OUTPUT(WORK(KINTEG),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI) 2003 END IF 2004C 2005C------------------------------------- 2006C Find irrep symmetry of the operator: 2007C------------------------------------- 2008C 2009 IDX = IDAMAX(N2BASX,WORK(KINTEG),1) 2010C 2011 IF ( ABS(WORK(KINTEG-1+IDX)) .GT. CKMXPR ) THEN 2012 IDXI = (IDX+NBAST-1) / NBAST 2013 IDXJ = IDX - (IDXI-1)*NBAST 2014 IRREP = MULD2H(ISAO(IDXI),ISAO(IDXJ)) 2015 ELSE 2016 IRREP = 0 2017 IERR = -1 2018C 2019 IF (IQM3PR .GT. 5) THEN 2020 WRITE (LUPRI,*) 2021 & 'WARNING: irrep symmetry can not be determined.' 2022 WRITE (LUPRI,*) 2023 & 'Irrep set to 1!!!! ' 2024 END IF 2025C 2026 IRREP = 1 2027 CALL FLSHFO(LUPRI) 2028 END IF 2029C 2030C--------------------------------------------- 2031C Resort integrals to coupled cluster storage: 2032C--------------------------------------------- 2033C 2034 DO ISYMA = 1, NSYM 2035 ISYMB = MULD2H(IRREP,ISYMA) 2036 DO A = 1, NBAS(ISYMA) 2037 DO B = 1, NBAS(ISYMB) 2038 IOLD = (IBAS(ISYMA)+A-1)*NBAST + (IBAS(ISYMB)+B) 2039 INEW = IAODIS(ISYMB,ISYMA) + NBAS(ISYMB)*(A-1) + B 2040 RAAOx(INEW) = WORK(KINTEG-1+IOLD) 2041 END DO 2042 END DO 2043 END DO 2044C 2045 RETURN 2046 END 2047C 2048C************************************************************** 2049C /* Deck hfint */ 2050 SUBROUTINE HFINT(EINTF,WORK,LWORK) 2051C************************************************************** 2052C 2053C This routine calculates the electrostatic HF/MM electronic 2054C interaction energy. Note that the interaction energy 2055C is calculated using the wave function optimized with 2056C the induced dipoles determined at the CC level of theory. 2057C Thus, this is only a pseudo HF result. 2058C 2059C************************************************************** 2060C 2061#include "implicit.h" 2062#include "maxorb.h" 2063#include "mxcent.h" 2064#include "nuclei.h" 2065#include "dummy.h" 2066#include "iratdef.h" 2067#include "priunit.h" 2068#include "ccorb.h" 2069#include "ccsdsym.h" 2070#include "ccsdio.h" 2071#include "ccsdinp.h" 2072#include "ccfield.h" 2073#include "exeinf.h" 2074#include "ccfdgeo.h" 2075#include "ccinftap.h" 2076#include "qm3.h" 2077C 2078 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, IZERO = 0 , TWO = 2.0D0) 2079C 2080 DIMENSION WORK(LWORK) 2081C 2082 CHARACTER*8 LABEL 2083C 2084C-------------------------- 2085C Dynamical allocation: 2086C-------------------------- 2087C 2088 KDENS = 1 2089 KNSAO = KDENS + N2BST(ISYMOP) 2090 KHFNS = KNSAO + N2BST(ISYMOP) 2091 KWRK2 = KHFNS + NUSITE 2092 LWRK2 = LWORK - KWRK2 2093C 2094 IF (LWRK2 .LT. 0) THEN 2095 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KWRK2 2096 CALL QUIT('Insufficient memory for CCS AO-density in '// 2097 & 'HFINT') 2098 ENDIF 2099C 2100 CALL DZERO(WORK(KDENS),N2BST(ISYMOP)) 2101 CALL DZERO(WORK(KNSAO),N2BST(ISYMOP)) 2102 CALL DZERO(WORK(KHFNS),NUSITE) 2103C 2104C------------------------------------------ 2105C First we calculate the HF AO density: 2106C------------------------------------------ 2107C 2108 CALL CCS_D1AO(WORK(KDENS),WORK(KWRK2),LWRK2) 2109C 2110 IF (IQM3PR .GT. 50) THEN 2111 CALL AROUND('CCS One electron density in HFINT') 2112 CALL CC_PRFCKAO(WORK(KDENS),1) 2113 ENDIF 2114C 2115C----------------------------------------- 2116C Read integrals from file 2117C and perform the calculation of HFNS: 2118C----------------------------------------- 2119C 2120 LUCCPO = -1 2121 CALL GPOPEN(LUCCPO,'POTMM','UNKNOWN',' ', 2122 & 'UNFORMATTED',IDUMMY,.FALSE.) 2123 REWIND(LUCCPO) 2124C 2125 LM = 0 2126 DO 950 I = 1, ISYTP 2127 IF (MDLWRD(I)(1:3) .EQ. 'SPC') THEN 2128 DO 951 J = NSYSBG(I), NSYSED(I) 2129 DO 952 K = 1,NSISY(I) 2130 LM = LM +1 2131 LABEL = 'READ INT' 2132 CALL CCMMINT(LABEL,LUCCPO,WORK(KNSAO),DUMMY, 2133 * IRREP, ISYM,IERR, 2134 * WORK(KWRK2),LWRK2) 2135 WORK(KHFNS+LM-1) = DDOT(N2BST(ISYMOP), 2136 * WORK(KNSAO),1, 2137 * WORK(KDENS),1) 2138 952 CONTINUE 2139 951 CONTINUE 2140 END IF 2141 950 CONTINUE 2142C 2143 CALL GPCLOSE(LUCCPO,'KEEP') 2144C 2145C------------------- 2146C Print section: 2147C------------------- 2148C 2149 IF (IQM3PR .GT. 5) THEN 2150 WRITE(LUPRI,'(//,A)') 2151 * ' +======================+' 2152 WRITE(LUPRI,'(A)') 2153 * ' |Site| <N_s> |' 2154 WRITE(LUPRI,'(1X,A)') 2155 * '+======================+' 2156 LS = 0 2157 DO 215 I = 1, ISYTP 2158 IF (MDLWRD(I)(1:3) .EQ. 'SPC') THEN 2159 DO 216 J = NSYSBG(I), NSYSED(I) 2160 DO 217 K = 1, NSISY(I) 2161 LS = LS + 1 2162 WRITE(LUPRI,'(1X,A,I3,A,F16.10,A)') 2163 * '| ', LS,'|', WORK(KHFNS+LS-1),' |' 2164 WRITE(LUPRI,'(1X,A)') 2165 * '+----------------------+' 2166 217 CONTINUE 2167 216 CONTINUE 2168 END IF 2169 215 CONTINUE 2170 WRITE(LUPRI,'(//,A)') 2171 END IF 2172C 2173C------------------------------- 2174C Adding each contribution : 2175C------------------------------- 2176C 2177 EINT = 0.0D0 2178 L = 0 2179C 2180 DO 130 I = 1, ISYTP 2181 IF (MDLWRD(I)(1:3) .EQ. 'SPC') THEN 2182 DO 140 J = NSYSBG(I), NSYSED(I) 2183 DO 150 K = 1, NSISY(I) 2184 L = L + 1 2185 EINT = EINT - WORK(KHFNS+L-1) 2186 150 CONTINUE 2187 140 CONTINUE 2188 END IF 2189 130 CONTINUE 2190C 2191 EINTF = EINT 2192C 2193 END 2194C 2195C************************************************************** 2196C /* Deck epert2 */ 2197 SUBROUTINE EPERT2(EPOL,WORK,LWORK) 2198C************************************************************** 2199C 2200#include "implicit.h" 2201#include "maxorb.h" 2202#include "mxcent.h" 2203#include "nuclei.h" 2204#include "dummy.h" 2205#include "iratdef.h" 2206#include "priunit.h" 2207#include "ccorb.h" 2208#include "ccsdsym.h" 2209#include "ccsdio.h" 2210#include "ccsdinp.h" 2211#include "ccfield.h" 2212#include "exeinf.h" 2213#include "ccfdgeo.h" 2214#include "ccinftap.h" 2215#include "qm3.h" 2216C 2217 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, IZERO = 0 , TWO = 2.0D0) 2218C 2219 DIMENSION WORK(LWORK) 2220 DIMENSION EELEC(3,MXQM3) 2221 DIMENSION FFJ(3) 2222 LOGICAL LOINDM 2223C 2224 DO 879 I = 1, MXQM3 2225 DO 880 J = 1, 3 2226 EELEC(J,I) = 0.0D0 2227 880 CONTINUE 2228 879 CONTINUE 2229C 2230C-------------------------- 2231C Dynamical allocation: 2232C-------------------------- 2233C 2234 KRAx = 1 2235 KRAy = KRAx + NCOMS 2236 KRAz = KRAy + NCOMS 2237 KOMMSx = KRAz + NCOMS 2238 KOMMSy = KOMMSx + NCOMS 2239 KOMMSz = KOMMSy + NCOMS 2240 KWRK1 = KOMMSz + NCOMS 2241 LWRK1 = LWORK - KWRK1 2242C 2243 CALL DZERO(WORK(KRAx),6*NCOMS) 2244C 2245 IF (LWRK1.LT.0) CALL QUIT( 'Too little work space in CC_QM3') 2246C 2247 CALL CC_GET31('CC_RA',NCOMS, 2248 * WORK(KRAx),WORK(KRAy),WORK(KRAz)) 2249C 2250 DO 111 I = 1,NCOMS 2251 EELEC(1,I) = WORK(KRAx + I - 1) 2252 EELEC(2,I) = WORK(KRAy + I - 1) 2253 EELEC(3,I) = WORK(KRAz + I - 1) 2254 111 CONTINUE 2255C 2256 FFJ(1) = 0.0D0 2257 FFJ(2) = 0.0D0 2258 FFJ(3) = 0.0D0 2259C 2260 IF (RELMOM) THEN 2261 DO 330 IF =1, NFIELD 2262 IF (LFIELD(IF) .EQ. 'XDIPLEN') FFJ(1) = FFJ(1) + EFIELD(IF) 2263 IF (LFIELD(IF) .EQ. 'YDIPLEN') FFJ(2) = FFJ(2) + EFIELD(IF) 2264 IF (LFIELD(IF) .EQ. 'XDIPLEN') FFJ(3) = FFJ(3) + EFIELD(IF) 2265 330 CONTINUE 2266C 2267 WRITE(LUPRI,*)'STATIC ELECTRIC FIELD ADDED (x,y,z)' 2268 WRITE(LUPRI,*) FFJ(1),FFJ(2),FFJ(3) 2269 END IF 2270C 2271 LOINDM = .TRUE. 2272 CALL INDMOM(EELEC,LOINDM,FFJ) 2273 CALL QM3_OBAR(FFJ) 2274 CALL QM3_OTILDE(OTILDE,FFJ) 2275C 2276 CALL CC_GET31('OBARFILE',NCOMS, 2277 * WORK(KOMMSx),WORK(KOMMSy),WORK(KOMMSz)) 2278C 2279 EQM3T = ZERO 2280 KS = 0 2281C 2282 DO 110 I = 1, ISYTP 2283 IF (MDLWRD(I)(1:5) .EQ. 'SPC_E') THEN 2284 DO 120 J = NSYSBG(I), NSYSED(I) 2285 DO 134 K = 1, NUALIS(I) 2286 KS = KS+1 2287 RAT = ZERO 2288 RAT = 0.5D0 * ALPIMM(I,K) * WORK(KRAx + KS -1) * 2289 & (WORK(KRAx + KS -1) + WORK(KOMMSx + KS - 1)) 2290 & + 0.5D0 * ALPIMM(I,K) * WORK(KRAy + KS -1) * 2291 & (WORK(KRAy + KS -1) + WORK(KOMMSy + KS - 1)) 2292 & + 0.5D0 * ALPIMM(I,K) * WORK(KRAz + KS -1) * 2293 & (WORK(KRAz + KS -1) + WORK(KOMMSz + KS - 1)) 2294 EQM3T = EQM3T - RAT 2295 134 CONTINUE 2296 120 CONTINUE 2297 END IF 2298 110 CONTINUE 2299C 2300 IF (IQM3PR .GT. 10) THEN 2301 DO 777 I=1,NCOMS 2302 WRITE(LUPRI,*)'WORK(KOMMSx) =',WORK(KOMMSx + I - 1) 2303 WRITE(LUPRI,*)'WORK(KOMMSy) =',WORK(KOMMSy + I - 1) 2304 WRITE(LUPRI,*)'WORK(KOMMSz) =',WORK(KOMMSz + I - 1) 2305 777 CONTINUE 2306 END IF 2307 2308 EPOL = EQM3T 2309 EPOL = EPOL + OTILDE 2310C 2311 END 2312C 2313C************************************************************** 2314C /* Deck cc_put31 */ 2315 SUBROUTINE CC_PUT31(FLNAME,NULOOP,OMMSx,OMMSy,OMMSz) 2316C************************************************************** 2317C 2318#include "implicit.h" 2319#include "dummy.h" 2320C 2321 CHARACTER*(*) FLNAME 2322 INTEGER NMBU,NULOOP 2323 DIMENSION OMMSx(*) , OMMSy(*) , OMMSz(*) 2324C 2325 NMBU = -1 2326 CALL GPOPEN(NMBU,FLNAME,'UNKNOWN',' ','FORMATTED',IDUMMY,.FALSE.) 2327C 2328 REWIND (NMBU) 2329 LM = 0 2330C 2331 DO 820 L = 1,NULOOP 2332 LM = LM + 1 2333 WRITE(NMBU,'(I5,3E25.15)') LM,OMMSx(LM),OMMSy(LM),OMMSz(LM) 2334 820 CONTINUE 2335C 2336 REWIND (NMBU) 2337 CALL GPCLOSE(NMBU,'KEEP') 2338C 2339 END 2340C 2341C************************************************************** 2342C /* Deck cc_get31 */ 2343 SUBROUTINE CC_GET31(FLNAME,NULOOP,OMMSx,OMMSy,OMMSz) 2344C************************************************************** 2345C 2346#include "implicit.h" 2347#include "dummy.h" 2348#include "priunit.h" 2349C 2350C 2351 INTEGER NMBU,NULOOP 2352 CHARACTER*(*) FLNAME 2353 DIMENSION OMMSx(*) , OMMSy(*) , OMMSz(*) 2354C 2355 LOGICAL FILE_EXIST 2356C 2357 INQUIRE(FILE=FLNAME,EXIST=FILE_EXIST) 2358 IF (.NOT. FILE_EXIST) THEN 2359 DO LM = 1,NULOOP 2360 OMMSx(LM) = 0.0D0 2361 OMMSy(LM) = 0.0D0 2362 OMMSz(LM) = 0.0D0 2363 END DO 2364 GO TO 9000 2365 END IF 2366C 2367 2368 NMBU = -1 2369 CALL GPOPEN(NMBU,FLNAME,'OLD',' ','FORMATTED',IDUMMY,.FALSE.) 2370 2371 DO 820 LM = 1,NULOOP 2372 2373 READ(NMBU,'(I5,3E25.15)',END=700,ERR=700) 2374 & LM1, OMMSx(LM), OMMSy(LM), OMMSz(LM) 2375 IF (LM.NE.LM1) THEN 2376 write (lupri,*) 'Error in CC_GET31 on ',FLNAME 2377 write (lupri,*) 'LM, LM1, NULOOP',LM,LM1,NULOOP 2378 CALL QUIT( 'Error in CC_GET31') 2379 END IF 2380 700 CONTINUE 2381 820 CONTINUE 2382 2383 CALL GPCLOSE(NMBU,'KEEP') 2384C 2385 9000 RETURN 2386C 2387 END 2388C 2389C************************************************************** 2390C 2391C /* Deck ccmm_pbtr */ 2392 SUBROUTINE CCMM_PBTR(RHO1,RHO2,ISYRES, 2393 * LISTL,IDLSTL,CTR1,ISYCTR, 2394 * LISTR,IDLSTR,BTR1,ISYBTR, 2395 * MODEL,WORK,LWORK) 2396C 2397C--------------------------------------------------------------- 2398C Calculates contributions from the T^gB , the ^CT^g 2399C and ^CT^gB operators. 2400C--------------------------------------------------------------- 2401C 2402#include "implicit.h" 2403#include "maxorb.h" 2404#include "mxcent.h" 2405 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, IZERO = 0 , TWO = 2.0D0) 2406 PARAMETER (HALF = 0.5D0) 2407#include "dummy.h" 2408#include "iratdef.h" 2409#include "priunit.h" 2410#include "ccorb.h" 2411#include "ccsdsym.h" 2412#include "ccsdio.h" 2413#include "ccsdinp.h" 2414#include "ccfield.h" 2415#include "exeinf.h" 2416#include "ccfdgeo.h" 2417#include "ccslvinf.h" 2418#include "qm3.h" 2419#include "ccinftap.h" 2420C 2421 DIMENSION WORK(LWORK),RHO1(*),RHO2(*),CTR1(*),BTR1(*) 2422C 2423 CHARACTER*(*) LISTL,LISTR,LIST*(2) 2424 CHARACTER*8 LABEL 2425 CHARACTER MODEL*10 2426 LOGICAL LEXIST 2427 INTEGER IDLSTL, IDLSTR, ISYCTR, ISYBTR, ISYBC 2428C 2429 INTEGER KDUM 2430 PARAMETER( KDUM = +99 999 999 ) ! dummy address for work space 2431C 2432 IF (IPRINT.GT.10) THEN 2433 WRITE(LUPRI,*)'CCMM_PBTR: ISYRES =',ISYRES 2434 WRITE(LUPRI,*)'CCMM_PBTR: ISYCTR =',ISYCTR 2435 WRITE(LUPRI,*)'CCMM_PBTR: ISYBTR =',ISYBTR 2436 WRITE(LUPRI,*)'CCMM_PBTR: LISTL =',LISTL 2437 WRITE(LUPRI,*)'CCMM_PBTR: LISTR =',LISTR 2438 WRITE(LUPRI,*)'CCMM_PBTR: IDLSTL =',IDLSTL 2439 WRITE(LUPRI,*)'CCMM_PBTR: IDLSTR =',IDLSTR 2440 CALL FLSHFO(LUPRI) 2441 ENDIF 2442C 2443C---------------------------------------------- 2444C If all models are SPC 2445C -> RETURN from CCMM_TGB: 2446C---------------------------------------------- 2447C 2448 IF (LOSPC) RETURN 2449C 2450C--------------------- 2451C Init parameters. 2452C--------------------- 2453C 2454C For the B (right) trial vector 2455C 2456 NAMP1B = NT1AM(ISYBTR) 2457 NAMP2B = NT2AM(ISYBTR) 2458 NAMPB = NAMP1B + NAMP2B 2459C 2460C For the C (left) trial vector 2461C 2462 NAMP1C = NT1AM(ISYCTR) 2463 NAMP2C = NT2AM(ISYCTR) 2464 NAMPC = NAMP1C + NAMP2C 2465C 2466C For the F = C X B vector 2467C 2468 ISYBC = MULD2H(ISYBTR,ISYCTR) 2469 NAMP1F = NT1AM(ISYBC) 2470 NAMP2F = NT2AM(ISYBC) 2471 NAMPF = NAMP1F + NAMP2F 2472C 2473 IF (CCS) THEN 2474 NT2AM(ISYBTR) = IZERO 2475 NT2AM(ISYCTR) = IZERO 2476 NT2AM(ISYBC) = IZERO 2477 END IF 2478C----------------------------------------------------------- 2479C Calculate contribution from the effective operators. 2480C----------------------------------------------------------- 2481C 2482C First we concentrate on the effective operator T^gB 2483C ----------------------------------------------------- 2484C 2485 KTGB = 1 2486 KWRK1 = KTGB + N2BST(ISYBTR) 2487 LWRK1 = LWORK - KWRK1 2488 IF (LWRK1.LT.0) CALL QUIT( 'Too little work in CCMM_PBTR, 1') 2489C 2490 CALL DZERO(WORK(KTGB),N2BST(ISYBTR)) 2491C 2492 CALL CCMM_TGB(BTR1,ISYBTR,LISTR,IDLSTR,WORK(KTGB),'ET', 2493 * MODEL,WORK(KWRK1),LWRK1) 2494C 2495C Symmetry: 2496C --------- 2497C 2498 ISYBC = MULD2H(ISYBTR,ISYCTR) 2499C 2500 IF (ISYBC .NE. ISYRES) THEN 2501 CALL QUIT( 'Symmetry problem in CCMM_PBTR') 2502 END IF 2503C 2504 KETA = KWRK1 2505 KWRK2 = KETA + NAMPF 2506 LWRK2 = LWORK - KWRK2 2507 IF (LWRK2.LT.0) CALL QUIT( 'Too little work in CCMM_PBTR, 2') 2508C 2509C Note, LISTL .EQ. LE/L1 so the HF part of the following 2510C eta-transformation is skipped 2511C 2512 LABEL = 'GIVE INT' 2513 CALL CC_ETAC(ISYBTR,LABEL,WORK(KETA), 2514 * LISTL,IDLSTL,0,WORK(KTGB),WORK(KWRK2),LWRK2) 2515C 2516 KETA1 = KETA 2517 KETA2 = KETA1 + NAMP1F 2518C 2519 CALL DAXPY(NAMP1F,ONE,WORK(KETA1),1,RHO1,1) 2520 CALL DAXPY(NAMP2F,ONE,WORK(KETA2),1,RHO2,1) 2521C 2522C Next, the effective operator ^CT^gB is considered 2523C---------------------------------------------------------- 2524C 2525C Symmetry of the operator: 2526C ------------------------- 2527C 2528 ISYBC = MULD2H(ISYBTR,ISYCTR) 2529C 2530 KTGB = 1 2531 KWRK1 = KTGB + N2BST(ISYBC) 2532 LWRK1 = LWORK - KWRK1 2533 IF (LWRK1.LT.0) CALL QUIT( 'Too little work in CCMM_PBTR, 3') 2534C 2535 CALL DZERO(WORK(KTGB),N2BST(ISYBC)) 2536C 2537 CALL CCMM_CTGB(LISTL,IDLSTL,CTR1,ISYCTR, 2538 * LISTR,IDLSTR,BTR1,ISYBTR, 2539 * WORK(KTGB),MODEL,WORK(KWRK1),LWRK1) 2540C 2541 KETA = KWRK1 2542 KWRK2 = KETA + NAMPF 2543 LWRK2 = LWORK - KWRK2 2544 IF (LWRK2.LT.0) CALL QUIT( 'Too little work in CCMM_PBTR, 4') 2545C 2546 CALL DZERO(WORK(KETA),NAMPF) 2547C 2548 LABEL = 'GIVE INT' 2549 LIST = 'L0' 2550 IDLINO = 1 2551C 2552 CALL CC_ETAC(ISYBC,LABEL,WORK(KETA), 2553 * LIST,IDLINO,0,WORK(KTGB),WORK(KWRK2),LWRK2) 2554C 2555 KETA1 = KETA 2556 KETA2 = KETA1 + NAMP1F 2557C 2558 CALL DAXPY(NT1AM(ISYBC),ONE,WORK(KETA1),1,RHO1,1) 2559 CALL DAXPY(NT2AM(ISYBC),ONE,WORK(KETA2),1,RHO2,1) 2560C 2561C Finally, we calculate the third contribution which is 2562C a contribution from a T^gB operator but calculated 2563C as a xi vector element. 2564C----------------------------------------------------------- 2565C 2566 KTGB = 1 2567 KWRK1 = KTGB + N2BST(ISYCTR) 2568 LWRK1 = LWORK - KWRK1 2569 IF (LWRK1.LT.0) CALL QUIT( 'Too little work in CCMM_PBTR, 5') 2570C 2571 CALL DZERO(WORK(KTGB),N2BST(ISYCTR)) 2572C 2573 CALL CCMM_TGB(CTR1,ISYCTR,LISTL,IDLSTL,WORK(KTGB),'XI', 2574 * MODEL,WORK(KWRK1),LWRK1) 2575C 2576C Symmetry: 2577C --------- 2578C 2579 ISYBC = MULD2H(ISYBTR,ISYCTR) 2580C 2581 IF (ISYBC .NE. ISYRES) THEN 2582 CALL QUIT( 'Symmetry problem in CCMM_PBTR') 2583 END IF 2584 2585 LABEL = 'GIVE INT' 2586 LIST = 'L0' 2587 LISTNO = 1 2588C 2589C (NB: Result vector from the following transformation is 2590C given at the beginning of WORK(KWRK1).) 2591C 2592 CALL CCLR_FA(LABEL,ISYCTR,LISTR,IDLSTR, 2593 & LIST,LISTNO,WORK(KTGB),WORK(KWRK1),LWRK1) 2594C 2595 CALL DAXPY(NT1AM(ISYBC),ONE,WORK(KWRK1),1,RHO1,1) 2596 CALL DAXPY(NT2AM(ISYBC),ONE, 2597 * WORK(KWRK1+NT1AM(ISYBC)),1,RHO2,1) 2598C 2599 END 2600******************************************************************* 2601C /* Deck ccmm_tgb */ 2602 SUBROUTINE CCMM_TGB(CTR1,ISYMTR,LISTIN,IDLIST,TGB,LR, 2603 * MODEL,WORK,LWORK) 2604C 2605C----------------------------------------------------------------------------- 2606C 2607C IF (LR.EQ.'ET') then the following effective operator is constructed 2608C T^gB = - Sum_a * Sum_sigma 2609C <Lambda|[Rr_a,Tau_sigma]|CC> t^B_sigma Rr_a 2610C 2611C IF (LR.EQ.'XI') then the following effective operator is constructed 2612C ^BT^g = - Sum_a * Sum_sigma 2613C t^B_sigma <bar sigma|Rr_a|CC> Rr_a 2614C 2615C JK+OC, marts 03 2616C KS, modifed to allow for J term needed from the left. I.e when cG operator is needed 2617C in the NYQMMM formulation. This routine is not used for calculating Gc term (LR=ET) - see 2618C instead CCMM_TRANSFORMER in cc_qmmm 2619C----------------------------------------------------------------------------- 2620C 2621 USE PELIB_INTERFACE, ONLY: USE_PELIB, PELIB_IFC_RESPONSE 2622#include "implicit.h" 2623#include "maxorb.h" 2624#include "mxcent.h" 2625 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, IZERO = 0 , TWO = 2.0D0) 2626 PARAMETER (HALF = 0.5D0) 2627#include "dummy.h" 2628#include "iratdef.h" 2629#include "priunit.h" 2630#include "ccorb.h" 2631#include "ccsdsym.h" 2632#include "ccsdio.h" 2633#include "ccsdinp.h" 2634#include "ccfield.h" 2635#include "exeinf.h" 2636#include "ccfdgeo.h" 2637#include "ccslvinf.h" 2638#include "qm3.h" 2639#include "ccinftap.h" 2640C 2641 DIMENSION WORK(LWORK),CTR1(*) 2642 DIMENSION TGB(*) 2643C 2644 INTEGER ISYMTR, IDLIST 2645 INTEGER KDUM 2646 PARAMETER( KDUM = +99 999 999 ) ! dummy address for work space 2647C 2648 CHARACTER*(*) LISTIN 2649 CHARACTER MODEL*10 2650C 2651 CHARACTER*8 LABEL,LR*(2),LIST 2652 CHARACTER*(8) FILMME, FILMMX 2653 2654 LOGICAL LEXIST, LOCDEB 2655 PARAMETER (LOCDEB=.FALSE.) 2656 2657 REAL*8, ALLOCATABLE :: DENMATS(:), JTERMTEMP(:) 2658C 2659C 2660 IF (IPRINT.GT.10) THEN 2661 WRITE(LUPRI,*)'CCMM_TGB : ISYMTR:', ISYMTR 2662 ENDIF 2663C 2664C---------------------------------------------- 2665C If all models are SPC 2666C -> RETURN from CCMM_TGB: 2667C---------------------------------------------- 2668C 2669 IF (LOSPC) RETURN 2670C 2671C--------------------- 2672C Init parameters. 2673C--------------------- 2674C 2675 NTAMP1 = NT1AM(ISYMTR) 2676 NTAMP2 = NT2AM(ISYMTR) 2677 NTAMP = NTAMP1 + NTAMP2 2678C 2679 IF (CCS) NT2AM(ISYMTR) = IZERO 2680C 2681 IF (DISCEX) THEN 2682 LUMMET = -1 2683 LUMMXI = -1 2684 FILMME = 'CCMM_ETA' 2685 FILMMX = 'CCMM__XI' 2686 CALL WOPEN2(LUMMET, FILMME, 64, 0) 2687 CALL WOPEN2(LUMMXI, FILMMX, 64, 0) 2688 END IF 2689C 2690C---------------------------------------- 2691C Readin trial vector from file. 2692C---------------------------------------- 2693C 2694 KT2AMPA = 1 2695 KWRK1 = KT2AMPA + NT2AM(ISYMTR) 2696 LWRK1 = LWORK - KWRK1 2697C 2698 IF (LWRK1 .LT. 0) THEN 2699 CALL QUIT('Insuff. work in CCMM_TGB 1') 2700 END IF 2701C 2702 IOPT = 2 2703 CALL CC_RDRSP(LISTIN,IDLIST,ISYMTR,IOPT,MODEL, 2704 * WORK(KDUM),WORK(KT2AMPA)) 2705 2706 IF (LOCDEB) THEN 2707 RHO2N = DDOT(NT2AM(ISYMTR),WORK(KT2AMPA),1,WORK(KT2AMPA),1) 2708 WRITE(LUPRI,*)'Norm af B2AM in CCMM_TGB on input:,1', RHO2N 2709 END IF 2710C 2711 IF (.NOT. (NYQMMM .OR. USE_PELIB())) THEN 2712C------------------------------------ 2713C Dynamical allocation for CCMM : 2714C------------------------------------ 2715C 2716 KRAAOx = KWRK1 2717 KRAAOy = KRAAOx + N2BST(ISYMTR) 2718 KRAAOz = KRAAOy + N2BST(ISYMTR) 2719 KWRK2 = KRAAOz + N2BST(ISYMTR) 2720 LWRK2 = LWORK - KWRK2 2721C 2722 IF (LWRK2.LT.0) CALL QUIT( 'Too little work in CCMM_TGB, 2') 2723C 2724 KXIx = KWRK2 2725 KXIy = KXIx + NTAMP 2726 KXIz = KXIy + NTAMP 2727 KWRK3 = KXIz + NTAMP 2728 LWRK3 = LWORK - KWRK3 2729C 2730 IF (LWRK3.LT.0) CALL QUIT( 'Too little work in CCMM_TGB, 3') 2731C 2732 CALL DZERO(WORK(KRAAOx),3*N2BST(ISYMTR)) 2733 CALL DZERO(WORK(KXIx),3*NTAMP) 2734C 2735C ---------------------------- 2736C Read Rra in AO basis: 2737C ---------------------------- 2738C 2739 LUCCEF = -1 2740 CALL GPOPEN(LUCCEF,'ELFDMM','OLD',' ', 2741 & 'UNFORMATTED',IDUMMY,.FALSE.) 2742 REWIND(LUCCEF) 2743C 2744 LM = 0 2745 IADR1 = 1 2746 IADR2 = 1 2747 LEN = NTAMP 2748C 2749 DO 700 I = 1, ISYTP 2750 IF (MDLWRD(I)(1:5) .EQ. 'SPC_E') THEN 2751 DO 701 J = NSYSBG(I), NSYSED(I) 2752 DO 702 K = 1, NUALIS(I) 2753C 2754 LM = LM + 1 2755 LABEL = 'READ INT' 2756C 2757 CALL CCMMINT(LABEL,LUCCEF,WORK(KRAAOx),DUMMY,IRREP, 2758 & ISYM,IERR,WORK(KWRK3),LWRK3) 2759 CALL CCMMINT(LABEL,LUCCEF,WORK(KRAAOy),DUMMY,IRREP, 2760 & ISYM,IERR,WORK(KWRK3),LWRK3) 2761 CALL CCMMINT(LABEL,LUCCEF,WORK(KRAAOz),DUMMY,IRREP, 2762 & ISYM,IERR,WORK(KWRK3),LWRK3) 2763C 2764C--------------------------------------------------------------------- 2765C (LR.EQ.'ET') Calculate eta^T_lm 2766C (LR.EQ.'XI') Calculate xi^T_lm 2767C Only if MDLWRD(I) = 'SPC_E01' are we going to calculate 2768C contributions from the T^g operator to the response 2769C equations. 2770C--------------------------------------------------------------------- 2771C 2772 IF (MDLWRD(I) .EQ. 'SPC_E01') THEN 2773 IF (.NOT. (DISCEX)) THEN 2774 IF (LR.EQ.'ET') THEN 2775 LABEL = 'GIVE INT' 2776 LIST = 'L0' 2777 IDLINO = 1 2778C 2779 CALL CC_ETAC(ISYMTR,LABEL,WORK(KXIx), 2780 * LIST,IDLINO,0,WORK(KRAAOx),WORK(KWRK3),LWRK3) 2781 CALL CC_ETAC(ISYMTR,LABEL,WORK(KXIy), 2782 * LIST,IDLINO,0,WORK(KRAAOy),WORK(KWRK3),LWRK3) 2783 CALL CC_ETAC(ISYMTR,LABEL,WORK(KXIz), 2784 * LIST,IDLINO,0,WORK(KRAAOz),WORK(KWRK3),LWRK3) 2785C 2786 ELSE IF (LR.EQ.'XI') THEN 2787 LABEL = 'GIVE INT' 2788C 2789 CALL CC_XKSI(WORK(KXIx),LABEL,ISYMTR,0,WORK(KRAAOx), 2790 * WORK(KWRK3),LWRK3) 2791 CALL CC_XKSI(WORK(KXIy),LABEL,ISYMTR,0,WORK(KRAAOy), 2792 * WORK(KWRK3),LWRK3) 2793 CALL CC_XKSI(WORK(KXIz),LABEL,ISYMTR,0,WORK(KRAAOz), 2794 * WORK(KWRK3),LWRK3) 2795C 2796 END IF 2797 ELSE 2798 IF (LR.EQ.'ET') THEN 2799 CALL GETWA2(LUMMET,FILMME,WORK(KXIx),IADR1,LEN) 2800 IADR1 = IADR1 + LEN 2801 CALL GETWA2(LUMMET,FILMME,WORK(KXIy),IADR1,LEN) 2802 IADR1 = IADR1 + LEN 2803 CALL GETWA2(LUMMET,FILMME,WORK(KXIz),IADR1,LEN) 2804 IADR1 = IADR1 + LEN 2805 ELSE IF (LR.EQ.'XI') THEN 2806 CALL GETWA2(LUMMXI,FILMMX,WORK(KXIx),IADR2,LEN) 2807 IADR2 = IADR2 + LEN 2808 CALL GETWA2(LUMMXI,FILMMX,WORK(KXIy),IADR2,LEN) 2809 IADR2 = IADR2 + LEN 2810 CALL GETWA2(LUMMXI,FILMMX,WORK(KXIz),IADR2,LEN) 2811 IADR2 = IADR2 + LEN 2812 END IF 2813 END IF 2814C 2815C-------------------------------------------------------------- 2816C Contract with trial vector: 2817C (If LR.EQ.'ET' then from the right hand side) 2818C (If LR.EQ.'XI' then from the left hand side) 2819C-------------------------------------------------------------- 2820C 2821 KXI1x = KXIx 2822 KXI2x = KXIx + NTAMP1 2823 BXILMD1x= DDOT(NTAMP1,CTR1,1,WORK(KXI1x),1) 2824 BXILMD2x= DDOT(NTAMP2,WORK(KT2AMPA),1,WORK(KXI2x),1) 2825 BXILMDx = BXILMD1x + BXILMD2x 2826C 2827 KXI1y = KXIy 2828 KXI2y = KXIy + NTAMP1 2829 BXILMD1y = DDOT(NTAMP1,CTR1,1,WORK(KXI1y),1) 2830 BXILMD2y = DDOT(NTAMP2,WORK(KT2AMPA),1,WORK(KXI2y),1) 2831 BXILMDy = BXILMD1y + BXILMD2y 2832C 2833 KXI1z = KXIz 2834 KXI2z = KXIz + NTAMP1 2835 BXILMD1z = DDOT(NTAMP1,CTR1,1,WORK(KXI1z),1) 2836 BXILMD2z = DDOT(NTAMP2,WORK(KT2AMPA),1,WORK(KXI2z),1) 2837 BXILMDz = BXILMD1z + BXILMD2z 2838C 2839 FACx = - ALPIMM(I,K) * (BXILMDx) 2840 FACy = - ALPIMM(I,K) * (BXILMDy) 2841 FACz = - ALPIMM(I,K) * (BXILMDz) 2842C 2843 IF (IQM3PR .GT. 5) THEN 2844 WRITE(LUPRI,*)'---------------------------' // 2845 * '-------------------' 2846 WRITE(LUPRI,*)'The factors calculated at cent.' // 2847 * ' of mass no.',LM 2848 WRITE(LUPRI,*)'---------------------------' // 2849 * '-------------------' 2850 WRITE(LUPRI,*)'FACx=',FACx 2851 WRITE(LUPRI,*)'FACy=',FACy 2852 WRITE(LUPRI,*)'FACz=',FACz 2853 WRITE(LUPRI,*)'---------------------------' // 2854 * '-------------------' 2855 END IF 2856C 2857 2858C -------------------------------------------------------------- 2859C Put factor on integrals 2860C -------------------------------------------------------------- 2861C 2862 CALL DAXPY(N2BST(ISYMTR),FACx,WORK(KRAAOx),1, 2863 * TGB,1) 2864 CALL DAXPY(N2BST(ISYMTR),FACy,WORK(KRAAOy),1, 2865 * TGB,1) 2866 CALL DAXPY(N2BST(ISYMTR),FACz,WORK(KRAAOz),1, 2867 * TGB,1) 2868C 2869 END IF 2870 702 CONTINUE 2871 701 CONTINUE 2872 END IF 2873 700 CONTINUE 2874C 2875 CALL GPCLOSE(LUCCEF,'KEEP') 2876C 2877 ELSE IF (NYQMMM) THEN 2878 2879 IF (LR .NE. 'XI' ) CALL QUIT('CCMM_TGB xi flag not allowed') 2880 KDENS = KWRK1 2881 KWRK2 = KDENS + N2BST(ISYMTR) 2882 LWRK2 = LWORK - KWRK2 2883 2884 CALL DZERO(WORK(KDENS),N2BST(ISYMTR)) 2885 2886C---------------------------- 2887C Construct auxiliary densities 2888C---------------------------- 2889 CALL CCMM_D1AO(WORK(KDENS),CTR1,WORK(KT2AMPA),MODEL,'L', 2890 * LISDUM,IDLDUM,ISYDUM, 2891 * WORK(KWRK2),LWRK2) 2892C---------------------------- 2893C Construct effective G operator 2894C---------------------------- 2895 2896 CALL CCMM_GEFF(WORK(KDENS),TGB(1),WORK(KWRK2),LWRK2) 2897 2898 ELSE IF (USE_PELIB()) THEN 2899 2900 IF (LR .NE. 'XI' ) CALL QUIT('CCMM_TGB XI flag now allowed!') 2901 KDENS = KWRK1 2902 KWRK2 = KDENS + N2BST(ISYMTR) 2903 LWRK2 = LWORK - KWRK2 2904 2905 CALL DZERO(WORK(KDENS),N2BST(ISYMTR)) 2906 2907 ALLOCATE(DENMATS(NNBASX),JTERMTEMP(NNBASX)) 2908 2909 JTERMTEMP = 0.0d0 2910 CALL CCMM_D1AO(WORK(KDENS),CTR1,WORK(KT2AMPA),MODEL,'L', 2911 & LISDUM,IDLDUM,ISYDUM,WORK(KWRK2),LWRK2) 2912 CALL DGEFSP(NBAS,WORK(KDENS),DENMATS) 2913 CALL PELIB_IFC_RESPONSE(DENMATS,JTERMTEMP,1) 2914 CALL DSPTSI(NBAS,JTERMTEMP,TGB(1)) 2915 2916 DEALLOCATE(DENMATS,JTERMTEMP) 2917 END IF 2918 2919 IF (DISCEX) THEN 2920 CALL WCLOSE2(LUMMET, FILMME, 'KEEP') 2921 CALL WCLOSE2(LUMMXI, FILMMX, 'KEEP') 2922 END IF 2923C 2924 END 2925************************************************************ 2926C /* Deck ccmm_ctgb */ 2927 SUBROUTINE CCMM_CTGB(LISTL,IDLSTL,CTR1,ISYCTR, 2928 * LISTR,IDLSTR,BTR1,ISYBTR, 2929 * TGB,MODEL,WORK,LWORK) 2930C 2931C----------------------------------------------------------------------------- 2932C 2933C Construcs effective operator equal to 2934C ^CT^gB = - Sum_a * Sum_sigma Sum_my 2935C t^C_my <bar my|[Rr_a,Tau_sigma]|CC> t^B_sigma * Rr_a 2936C 2937C JK+OC, marts 03 2938C----------------------------------------------------------------------------- 2939C 2940#include "implicit.h" 2941#include "maxorb.h" 2942#include "mxcent.h" 2943 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, IZERO = 0 , TWO = 2.0D0) 2944 PARAMETER (HALF = 0.5D0) 2945#include "dummy.h" 2946#include "iratdef.h" 2947#include "priunit.h" 2948#include "ccorb.h" 2949#include "ccsdsym.h" 2950#include "ccsdio.h" 2951#include "ccsdinp.h" 2952#include "ccfield.h" 2953#include "exeinf.h" 2954#include "ccfdgeo.h" 2955#include "ccslvinf.h" 2956#include "qm3.h" 2957#include "ccinftap.h" 2958C 2959 DIMENSION WORK(LWORK),BTR1(*),CTR1(*) 2960 DIMENSION TGB(*) 2961C 2962 CHARACTER*8 LABEL, LIST*(2) 2963 CHARACTER*(*) LISTL, LISTR 2964 CHARACTER MODEL*10 2965 INTEGER IDLSTL, IDLSTR, ISYCTR, ISYBTR 2966 LOGICAL LEXIST 2967C 2968 IF (IPRINT.GT.10) THEN 2969 WRITE(LUPRI,*)'CCMM_CTGB: ISYCTR =',ISYCTR 2970 WRITE(LUPRI,*)'CCMM_CTGB: ISYBTR =',ISYBTR 2971 WRITE(LUPRI,*)'CCMM_CTGB: LISTL =',LISTL 2972 WRITE(LUPRI,*)'CCMM_CTGB: LISTR =',LISTR 2973 WRITE(LUPRI,*)'CCMM_CTGB: IDLSTL =',IDLSTL 2974 WRITE(LUPRI,*)'CCMM_CTGB: IDLSTR =',IDLSTR 2975 CALL FLSHFO(LUPRI) 2976 ENDIF 2977C 2978C---------------------------------------------- 2979C If all models are SPC 2980C -> RETURN from CCMM_TGB: 2981C---------------------------------------------- 2982C 2983 IF (LOSPC) RETURN 2984C 2985C--------------------- 2986C Init parameters. 2987C--------------------- 2988C 2989C For the B (right) trial vector 2990C 2991 NAMP1B = NT1AM(ISYBTR) 2992 NAMP2B = NT2AM(ISYBTR) 2993 NAMPB = NAMP1B + NAMP2B 2994C 2995C For the C (left) trial vector 2996C 2997 NAMP1C = NT1AM(ISYCTR) 2998 NAMP2C = NT2AM(ISYCTR) 2999 NAMPC = NAMP1C + NAMP2C 3000C 3001C For the F = B x C vector 3002C 3003 ISYBC = MULD2H(ISYBTR,ISYCTR) 3004 NAMP1F = NT1AM(ISYBC) 3005 NAMP2F = NT2AM(ISYBC) 3006 NAMPF = NAMP1F + NAMP2F 3007C 3008 IF (CCS) THEN 3009 NT2AM(ISYBTR) = IZERO 3010 NT2AM(ISYCTR) = IZERO 3011 NT2AM(ISYBC) = IZERO 3012 END IF 3013 3014C----------------------------------------------- 3015C Readin response vector (B, right) from file 3016C----------------------------------------------- 3017C 3018 KT2AMPA = 1 3019 KWRK1 = KT2AMPA + NT2AM(ISYBTR) 3020 LWRK1 = LWORK - KWRK1 3021C 3022 IF (LWRK1 .LT. 0) THEN 3023 CALL QUIT('Insuff. work in CCMM_TGB') 3024 END IF 3025C 3026 IOPT = 2 3027 CALL CC_RDRSP(LISTR,IDLSTR,ISYBTR,IOPT,MODEL, 3028 * WORK(KWRK1),WORK(KT2AMPA)) 3029C 3030C------------------------------------ 3031C Dynamical allocation for CCMM : 3032C------------------------------------ 3033C 3034 KRAAOx = KWRK1 3035 KRAAOy = KRAAOx + N2BST(ISYBC) 3036 KRAAOz = KRAAOy + N2BST(ISYBC) 3037 KWRK2 = KRAAOz + N2BST(ISYBC) 3038 LWRK2 = LWORK - KWRK2 3039C 3040 IF (LWRK2.LT.0) CALL QUIT( 'Too little work in CCMM_CTGB, 2') 3041C 3042 KXIx = KWRK2 3043 KXIy = KXIx + NAMPB 3044 KXIz = KXIy + NAMPB 3045 KWRK3 = KXIz + NAMPB 3046 LWRK3 = LWORK - KWRK3 3047C 3048 IF (LWRK3.LT.0) CALL QUIT( 'Too little work in CCMM_CTGB, 3') 3049C 3050 CALL DZERO(WORK(KRAAOx),3*N2BST(ISYBC)) 3051 CALL DZERO(WORK(KXIx),3*NAMPB) 3052C 3053C---------------------------- 3054C Read Rra in AO basis: 3055C---------------------------- 3056C 3057 LUCCPO = -1 3058 CALL GPOPEN(LUCCEF,'ELFDMM','OLD',' ', 3059 & 'UNFORMATTED',IDUMMY,.FALSE.) 3060 REWIND(LUCCEF) 3061C 3062 LM = 0 3063C 3064 DO 700 I = 1, ISYTP 3065 IF (MDLWRD(I)(1:5) .EQ. 'SPC_E') THEN 3066 DO 701 J = NSYSBG(I), NSYSED(I) 3067 DO 702 K = 1, NUALIS(I) 3068C 3069 LM = LM + 1 3070 LABEL = 'READ INT' 3071C 3072 CALL CCMMINT(LABEL,LUCCEF,WORK(KRAAOx),DUMMY,IRREP, 3073 & ISYM,IERR,WORK(KWRK3),LWRK3) 3074 CALL CCMMINT(LABEL,LUCCEF,WORK(KRAAOy),DUMMY,IRREP, 3075 & ISYM,IERR,WORK(KWRK3),LWRK3) 3076 CALL CCMMINT(LABEL,LUCCEF,WORK(KRAAOz),DUMMY,IRREP, 3077 & ISYM,IERR,WORK(KWRK3),LWRK3) 3078C 3079C--------------------------------------------------------------------- 3080C Calculate eta^Rr_a 3081C (LIST = 'LE' ensures that the HF part of eta^Rr_a is skipped. 3082C--------------------------------------------------------------------- 3083C 3084 IF (MDLWRD(I) .EQ. 'SPC_E01') THEN 3085 LABEL = 'GIVE INT' 3086 CALL CC_ETAC(ISYBC,LABEL,WORK(KXIx), 3087 * LISTL,IDLSTL,0,WORK(KRAAOx),WORK(KWRK3),LWRK3) 3088 CALL CC_ETAC(ISYBC,LABEL,WORK(KXIy), 3089 * LISTL,IDLSTL,0,WORK(KRAAOy),WORK(KWRK3),LWRK3) 3090 CALL CC_ETAC(ISYBC,LABEL,WORK(KXIz), 3091 * LISTL,IDLSTL,0,WORK(KRAAOz),WORK(KWRK3),LWRK3) 3092C 3093C--------------------------------------------------------------- 3094C Contract with trial vector B (from the right hand side) 3095C--------------------------------------------------------------- 3096C 3097 KXI1x = KXIx 3098 KXI2x = KXIx + NT1AM(ISYBTR) 3099 BXILMD1x= DDOT(NT1AM(ISYBTR),BTR1,1,WORK(KXI1x),1) 3100 BXILMD2x= DDOT(NT2AM(ISYBTR),WORK(KT2AMPA), 3101 * 1,WORK(KXI2x),1) 3102 BXILMDx = BXILMD1x + BXILMD2x 3103C 3104 KXI1y = KXIy 3105 KXI2y = KXIy + NT1AM(ISYBTR) 3106 BXILMD1y= DDOT(NT1AM(ISYBTR),BTR1,1,WORK(KXI1y),1) 3107 BXILMD2y= DDOT(NT2AM(ISYBTR),WORK(KT2AMPA), 3108 * 1,WORK(KXI2y),1) 3109 BXILMDy = BXILMD1y + BXILMD2y 3110C 3111 KXI1z = KXIz 3112 KXI2z = KXIz + NT1AM(ISYBTR) 3113 BXILMD1z= DDOT(NT1AM(ISYBTR),BTR1,1,WORK(KXI1z),1) 3114 BXILMD2z= DDOT(NT2AM(ISYBTR),WORK(KT2AMPA), 3115 * 1,WORK(KXI2z),1) 3116 BXILMDz = BXILMD1z + BXILMD2z 3117C 3118 FACx = - ALPIMM(I,K) * (BXILMDx) 3119 FACy = - ALPIMM(I,K) * (BXILMDy) 3120 FACz = - ALPIMM(I,K) * (BXILMDz) 3121C 3122 IF (IQM3PR .GT. 5) THEN 3123 WRITE(LUPRI,*)'---------------------------' // 3124 * '-------------------' 3125 WRITE(LUPRI,*)'The factors calculated at cent.' // 3126 * ' of mass no.',LM 3127 WRITE(LUPRI,*)'---------------------------' // 3128 * '-------------------' 3129 WRITE(LUPRI,*)'FACx=',FACx 3130 WRITE(LUPRI,*)'FACy=',FACy 3131 WRITE(LUPRI,*)'FACz=',FACz 3132 WRITE(LUPRI,*)'---------------------------' // 3133 * '-------------------' 3134 END IF 3135C 3136C-------------------------------------------------------------- 3137C Put factor on integrals 3138C-------------------------------------------------------------- 3139C 3140 CALL DAXPY(N2BST(ISYBC),FACx,WORK(KRAAOx),1, 3141 * TGB,1) 3142 CALL DAXPY(N2BST(ISYBC),FACy,WORK(KRAAOy),1, 3143 * TGB,1) 3144 CALL DAXPY(N2BST(ISYBC),FACz,WORK(KRAAOz),1, 3145 * TGB,1) 3146C 3147C--------------------------------------------------------------- 3148C 3149 END IF 3150 702 CONTINUE 3151 701 CONTINUE 3152 END IF 3153 700 CONTINUE 3154C 3155 CALL GPCLOSE(LUCCEF,'KEEP') 3156 END 3157C******************************************************************* 3158C /* Deck ccmm_btr */ 3159 SUBROUTINE CCMM_BTR(RHO1,RHO2,ISYMAB, 3160 * LISTA,IDLSTA,ISYMA, 3161 * LISTB,IDLSTB,ISYMB, 3162 * MODEL,WORK,LWORK) 3163C 3164C Calculates the CC/MM contribution to the B-matrix 3165C transformation. 3166C JK+OC, marts 03 3167C--------------------------------------------------------------- 3168C 3169#include "implicit.h" 3170#include "maxorb.h" 3171#include "mxcent.h" 3172 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, IZERO = 0 , TWO = 2.0D0) 3173 PARAMETER (HALF = 0.5D0) 3174#include "dummy.h" 3175#include "iratdef.h" 3176#include "priunit.h" 3177#include "ccorb.h" 3178#include "ccsdsym.h" 3179#include "ccsdio.h" 3180#include "ccsdinp.h" 3181#include "ccfield.h" 3182#include "exeinf.h" 3183#include "ccfdgeo.h" 3184#include "ccslvinf.h" 3185#include "qm3.h" 3186#include "ccinftap.h" 3187C 3188 DIMENSION WORK(LWORK),RHO1(*),RHO2(*) 3189C 3190 CHARACTER*(*) LISTA,LISTB 3191 CHARACTER*(3) LISTC 3192 CHARACTER*8 LABEL 3193 CHARACTER MODEL*10 3194 LOGICAL LEXIST, LSAME, LOCDEB 3195 PARAMETER (LOCDEB=.FALSE.) 3196 INTEGER IDLSTA, IDLSTB, ISYMAB, ISYMA, ISYMB 3197 PARAMETER( KDUM = +99 999 999 ) ! dummy address for work space 3198C 3199C---------------------------------------------- 3200C If all models are SPC 3201C -> RETURN from CCMM_TGB: 3202C---------------------------------------------- 3203C 3204 IF (LOSPC) RETURN 3205C 3206C 3207 IF (IPRINT.GT.10) THEN 3208 WRITE(LUPRI,*)'CCMM_BTR: ISYMAB =',ISYMAB 3209 WRITE(LUPRI,*)'CCMM_BTR: ISYMA =',ISYMA 3210 WRITE(LUPRI,*)'CCMM_BTR: ISYMB =',ISYMB 3211 WRITE(LUPRI,*)'CCMM_BTR: LISTA =',LISTA 3212 WRITE(LUPRI,*)'CCMM_BTR: LISTB =',LISTB 3213 WRITE(LUPRI,*)'CCMM_BTR: IDLSTA =',IDLSTA 3214 WRITE(LUPRI,*)'CCMM_BTR: IDLSTB =',IDLSTB 3215 CALL FLSHFO(LUPRI) 3216 ENDIF 3217C 3218 IF (CCS) THEN 3219 NT2AM(ISYMAB) = IZERO 3220 NT2AM(ISYMA) = IZERO 3221 NT2AM(ISYMB) = IZERO 3222 END IF 3223C 3224C First we concentrate on the effective operator T^gA 3225C ----------------------------------------------------- 3226C 3227 KT1AMPA = 1 3228 KTGA = KT1AMPA + NT1AM(ISYMA) 3229 KWRK1 = KTGA + N2BST(ISYMA) 3230 LWRK1 = LWORK - KWRK1 3231C 3232 IF (LWRK1 .LT. 0) THEN 3233 CALL QUIT('Insuff. work in CCMM_BTR 1') 3234 END IF 3235C 3236 CALL DZERO(WORK(KT1AMPA),NT1AM(ISYMA)) 3237 CALL DZERO(WORK(KTGA),N2BST(ISYMA)) 3238C 3239C Read one electron excitation part of response vector 3240C 3241 IOPT = 1 3242 CALL CC_RDRSP(LISTA,IDLSTA,ISYMA,IOPT,MODEL, 3243 * WORK(KT1AMPA),WORK(KDUM)) 3244 3245 IF (LOCDEB) THEN 3246 RHO1N = DDOT(NT1AM(ISYMAB),RHO1,1,RHO1,1) 3247 RHO2N = DDOT(NT2AM(ISYMAB),RHO2,1,RHO2,1) 3248 WRITE(LUPRI,*)'Norm of RHO1 in CCMM_QRTRANS on input:,1', RHO1N 3249 WRITE(LUPRI,*)'Norm of RHO2 in CCMM_QRTRANS on input:,1', RHO2N 3250 RHO1N = DDOT(NT1AM(ISYMA),WORK(KT1AMPA),1,WORK(KT1AMPA),1) 3251 WRITE(LUPRI,*)'Norm af B1AM in CCMM_QRTRANS on input:,1', RHO1N 3252 END IF 3253C 3254 CALL CCMM_TGB(WORK(KT1AMPA),ISYMA,LISTA,IDLSTA,WORK(KTGA),'ET', 3255 * MODEL,WORK(KWRK1),LWRK1) 3256C 3257 LABEL = 'GIVE INT' 3258 CALL CCCR_AA(LABEL,ISYMA,LISTB,IDLSTB, 3259 & WORK(KTGA),WORK(KWRK1),LWRK1) 3260 CALL CCLR_DIASCL(WORK(KWRK1+NT1AM(ISYMAB)),2.0d0,ISYMAB) 3261C 3262C 3263 LSAME = (LISTA.EQ.LISTB .AND. IDLSTA.EQ.IDLSTB) 3264 IF (LSAME) THEN 3265 FACTSLV = TWO 3266 ELSE 3267 FACTSLV = ONE 3268 END IF 3269C 3270 CALL DAXPY(NT1AM(ISYMAB),FACTSLV,WORK(KWRK1),1,RHO1,1) 3271 CALL DAXPY(NT2AM(ISYMAB),FACTSLV, 3272 * WORK(KWRK1+NT1AM(ISYMAB)),1,RHO2,1) 3273 3274 IF (LOCDEB) THEN 3275 TAL1= DDOT(NT1AM(ISYMAB),WORK(KWRK1),1,WORK(KWRK1),1) 3276 TAL2= DDOT(NT2AM(ISYMAB),WORK(KWRK1+NT1AM(ISYMAB)),1, 3277 * WORK(KWRK1+NT1AM(ISYMAB)),1) 3278 WRITE(LUPRI,*) 'Printing first contribution. 3279 & Norm2 of singles: ', TAL1, 3280 & 'Norm2 of doubles: ', TAL2 3281 TAL1= DDOT(NT1AM(ISYMAB),RHO1,1,RHO1,1) 3282 TAL2= DDOT(NT2AM(ISYMAB),RHO2,1,RHO2,1) 3283 WRITE(LUPRI,*) 'Printing RHO total. 3284 & Norm2 of singles: ', TAL1, 3285 & 'Norm2 of doubles: ', TAL2 3286 END IF 3287C 3288C 3289 IF (.NOT.(LSAME)) THEN 3290C 3291C Next we concentrate on the effective operator T^gB 3292C ----------------------------------------------------- 3293C 3294 KT1AMPB = 1 3295 KTGB = KT1AMPB + NT1AM(ISYMB) 3296 KWRK1 = KTGB + N2BST(ISYMB) 3297 LWRK1 = LWORK - KWRK1 3298C 3299 IF (LWRK1 .LT. 0) THEN 3300 CALL QUIT('Insuff. work in CCMM_BTR 2') 3301 END IF 3302C 3303 CALL DZERO(WORK(KT1AMPB),NT1AM(ISYMB)) 3304 CALL DZERO(WORK(KTGB),N2BST(ISYMB)) 3305C 3306C Read one electron excitation part of response vector 3307C 3308 IOPT = 1 3309 CALL CC_RDRSP(LISTB,IDLSTB,ISYMB,IOPT,MODEL, 3310 * WORK(KT1AMPB),WORK(KDUM)) 3311 3312 IF (LOCDEB) THEN 3313 RHO1N = DDOT(NT1AM(ISYMAB),RHO1,1,RHO1,1) 3314 RHO2N = DDOT(NT2AM(ISYMAB),RHO2,1,RHO2,1) 3315 WRITE(LUPRI,*)'Norm of RHO1 in CCMM_QRTRANS on input:,1', RHO1N 3316 WRITE(LUPRI,*)'Norm of RHO2 in CCMM_QRTRANS on input:,1', RHO2N 3317 RHO1N = DDOT(NT1AM(ISYMB),WORK(KT1AMPA),1,WORK(KT1AMPA),1) 3318 WRITE(LUPRI,*)'Norm af B1AM in CCMM_QRTRANS on input:,1', RHO1N 3319 END IF 3320 3321 CALL CCMM_TGB(WORK(KT1AMPB),ISYMB,LISTB,IDLSTB,WORK(KTGB),'ET', 3322 * MODEL,WORK(KWRK1),LWRK1) 3323C 3324 LABEL = 'GIVE INT' 3325 CALL CCCR_AA(LABEL,ISYMB,LISTA,IDLSTA, 3326 & WORK(KTGB),WORK(KWRK1),LWRK1) 3327 CALL CCLR_DIASCL(WORK(KWRK1+NT1AM(ISYMAB)),2.0d0,ISYMAB) 3328C 3329C 3330 CALL DAXPY(NT1AM(ISYMAB),ONE,WORK(KWRK1),1,RHO1,1) 3331 CALL DAXPY(NT2AM(ISYMAB),ONE, 3332 * WORK(KWRK1+NT1AM(ISYMAB)),1,RHO2,1) 3333 3334 IF (LOCDEB) THEN 3335 TAL1= DDOT(NT1AM(ISYMAB),WORK(KWRK1),1,WORK(KWRK1),1) 3336 TAL2= DDOT(NT2AM(ISYMAB),WORK(KWRK1+NT1AM(ISYMAB)),1, 3337 * WORK(KWRK1+NT1AM(ISYMAB)),1) 3338 WRITE(LUPRI,*) 'Printing second contribution. 3339 & Norm2 of singles: ', TAL1, 3340 & 'Norm2 of doubles: ', TAL2 3341 3342 TAL1= DDOT(NT1AM(ISYMAB),RHO1,1,RHO1,1) 3343 TAL2= DDOT(NT2AM(ISYMAB),RHO2,1,RHO2,1) 3344 WRITE(LUPRI,*) 'Printing RHO total. 3345 & Norm2 of singles: ', TAL1, 3346 & 'Norm2 of doubles: ', TAL2 3347 END IF 3348C 3349 END IF 3350C 3351C Finally, we concentrate on the effective operator T^gAB 3352C ------------------------------------------------------- 3353C 3354 KTGAB = 1 3355 KWRK1 = KTGAB + N2BST(ISYMAB) 3356 LWRK1 = LWORK - KWRK1 3357 IF (LWRK1.LT.0) CALL QUIT( 'Too little work in CCSL_BTR, 3') 3358C 3359 CALL DZERO(WORK(KTGAB),N2BST(ISYMAB)) 3360C 3361 LISTC = 'L0' 3362 IDLSTC = 1 3363 ISYMC = ILSTSYM(LISTC,IDLSTC) 3364C 3365C TAL2= DDOT(N2BST(ISYMAB),WORK(KTGAB),1,WORK(KTGAB),1) 3366C WRITE(LUPRI,*) 'Printing TGAB before build. 3367C & Norm2 of operator: ', TAL2 3368 CALL CCMM_TGCB(ISYMAB,LISTC,LISTA,LISTB, 3369 * IDLSTC,IDLSTA,IDLSTB, 3370 * ISYMC,ISYMA,ISYMB,WORK(KTGAB), 3371 * MODEL,WORK(KWRK1),LWRK1) 3372C 3373C TAL2= DDOT(N2BST(ISYMAB),WORK(KTGAB),1,WORK(KTGAB),1) 3374C WRITE(LUPRI,*) 'Printing TGAB after build. 3375C & Norm2 of operator: ', TAL2 3376 3377 NAMPF = NT1AM(ISYMAB) + NT2AM(ISYMAB) 3378C 3379 KXI = KWRK1 3380 KWRK2 = KXI + NAMPF 3381 LWRK2 = LWORK - KWRK2 3382 IF (LWRK2.LT.0) CALL QUIT( 'Too little work in CCMM_BTR, 4') 3383C 3384 CALL DZERO(WORK(KXI),NAMPF) 3385C 3386 LABEL = 'GIVE INT' 3387 CALL CC_XKSI(WORK(KXI),LABEL,ISYMAB,0,WORK(KTGAB), 3388 * WORK(KWRK2),LWRK2) 3389C 3390 KXI1 = KXI 3391 KXI2 = KXI1 + NT1AM(ISYMAB) 3392C 3393 CALL CCLR_DIASCL(WORK(KXI2),2.0d0,ISYMAB) 3394 3395 IF (LOCDEB) THEN 3396 TAL1= DDOT(NT1AM(ISYMAB),WORK(KXI1),1,WORK(KXI1),1) 3397 TAL2= DDOT(NT2AM(ISYMAB),WORK(KXI2),1,WORK(KXI2),1) 3398 WRITE(LUPRI,*) 'Printing XI^TGCB contribution. 3399 & Norm2 of singles: ', TAL1, 3400 & 'Norm2 of doubles: ', TAL2 3401 END IF 3402C 3403 CALL DAXPY(NT1AM(ISYMAB),ONE,WORK(KXI1),1,RHO1,1) 3404 CALL DAXPY(NT2AM(ISYMAB),ONE,WORK(KXI2),1,RHO2,1) 3405 3406 END 3407*********************************************************** 3408C /* Deck ccmm_tgcb */ 3409 SUBROUTINE CCMM_TGCB(ISYRES,LISTA,LISTB,LISTC, 3410 * IDLSTA,IDLSTB,IDLSTC, 3411 * ISYMTA,ISYMTB,ISYMTC,TGCB, 3412 * MODEL,WORK,LWORK) 3413C 3414C Constructs effective operator equal to 3415C T^gCB = - sum_a <Lambda|[[Rr_a,C],B]|CC> Rr_a 3416C JK+OC, marts 03 3417C----------------------------------------------------------------------------- 3418C 3419#include "implicit.h" 3420#include "maxorb.h" 3421#include "mxcent.h" 3422 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, IZERO = 0 , TWO = 2.0D0) 3423 PARAMETER (HALF = 0.5D0) 3424#include "dummy.h" 3425#include "iratdef.h" 3426#include "priunit.h" 3427#include "ccorb.h" 3428#include "ccsdsym.h" 3429#include "ccsdio.h" 3430#include "ccsdinp.h" 3431#include "ccfield.h" 3432#include "exeinf.h" 3433#include "ccfdgeo.h" 3434#include "ccslvinf.h" 3435#include "qm3.h" 3436#include "ccinftap.h" 3437C 3438 DIMENSION WORK(LWORK) 3439 DIMENSION TGCB(*) 3440C 3441 INTEGER KDUM 3442 INTEGER IDLSTA,IDLSTB,IDLSTC,ISYMTC,ISYMTA,ISYMTB,ISYRES 3443 PARAMETER( KDUM = +99 999 999 ) ! dummy address for work space 3444C 3445 CHARACTER*(*) LISTA,LISTB,LISTC 3446 CHARACTER MODEL*10 3447C 3448 CHARACTER*8 LABEL 3449 LOGICAL LEXIST 3450C 3451C---------------------------------------------- 3452C If all models are SPC 3453C -> RETURN from CCMM_TGB: 3454C---------------------------------------------- 3455C 3456 IF (LOSPC) RETURN 3457C 3458C-------------------------------------- 3459C Readin trial vector (B) from file 3460C-------------------------------------- 3461C 3462 KT1AMB = 1 3463 KT2AMB = KT1AMB + NT1AM(ISYMTB) 3464 KWRK1 = KT2AMB + NT2AM(ISYMTB) 3465 LWRK1 = LWORK - KWRK1 3466C 3467 IF (LWRK1 .LT. 0) THEN 3468 CALL QUIT('Insuff. work in CCMM_TGCB 1') 3469 END IF 3470C 3471 IOPT = 3 3472 CALL CC_RDRSP(LISTB,IDLSTB,ISYMTB,IOPT,MODEL, 3473 * WORK(KT1AMB),WORK(KT2AMB)) 3474C 3475C------------------ 3476C Symmetry 1 3477C------------------ 3478C 3479 ISYCB = MULD2H(ISYMTB,ISYMTC) 3480C 3481 IF (ISYCB .NE. ISYRES) THEN 3482 CALL QUIT( 'Symmetry problem in CCMM_TGCB') 3483 END IF 3484C 3485C------------------------------------ 3486C Dynamical allocation for CCMM : 3487C------------------------------------ 3488C 3489 KRAAOx = KWRK1 3490 KRAAOy = KRAAOx + N2BST(ISYCB) 3491 KRAAOz = KRAAOy + N2BST(ISYCB) 3492 KWRK2 = KRAAOz + N2BST(ISYCB) 3493 LWRK2 = LWORK - KWRK2 3494C 3495 IF (LWRK2.LT.0) CALL QUIT( 'Too little work in CCMM_CTGB, 2') 3496C 3497 CALL DZERO(WORK(KRAAOx),3*N2BST(ISYCB)) 3498C 3499C---------------------------- 3500C Read Rra in AO basis: 3501C---------------------------- 3502C 3503 LUCCPO = -1 3504 CALL GPOPEN(LUCCEF,'ELFDMM','OLD',' ', 3505 & 'UNFORMATTED',IDUMMY,.FALSE.) 3506 REWIND(LUCCEF) 3507C 3508 LM = 0 3509C 3510 DO 700 I = 1, ISYTP 3511 IF (MDLWRD(I)(1:5) .EQ. 'SPC_E') THEN 3512 DO 701 J = NSYSBG(I), NSYSED(I) 3513 DO 702 K = 1, NUALIS(I) 3514C 3515 LM = LM + 1 3516 LABEL = 'READ INT' 3517C 3518 CALL CCMMINT(LABEL,LUCCEF,WORK(KRAAOx),DUMMY,IRREP, 3519 & ISYM,IERR,WORK(KWRK2),LWRK2) 3520 CALL CCMMINT(LABEL,LUCCEF,WORK(KRAAOy),DUMMY,IRREP, 3521 & ISYM,IERR,WORK(KWRK2),LWRK2) 3522 CALL CCMMINT(LABEL,LUCCEF,WORK(KRAAOz),DUMMY,IRREP, 3523 & ISYM,IERR,WORK(KWRK2),LWRK2) 3524C 3525C------------------------------------------------------------------- 3526C 3527 IF (MDLWRD(I) .EQ. 'SPC_E01') THEN 3528C 3529 LABEL = 'GIVE INT' 3530 CALL CCLR_FA(LABEL,ISYCB,LISTC,IDLSTC, 3531 & LISTA,IDLSTA,WORK(KRAAOx),WORK(KWRK2),LWRK2) 3532C 3533 CBDOT1x = DDOT(NT1AM(ISYMTB),WORK(KWRK2), 3534 & 1,WORK(KT1AMB),1) 3535 CBDOT2x = DDOT(NT2AM(ISYMTB),WORK(KWRK2+NT1AM(ISYMTB)), 3536 & 1,WORK(KT2AMB),1) 3537 CBDOTx = CBDOT1x + CBDOT2x 3538 FACx = - ALPIMM(I,K) * (CBDOTx) 3539 CALL DAXPY(N2BST(ISYCB),FACx,WORK(KRAAOx),1, 3540 * TGCB,1) 3541C 3542 LABEL = 'GIVE INT' 3543 CALL CCLR_FA(LABEL,ISYCB,LISTC,IDLSTC, 3544 & LISTA,IDLSTA,WORK(KRAAOy),WORK(KWRK2),LWRK2) 3545C 3546 CBDOT1y = DDOT(NT1AM(ISYMTB),WORK(KWRK2), 3547 & 1,WORK(KT1AMB),1) 3548 CBDOT2y = DDOT(NT2AM(ISYMTB),WORK(KWRK2+NT1AM(ISYMTB)), 3549 & 1,WORK(KT2AMB),1) 3550 CBDOTy = CBDOT1y + CBDOT2y 3551 FACy = - ALPIMM(I,K) * (CBDOTy) 3552 CALL DAXPY(N2BST(ISYCB),FACy,WORK(KRAAOy),1, 3553 * TGCB,1) 3554C 3555 LABEL = 'GIVE INT' 3556 CALL CCLR_FA(LABEL,ISYCB,LISTC,IDLSTC, 3557 & LISTA,IDLSTA,WORK(KRAAOz),WORK(KWRK2),LWRK2) 3558C 3559 CBDOT1z = DDOT(NT1AM(ISYMTB),WORK(KWRK2), 3560 & 1,WORK(KT1AMB),1) 3561 CBDOT2z = DDOT(NT2AM(ISYMTB),WORK(KWRK2+NT1AM(ISYMTB)), 3562 & 1,WORK(KT2AMB),1) 3563 CBDOTz = CBDOT1z + CBDOT2z 3564 FACz = - ALPIMM(I,K) * (CBDOTz) 3565 CALL DAXPY(N2BST(ISYCB),FACz,WORK(KRAAOz),1, 3566 * TGCB,1) 3567C 3568 END IF 3569 702 CONTINUE 3570 701 CONTINUE 3571 END IF 3572 700 CONTINUE 3573C 3574 CALL GPCLOSE(LUCCEF,'KEEP') 3575 END 3576C**************************************************************** 3577C /* Deck ccsl_gtr */ 3578 SUBROUTINE CCMM_GTR(RHO1,RHO2,ISYRES, 3579 * LISTA,IDLSTA,ISYMTA, 3580 * LISTB,IDLSTB,ISYMTB, 3581 * LISTC,IDLSTC,ISYMTC, 3582 * MODEL,WORK,LWORK) 3583C 3584C----------------------------------------------------------------- 3585C JK+OC, marts.03 3586C Calculates the contribution to the CC/MM G matrx 3587C transformations. 3588C 3589C G_my,ny,sigma = 1/2 P(my,ny,sigma) 3590C <Lambda|[[T^g,mu,tau_ny],tau_sigma]|CC > 3591C 3592C The permutation creates 3 different contributions. 3593C 3594C------------------------------------------------------------------ 3595#include "implicit.h" 3596#include "maxorb.h" 3597#include "mxcent.h" 3598 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, IZERO = 0 , TWO = 2.0D0) 3599 PARAMETER (HALF = 0.5D0) 3600#include "dummy.h" 3601#include "iratdef.h" 3602#include "priunit.h" 3603#include "ccorb.h" 3604#include "ccsdsym.h" 3605#include "ccsdio.h" 3606#include "ccsdinp.h" 3607#include "ccfield.h" 3608#include "exeinf.h" 3609#include "ccfdgeo.h" 3610#include "ccslvinf.h" 3611#include "qm3.h" 3612#include "ccinftap.h" 3613C 3614 DIMENSION WORK(LWORK),RHO1(*),RHO2(*) 3615C 3616 CHARACTER*(*) LISTA, LISTB, LISTC 3617 CHARACTER*8 LABEL 3618 CHARACTER MODEL*10 3619 LOGICAL LEXIST, LSAME, LOCDEB 3620 PARAMETER (LOCDEB=.FALSE.) 3621 INTEGER ISYCB,ISYRES,ISYMTA,ISYMTB,ISYMTC 3622 INTEGER IDLSTA,IDLSTB,IDLSTC 3623 INTEGER FACTSL 3624C 3625 INTEGER KDUM 3626 PARAMETER( KDUM = +99 999 999 ) ! dummy address for work space 3627C 3628C---------------------------------------------- 3629C If all models are SPC 3630C -> RETURN from CCMM_TGB: 3631C---------------------------------------------- 3632C 3633 IF (LOSPC) RETURN 3634C 3635 IF (IPRINT.GT.10) THEN 3636 WRITE(LUPRI,*)'CCMM_GTR: ISYRES =',ISYRES 3637 WRITE(LUPRI,*)'CCMM_GTR: ISYMTA =',ISYMTA 3638 WRITE(LUPRI,*)'CCMM_GTR: ISYMTB =',ISYMTB 3639 WRITE(LUPRI,*)'CCMM_GTR: ISYMTC =',ISYMTC 3640 WRITE(LUPRI,*)'CCMM_GTR: LISTA =',LISTA 3641 WRITE(LUPRI,*)'CCMM_GTR: LISTB =',LISTB 3642 WRITE(LUPRI,*)'CCMM_GTR: LISTC =',LISTC 3643 WRITE(LUPRI,*)'CCMM_GTR: IDLSTA =',IDLSTA 3644 WRITE(LUPRI,*)'CCMM_GTR: IDLSTB =',IDLSTB 3645 WRITE(LUPRI,*)'CCMM_GTR: IDLSTC =',IDLSTC 3646 CALL FLSHFO(LUPRI) 3647 ENDIF 3648C 3649C Symmetry: 3650C --------- 3651C 3652 ISYCB = MULD2H(ISYMTB,ISYMTC) 3653C 3654 IF (ISYCB .NE. ISYRES) THEN 3655 CALL QUIT( 'Symmetry problem in CCMM_GTR 1') 3656 END IF 3657C 3658C 3659C First we concentrate on the effective operator T^gB 3660C ----------------------------------------------------- 3661C 3662 KT1AMPA = 1 3663 KTGB = KT1AMPA + NT1AM(ISYMTB) 3664 KWRK1 = KTGB + N2BST(ISYMTB) 3665 LWRK1 = LWORK - KWRK1 3666C 3667 IF (LWRK1 .LT. 0) THEN 3668 CALL QUIT('Insuff. work in CCMM_GTR 1') 3669 END IF 3670C 3671 CALL DZERO(WORK(KT1AMPA),NT1AM(ISYMTB)) 3672 CALL DZERO(WORK(KTGB),N2BST(ISYMTB)) 3673C 3674C Read one electron excitation part of response vector 3675C 3676 IOPT = 1 3677 CALL CC_RDRSP(LISTB,IDLSTB,ISYMTB,IOPT,MODEL, 3678 * WORK(KT1AMPA),WORK(KDUM)) 3679 3680 IF (LOCDEB) THEN 3681 RHO1N = DDOT(NT1AM(ISYRES),RHO1,1,RHO1,1) 3682 RHO2N = DDOT(NT2AM(ISYRES),RHO2,1,RHO2,1) 3683 WRITE(LUPRI,*)'Norm of RHO1 in CCMM_QRTRANS on input:,1', RHO1N 3684 WRITE(LUPRI,*)'Norm of RHO2 in CCMM_QRTRANS on input:,1', RHO2N 3685 RHO1N = DDOT(NT1AM(ISYMTB),WORK(KT1AMPA),1,WORK(KT1AMPA),1) 3686 WRITE(LUPRI,*)'Norm af B1AM in CCMM_QRTRANS on input:,1', RHO1N 3687 END IF 3688 3689 CALL CCMM_TGB(WORK(KT1AMPA),ISYMTB,LISTB,IDLSTB,WORK(KTGB),'ET', 3690 * MODEL,WORK(KWRK1),LWRK1) 3691C 3692 LABEL = 'GIVE INT' 3693 CALL CCLR_FA(LABEL,ISYMTB,LISTC,IDLSTC, 3694 & LISTA,IDLSTA,WORK(KTGB),WORK(KWRK1),LWRK1) 3695C 3696 LSAME = (LISTC.EQ.LISTB .AND. IDLSTC.EQ.IDLSTB) 3697 IF (LSAME) THEN 3698 FACTSLV = TWO 3699 ELSE 3700 FACTSLV = ONE 3701 END IF 3702C 3703 IF (LOCDEB) THEN 3704 TAL1= DDOT(NT1AM(ISYCB),WORK(KWRK1),1,WORK(KWRK1),1) 3705 TAL2= DDOT(NT2AM(ISYCB),WORK(KWRK1+NT1AM(ISYCB)),1, 3706 * WORK(KWRK1+NT1AM(ISYCB)),1) 3707 WRITE(LUPRI,*) 'Printing transformed first contribution. 3708 & Norm2 of singles: ', TAL1, 3709 & 'Norm2 of doubles: ', TAL2 3710 END IF 3711C 3712 CALL DAXPY(NT1AM(ISYCB),FACTSLV,WORK(KWRK1),1,RHO1,1) 3713 CALL DAXPY(NT2AM(ISYCB),FACTSLV, 3714 * WORK(KWRK1+NT1AM(ISYCB)),1,RHO2,1) 3715 3716 IF (LOCDEB) THEN 3717 TAL1= DDOT(NT1AM(ISYCB),RHO1,1,RHO1,1) 3718 TAL2= DDOT(NT2AM(ISYCB),RHO2,1,RHO2,1) 3719 WRITE(LUPRI,*) 'Printing RHO total. 3720 & Norm2 of singles: ', TAL1, 3721 & 'Norm2 of doubles: ', TAL2 3722 END IF 3723C 3724 IF (.NOT. (LSAME)) THEN 3725C 3726C second we consider the effective operator T^gC 3727C ----------------------------------------------------- 3728C 3729 KT1AMPA = 1 3730 KTGC = KT1AMPA + NT1AM(ISYMTC) 3731 KWRK1 = KTGC + N2BST(ISYMTC) 3732 LWRK1 = LWORK - KWRK1 3733C 3734 IF (LWRK1 .LT. 0) CALL QUIT('Insuff. work in CCMM_GTR 2') 3735C 3736 CALL DZERO(WORK(KT1AMPA),NT1AM(ISYMTC)) 3737 CALL DZERO(WORK(KTGC),N2BST(ISYMTC)) 3738C 3739C Read one electron excitation part of response vector 3740C 3741 IOPT = 1 3742 CALL CC_RDRSP(LISTC,IDLSTC,ISYMTC,IOPT,MODEL, 3743 * WORK(KT1AMPA),WORK(KDUM)) 3744 3745 IF (LOCDEB) THEN 3746 RHO1N = DDOT(NT1AM(ISYRES),RHO1,1,RHO1,1) 3747 RHO2N = DDOT(NT2AM(ISYRES),RHO2,1,RHO2,1) 3748 WRITE(LUPRI,*)'Norm of RHO1 in CCMM_QRTRANS on input:,1', RHO1N 3749 WRITE(LUPRI,*)'Norm of RHO2 in CCMM_QRTRANS on input:,1', RHO2N 3750 RHO1N = DDOT(NT1AM(ISYMTB),WORK(KT1AMPA),1,WORK(KT1AMPA),1) 3751 WRITE(LUPRI,*)'Norm af B1AM in CCMM_QRTRANS on input:,1', RHO1N 3752 END IF 3753 3754 CALL CCMM_TGB(WORK(KT1AMPA),ISYMTC,LISTC,IDLSTC,WORK(KTGC),'ET', 3755 * MODEL,WORK(KWRK1),LWRK1) 3756C 3757 LABEL = 'GIVE INT' 3758 CALL CCLR_FA(LABEL,ISYMTC,LISTB,IDLSTB, 3759 & LISTA,IDLSTA,WORK(KTGC),WORK(KWRK1),LWRK1) 3760C 3761 IF (LOCDEB) THEN 3762 TAL1= DDOT(NT1AM(ISYCB),WORK(KWRK1),1,WORK(KWRK1),1) 3763 TAL2= DDOT(NT2AM(ISYCB),WORK(KWRK1+NT1AM(ISYCB)),1, 3764 * WORK(KWRK1+NT1AM(ISYCB)),1) 3765 WRITE(LUPRI,*) 'Printing transformed second contribution. 3766 & Norm2 of singles: ', TAL1, 3767 & 'Norm2 of doubles: ', TAL2 3768 END IF 3769CC 3770 CALL DAXPY(NT1AM(ISYCB),ONE,WORK(KWRK1),1,RHO1,1) 3771 CALL DAXPY(NT2AM(ISYCB),ONE, 3772 * WORK(KWRK1+NT1AM(ISYCB)),1,RHO2,1) 3773 3774 IF (LOCDEB) THEN 3775 TAL1= DDOT(NT1AM(ISYCB),RHO1,1,RHO1,1) 3776 TAL2= DDOT(NT2AM(ISYCB),RHO2,1,RHO2,1) 3777 WRITE(LUPRI,*) 'Printing RHO total. 3778 & Norm2 of singles: ', TAL1, 3779 & 'Norm2 of doubles: ', TAL2 3780 END IF 3781C 3782 END IF 3783C finally, we consider the effective operator T^g sigma 3784C ----------------------------------------------------- 3785C 3786 KTGCB = 1 3787 KWRK1 = KTGCB + N2BST(ISYCB) 3788 LWRK1 = LWORK - KWRK1 3789 IF (LWRK1.LT.0) CALL QUIT( 'Too little work in CCMM_GTR, 3') 3790C 3791 CALL DZERO(WORK(KTGCB),N2BST(ISYCB)) 3792C 3793 CALL CCMM_TGCB(ISYRES,LISTA,LISTB,LISTC, 3794 * IDLSTA,IDLSTB,IDLSTC, 3795 * ISYMTA,ISYMTB,ISYMTC,WORK(KTGCB), 3796 * MODEL,WORK(KWRK1),LWRK1) 3797C 3798C TAL2= DDOT(N2BST(ISYRES),WORK(KTGCB),1,WORK(KTGCB),1) 3799C WRITE(LUPRI,*) 'Printing TGAB after build. 3800C & Norm2 of operator: ', TAL2 3801 NAMPF = NT1AM(ISYCB) + NT2AM(ISYCB) 3802C 3803 KETA = KWRK1 3804 KWRK2 = KETA + NAMPF 3805 LWRK2 = LWORK - KWRK2 3806 IF (LWRK2.LT.0) CALL QUIT( 'Too little work in CCMM_GTR, 4') 3807C 3808 CALL DZERO(WORK(KETA),NAMPF) 3809C 3810 LABEL = 'GIVE INT' 3811 CALL CC_ETAC(ISYCB,LABEL,WORK(KETA), 3812 * LISTA,IDLSTA,0,WORK(KTGCB),WORK(KWRK2),LWRK2) 3813C 3814 KETA1 = KETA 3815 KETA2 = KETA1 + NT1AM(ISYCB) 3816 3817 IF (LOCDEB) THEN 3818 TAL1= DDOT(NT1AM(ISYCB),WORK(KETA1),1,WORK(KETA1),1) 3819 TAL2= DDOT(NT2AM(ISYCB),WORK(KETA2),1, 3820 * WORK(KETA2),1) 3821 WRITE(LUPRI,*) 'Printing third GTR contribution. 3822 & Norm2 of singles: ', TAL1, 3823 & 'Norm2 of doubles: ', TAL2 3824 END IF 3825C 3826C 3827 CALL DAXPY(NT1AM(ISYCB),ONE,WORK(KETA1),1,RHO1,1) 3828 CALL DAXPY(NT2AM(ISYCB),ONE,WORK(KETA2),1,RHO2,1) 3829C 3830 END 3831************************************************************ 3832C /* Deck ccsl_gtr */ 3833 SUBROUTINE CCMM_XIETA(WORK,LWORK) 3834C 3835C----------------------------------------------------------------- 3836C JK, april .03 3837C 3838C Calculates the CCMM ETA and XI vectors and write them to files. 3839C 3840C------------------------------------------------------------------ 3841#include "implicit.h" 3842#include "maxorb.h" 3843#include "mxcent.h" 3844 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, IZERO = 0 , TWO = 2.0D0) 3845 PARAMETER (HALF = 0.5D0) 3846#include "dummy.h" 3847#include "iratdef.h" 3848#include "priunit.h" 3849#include "ccorb.h" 3850#include "ccsdsym.h" 3851#include "ccsdio.h" 3852#include "ccsdinp.h" 3853#include "ccfield.h" 3854#include "exeinf.h" 3855#include "ccfdgeo.h" 3856#include "ccslvinf.h" 3857#include "qm3.h" 3858#include "ccinftap.h" 3859C 3860 DIMENSION WORK(LWORK) 3861C 3862 CHARACTER*(8) FILMME, FILMMX, LABEL 3863 CHARACTER*2 LIST 3864 INTEGER LUMMET, LUMMXI, ISYMTR 3865 INTEGER KDUM 3866 PARAMETER( KDUM = +99 999 999 ) ! dummy address for work space 3867C 3868C--------------------- 3869C Init parameters. 3870C--------------------- 3871C 3872 ISYMTR = ISYMOP 3873 NTAMP1 = NT1AM(ISYMTR) 3874 NTAMP2 = NT2AM(ISYMTR) 3875 IF (CCS) NTAMP2 = 0 3876 NTAMP = NTAMP1 + NTAMP2 3877C 3878C---------------------------------------------- 3879C If all models are SPC 3880C -> RETURN from CCMM_XIETA 3881C---------------------------------------------- 3882C 3883 IF (LOSPC) RETURN 3884C 3885 IF (IPRINT.GT.10) THEN 3886 CALL FLSHFO(LUPRI) 3887 ENDIF 3888C 3889 LUMMET = -1 3890 LUMMXI = -1 3891 FILMME = 'CCMM_ETA' 3892 FILMMX = 'CCMM__XI' 3893C 3894 CALL WOPEN2(LUMMET, FILMME, 64, 0) 3895 CALL WOPEN2(LUMMXI, FILMMX, 64, 0) 3896C 3897C------------------------------------ 3898C Dynamical allocation for CCMM : 3899C------------------------------------ 3900C 3901 KRAAOx = 1 3902 KRAAOy = KRAAOx + N2BST(ISYMTR) 3903 KRAAOz = KRAAOy + N2BST(ISYMTR) 3904 KWRK1 = KRAAOz + N2BST(ISYMTR) 3905 LWRK1 = LWORK - KWRK1 3906C 3907 IF (LWRK1.LT.0) CALL QUIT( 'Too little work in CCMM_XIETA, 1') 3908C 3909 KXIx = KWRK1 3910 KXIy = KXIx + NTAMP 3911 KXIz = KXIy + NTAMP 3912 KWRK2 = KXIz + NTAMP 3913 LWRK2 = LWORK - KWRK2 3914C 3915 IF (LWRK2.LT.0) CALL QUIT( 'Too little work in CCMM_XIETA, 2') 3916C 3917 CALL DZERO(WORK(KRAAOx),3*N2BST(ISYMTR)) 3918 CALL DZERO(WORK(KXIx),3*NTAMP) 3919C 3920C---------------------------- 3921C Read Rra in AO basis: 3922C---------------------------- 3923C 3924 LUCCPO = -1 3925 CALL GPOPEN(LUCCEF,'ELFDMM','OLD',' ', 3926 & 'UNFORMATTED',IDUMMY,.FALSE.) 3927 REWIND(LUCCEF) 3928C 3929 LM = 0 3930 IADR1 = 1 3931 IADR2 = 1 3932 LEN = NTAMP 3933C 3934 DO 700 I = 1, ISYTP 3935 IF (MDLWRD(I)(1:5) .EQ. 'SPC_E') THEN 3936 DO 701 J = NSYSBG(I), NSYSED(I) 3937 DO 702 K = 1, NUALIS(I) 3938C 3939 LM = LM + 1 3940 LABEL = 'READ INT' 3941C 3942 CALL CCMMINT(LABEL,LUCCEF,WORK(KRAAOx),DUMMY,IRREP, 3943 & ISYM,IERR,WORK(KWRK2),LWRK2) 3944 CALL CCMMINT(LABEL,LUCCEF,WORK(KRAAOy),DUMMY,IRREP, 3945 & ISYM,IERR,WORK(KWRK2),LWRK2) 3946 CALL CCMMINT(LABEL,LUCCEF,WORK(KRAAOz),DUMMY,IRREP, 3947 & ISYM,IERR,WORK(KWRK2),LWRK2) 3948C 3949 IF (MDLWRD(I) .EQ. 'SPC_E01') THEN 3950 LABEL = 'GIVE INT' 3951 LIST = 'L0' 3952 IDLINO = 1 3953C 3954 CALL CC_ETAC(ISYMTR,LABEL,WORK(KXIx), 3955 * LIST,IDLINO,0,WORK(KRAAOx),WORK(KWRK2),LWRK2) 3956 CALL CC_ETAC(ISYMTR,LABEL,WORK(KXIy), 3957 * LIST,IDLINO,0,WORK(KRAAOy),WORK(KWRK2),LWRK2) 3958 CALL CC_ETAC(ISYMTR,LABEL,WORK(KXIz), 3959 * LIST,IDLINO,0,WORK(KRAAOz),WORK(KWRK2),LWRK2) 3960C 3961 CALL PUTWA2(LUMMET,FILMME,WORK(KXIx),IADR1,LEN) 3962 IADR1 = IADR1 + LEN 3963 CALL PUTWA2(LUMMET,FILMME,WORK(KXIy),IADR1,LEN) 3964 IADR1 = IADR1 + LEN 3965 CALL PUTWA2(LUMMET,FILMME,WORK(KXIz),IADR1,LEN) 3966 IADR1 = IADR1 + LEN 3967C 3968C 3969 CALL CC_XKSI(WORK(KXIx),LABEL,ISYMTR,0,WORK(KRAAOx), 3970 * WORK(KWRK2),LWRK2) 3971 CALL CC_XKSI(WORK(KXIy),LABEL,ISYMTR,0,WORK(KRAAOy), 3972 * WORK(KWRK2),LWRK2) 3973 CALL CC_XKSI(WORK(KXIz),LABEL,ISYMTR,0,WORK(KRAAOz), 3974 * WORK(KWRK2),LWRK2) 3975C 3976 CALL PUTWA2(LUMMXI,FILMMX,WORK(KXIx),IADR2,LEN) 3977 IADR2 = IADR2 + LEN 3978 CALL PUTWA2(LUMMXI,FILMMX,WORK(KXIy),IADR2,LEN) 3979 IADR2 = IADR2 + LEN 3980 CALL PUTWA2(LUMMXI,FILMMX,WORK(KXIz),IADR2,LEN) 3981 IADR2 = IADR2 + LEN 3982C 3983 END IF 3984 702 CONTINUE 3985 701 CONTINUE 3986 END IF 3987 700 CONTINUE 3988C 3989 CALL GPCLOSE(LUCCEF,'KEEP') 3990C 3991 CALL WCLOSE2(LUMMET, FILMME, 'KEEP') 3992 CALL WCLOSE2(LUMMXI, FILMMX, 'KEEP') 3993C 3994 END 3995**************************************************************** 3996C /* deck ccsl_gtr */ 3997 SUBROUTINE CCMM_REP1(WORK,LWORK) 3998C 3999C---------------------------------------------------------------- 4000C JK+OC, mai .03 4001C 4002C 4003C----------------------------------------------------------------- 4004#include "implicit.h" 4005#include "priunit.h" 4006#include "ccorb.h" 4007#include "ccsdsym.h" 4008#include "inftap.h" 4009#include "dummy.h" 4010#include "ccsdinp.h" 4011#include "ccsections.h" 4012#include "ccslvinf.h" 4013#include "qm3.h" 4014C 4015 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, IZERO = 0 , TWO = 2.0D0) 4016 DIMENSION WORK(LWORK) 4017C 4018 CHARACTER*5 FILMM 4019C 4020 KCMO = 1 4021 KEND1 = KCMO + NLAMDS 4022 LWRK1 = LWORK - KEND1 4023C 4024 IF (LWRK1 .LT. 0) THEN 4025 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 4026 CALL QUIT('Insufficient memory for allocation in CCMM_REP1 1') 4027 ENDIF 4028C 4029 IF (NSYM .NE. 1) THEN 4030 WRITE(LUPRI,*)'ONLYMO is not implemented for other than NSYM = 1' 4031 CALL QUIT('Symmetry problem in CCMM_REP1') 4032 END IF 4033C 4034C First we write to the INFO file 4035C 4036 NTAL = NRHFS(1)*NBAS(1) 4037C 4038 LUMMMO = -1 4039 CALL GPOPEN(LUMMMO,'MMMO_INFO','UNKNOWN',' ', 4040 & 'FORMATTED',IDUMMY,.FALSE.) 4041C 4042 WRITE(LUMMMO,'(I10,5X,I10,5X,I10)') NTAL, NRHFS(1), NBAS(1) 4043 CALL GPCLOSE(LUMMMO,'KEEP') 4044 4045 4046 CALL CC_GET_CMO(WORK(KCMO)) 4047 CALL CMO_REORDER(WORK(KCMO),WORK(KEND1),LWRK1) 4048C 4049 FILMM = 'MMMOS' 4050 LUMMMO = -1 4051 IADRFIL = 1 4052 NDATA = NTAL 4053C 4054 CALL WOPEN2(LUMMMO,FILMM,64,0) 4055 CALL PUTWA2(LUMMMO,FILMM,WORK(KCMO), 4056 * IADRFIL,NDATA) 4057 CALL WCLOSE2(LUMMMO,FILMM,'KEEP') 4058C 4059 KORBE = 1 4060 KEND1 = KORBE + NORBTS 4061 LWRK1 = LWORK - KEND1 4062C 4063 IF (LWRK1 .LT. 0) THEN 4064 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 4065 CALL QUIT('Insufficient memory for allocation in CCMM_REP1 2') 4066 ENDIF 4067C 4068 LUSIFC = -1 4069 CALL GPOPEN(LUSIFC,'SIRIFC','UNKNOWN',' ','UNFORMATTED',IDUMMY, 4070 & .FALSE.) 4071 REWIND LUSIFC 4072C 4073 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 4074 READ (LUSIFC) 4075 READ (LUSIFC) (WORK(KORBE+I-1), I=1,NORBTS) 4076C 4077 CALL GPCLOSE(LUSIFC,'KEEP') 4078C 4079 IF (IPRINT .GT. 20) THEN 4080 CALL AROUND('HF - Orbital energies. ') 4081 WRITE(LUPRI,*) (WORK(KORBE+I-1), I=1,NORBTS) 4082 ENDIF 4083C 4084 IF (IPRINT .GT. 5) THEN 4085 CALL AROUND('HF - Occupied Orbital energies. ') 4086 WRITE(LUPRI,*) (WORK(KORBE+I-1), I=1,NRHFS(1)) 4087 ENDIF 4088C 4089 FILMM = 'MMEPS' 4090 LUMMEP = -1 4091 IADRFIL = 1 4092 NDATA = NRHFS(1) 4093C 4094 CALL WOPEN2(LUMMEP,FILMM,64,0) 4095 CALL PUTWA2(LUMMEP,FILMM,WORK(KORBE), 4096 * IADRFIL,NDATA) 4097 CALL WCLOSE2(LUMMEP,FILMM,'KEEP') 4098C 4099C 4100 END 4101C************************************************************** 4102C /* Deck cc_qm3nsint */ 4103 SUBROUTINE CCMM_REP2(EREP,AODEN,WORK,LWORK) 4104C************************************************************** 4105C 4106C Calculates the non-coulombic repulsive contribution to the 4107C QM/MM interaction energy. 4108C JK+OC, mai 03 4109C************************************************************** 4110C 4111#include "implicit.h" 4112#include "maxorb.h" 4113#include "mxcent.h" 4114#include "nuclei.h" 4115#include "dummy.h" 4116#include "iratdef.h" 4117#include "priunit.h" 4118#include "ccorb.h" 4119#include "ccsdsym.h" 4120#include "ccsdio.h" 4121#include "ccsdinp.h" 4122#include "ccfield.h" 4123#include "exeinf.h" 4124#include "ccfdgeo.h" 4125#include "ccinftap.h" 4126#include "qm3.h" 4127C 4128 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, IZERO = 0 , TWO = 2.0D0) 4129C 4130 4131 DIMENSION WORK(LWORK) 4132 DIMENSION AODEN(*) 4133C 4134 CHARACTER*9 FILMMR 4135C 4136C----------------------------------- 4137C 4138 EREP = 0.0D0 4139C 4140 IF (.NOT. (LOSHAW)) RETURN 4141C 4142 KINT = 1 4143 KWRK1 = KINT + N2BST(ISYMOP) 4144 LWRK1 = LWORK - KWRK1 4145C 4146 IF (LWRK1.LT.0) CALL QUIT( 'Too little work in CCMM_REP2 1' ) 4147C 4148 CALL DZERO(WORK(KINT),N2BST(ISYMOP)) 4149C 4150 FILMMR = 'MMREPOP_0' 4151 LUMMRE = -1 4152 IADRFIL = 1 4153 NDATA = N2BST(ISYMOP) 4154C 4155 CALL WOPEN2(LUMMRE,FILMMR,64,0) 4156C 4157 CALL GETWA2(LUMMRE,FILMMR,WORK(KINT),IADRFIL,NDATA) 4158C 4159 IF (REPTST) THEN 4160 NBASLO = MAX(NBAST,1) 4161C 4162 KINT1 = KWRK1 4163 KINT2 = KINT1 + N2BST(ISYMOP) 4164 KWRK2 = KINT2 + N2BST(ISYMOP) 4165 LWRK2 = LWORK - KWRK2 4166C 4167 IF (LWRK2 .LT. 0) THEN 4168 CALL QUIT( 'Too litle work in CCMM_REP2 Rep. 2') 4169 END IF 4170C 4171 CALL DZERO(WORK(KINT1),2*N2BST(ISYMOP)) 4172 CALL DCOPY(N2BST(ISYMOP),WORK(KINT),1,WORK(KINT1),1) 4173C 4174 DO 668 KR=1,NREPMT 4175 CALL DGEMM('N','N',NBAST,NBAST,NBAST, 4176 * ONE,WORK(KINT),NBASLO, 4177 * WORK(KINT1),NBASLO, 4178 * ZERO,WORK(KINT2),NBASLO) 4179C 4180 CALL DCOPY(N2BST(ISYMOP),WORK(KINT2),1,WORK(KINT),1) 4181 668 CONTINUE 4182 END IF 4183C 4184 IF (IQM3PR .GT.14) THEN 4185 TAL1 = DDOT(N2BST(ISYMOP),WORK(KINT),1,WORK(KINT),1) 4186 TAL2 = SQRT(TAL1) 4187 END IF 4188C 4189 EREP = DDOT(N2BST(ISYMOP),WORK(KINT),1,AODEN,1) 4190C 4191 CALL WCLOSE2(LUMMRE,FILMMR, 'KEEP') 4192C 4193 END 4194