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 20C*********************************************************************** 21C /* Deck ccmm_addg */ 22 SUBROUTINE CCMM_ADDG(FOCK,WORK,LWORK) 23C************************************************************** 24C 25C Purpose: Construct the EFFective G solvent operator 26C and add to Fock 27C In the new QMMM implementation it replaces CCMM_RHSTG 28C 29C The general strategy of the routine is the following: 30C IF (LGFILE==FALSE) 31C 1) Get induced dipoles from file (preceeding HF/MM or CC/MM 32C calc depending on first (HF/MM) or later(CC/MM) outer iteration 33C number. 34C 2) Construct effective operator G^pol=-sum_a\mu^ind_aF^el_a (polarization 35C contribution) 36C 3) Get G^elecstat from file (multipole integrals calculated in 37C preceeding HF/MM calc) 38C ELSE IF (LFGILE==TRUE) 39C Get G operator from file dumped in the first iteration. 40C ENDIF 41C Add effective operator G to FOCK 42C 43C The LGFILE is reset in TINE CC_QMMME 44C KS oct 09 45C************************************************************** 46C 47#include "implicit.h" 48#include "maxorb.h" 49#include "mxcent.h" 50#include "qmmm.h" 51#include "dummy.h" 52#include "iratdef.h" 53#include "priunit.h" 54#include "ccorb.h" 55#include "ccsdsym.h" 56#include "ccsdio.h" 57#include "ccsdinp.h" 58#include "ccfield.h" 59#include "exeinf.h" 60#include "ccfdgeo.h" 61#include "qm3.h" 62#include "ccinftap.h" 63#include "nuclei.h" 64#include "orgcom.h" 65C 66 PARAMETER (ZERO = 0.0D0, MONE=-1.0D0 ,ONE = 1.0D0, 67 & IZERO = 0 , TWO = 2.0D0) 68 DIMENSION WORK(LWORK),FOCK(*) 69 DIMENSION INTREP(9*MXCENT), INTADR(9*MXCENT) 70 CHARACTER*8 LABINT(9*MXCENT) 71 LOGICAL TOFILE,TRIMAT,EXP1VL,LOCDEB 72 PARAMETER (LOCDEB=.FALSE.) 73C 74 CHARACTER*8 LABEL 75C 76C 77 CALL QENTER('CCMM_ADDG') 78C------------------------- 79C Init parameters 80C------------------------- 81 NDIM=3*NNZAL 82C-------------------------- 83C Dynamical allocation 84C-------------------------- 85C 86 KTAO = 1 87 KINT = KTAO + NNBASX 88 KWRK1 = KINT + N2BST(ISYMOP) 89 LWRK1 = LWORK - KWRK1 90 91 IF (LWRK1.LT.0) CALL QUIT( 'Too little work in CCMM_ADDG' ) 92 93 CALL DZERO(WORK(KTAO),NNBASX) 94 CALL DZERO(WORK(KINT),N2BST(ISYMOP)) 95C 96C 97C--------------------------------------- 98C First the polarization contribution 99C--------------------------------------- 100 101 IF (.NOT. LGFILE) THEN 102 IF (IPOLTP .GT. 0 .AND. NNZAL.NE.0) THEN 103C 104 CALL CCMM_POL(WORK(KTAO),WORK(KWRK1),LWRK1) 105C 106 END IF 107C 108C-------------------------------------- 109C Then the electrostatic contribution 110C-------------------------------------- 111 CALL CCMM_ELSTAT(WORK(KTAO),WORK(KWRK1),LWRK1) 112C 113C 114C 115 IF (IPQMMM .GT. 14) THEN 116 WRITE(LUPRI,*) 'NORM of G after electrostatic cont: ', 117 * DDOT(NNBASX,WORK(KTAO),1,WORK(KTAO),1) 118 WRITE (LUPRI,*) ' G matrix: ' 119 CALL OUTPAK(WORK(KTAO),NBAST,1,LUPRI) 120 ENDIF 121C 122C 123C-------------------------------------------------- 124C Add total contribution to effective AO fock matrix: 125C-------------------------------------------------- 126 127 LABEL= 'GIVE INT' 128 CALL CCMMINT(LABEL,-1,WORK(KINT),WORK(KTAO), 129 & IRREP,ISYM,IERR,WORK(KWRK1),LWRK1) 130C 131C Print G matrix to file 132 IF (FFIRST ) THEN !needed for PECC2 to ensure that we dump also the HF G to fil 133 CALL PUT_TO_FILE('G_MATRIXHF',N2BST(ISYMOP),WORK(KINT)) 134 FFIRST=.FALSE. 135 ELSE 136 CALL PUT_TO_FILE('G_MATRIX',N2BST(ISYMOP),WORK(KINT)) 137 138C Put LGFILE to TRUE to indicate that we read it from file in future 139C inner iterations. This keyword is reset in TINE QMMME 140 LGFILE=.TRUE. 141 ENDIF 142 143 ELSE IF (LGFILE) THEN 144 145 CALL GET_FROM_FILE('G_MATRIX',N2BST(ISYMOP),WORK(KINT)) 146 147 END IF 148 149 IF (IPQMMM .GT.14) THEN 150 CALL AROUND( 'CCMM cont. to AO matrix' ) 151 CALL OUTPUT(WORK(KINT),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI) 152 ENDIF 153C 154 IF (IPQMMM .GT.14) THEN 155 CALL AROUND( 'Usual Fock AO matrix' ) 156 CALL CC_PRFCKAO(FOCK,ISYMOP) 157 ENDIF 158C 159 CALL DAXPY(N2BST(ISYMOP),ONE,WORK(KINT),1,FOCK,1) 160C 161 IF (IPQMMM .GT.14) THEN 162 CALL AROUND( 'Modified Fock AO matrix when leaving ADDG') 163 CALL CC_PRFCKAO(FOCK,ISYMOP) 164 ENDIF 165C 166 CALL QEXIT('CCMM_ADDG') 167 END 168C 169C*********************************************************************** 170C /* Deck ccmm_addg */ 171 SUBROUTINE CCMM_ADDGHF(FOCK,WORK,LWORK) 172C************************************************************** 173C 174C Purpose: Construct the Effective G solvent operator FROM 175C THE HF INDUCED DIPOLES and add to Fock 176C 177C The general strategy of the routine is the following: 178C IF (LGFILE==FALSE) 179C 1) Get induced dipoles from file (preceeding HF/MM) 180C calc 181C 2) Construct effective operator G^pol=-sum_a\mu^ind_aF^el_a (polarization 182C contribution) 183C 3) Get G^elecstat from file (multipole integrals calculated in 184C preceeding HF/MM calc) 185C ELSE IF (LFGILE==TRUE) 186C Get G operator from file dumped in the first iteration. 187C ENDIF 188C Add effective operator G to FOCK 189C 190C The LGFILE is reset in TINE CC_QMMME 191C KS oct 09 192C************************************************************** 193C 194#include "implicit.h" 195#include "maxorb.h" 196#include "mxcent.h" 197#include "qmmm.h" 198#include "dummy.h" 199#include "iratdef.h" 200#include "priunit.h" 201#include "ccorb.h" 202#include "ccsdsym.h" 203#include "ccsdio.h" 204#include "ccsdinp.h" 205#include "ccfield.h" 206#include "exeinf.h" 207#include "ccfdgeo.h" 208#include "qm3.h" 209#include "ccinftap.h" 210#include "nuclei.h" 211#include "orgcom.h" 212C 213 PARAMETER (ZERO = 0.0D0, MONE=-1.0D0 ,ONE = 1.0D0, 214 & IZERO = 0 , TWO = 2.0D0) 215 PARAMETER (D3I = 1.0D0/3.0D0, D6I = 1.0D0/6.0D0) 216 DIMENSION WORK(LWORK),FOCK(*) 217 DIMENSION INTREP(9*MXCENT), INTADR(9*MXCENT) 218 CHARACTER*8 LABINT(9*MXCENT) 219 LOGICAL TOFILE,TRIMAT,EXP1VL,LOCDEB 220 PARAMETER (LOCDEB=.FALSE.) 221C 222 CHARACTER*8 LABEL 223C 224C 225 CALL QENTER('CCMM_ADDGHF') 226C------------------------- 227C Init parameters 228C------------------------- 229 NDIM=3*NNZAL 230C-------------------------- 231C Dynamical allocation 232C-------------------------- 233C 234 KTAO = 1 235 KINT = KTAO + NNBASX 236 KWRK1 = KINT + N2BST(ISYMOP) 237 LWRK1 = LWORK - KWRK1 238 239 IF (LWRK1.LT.0) CALL QUIT( 'Too little work in CCMM_ADDGHF' ) 240 241 CALL DZERO(WORK(KTAO),NNBASX) 242 CALL DZERO(WORK(KINT),N2BST(ISYMOP)) 243C 244C--------------------------------------- 245C First the polarization contribution 246C--------------------------------------- 247 248 IF (.NOT. LHFFIL ) THEN 249 IF (IPOLTP .GT. 0 .AND. NNZAL.NE.0) THEN 250C 251 CALL CCMM_POL(WORK(KTAO),WORK(KWRK1),LWRK1) 252C 253 END IF 254C 255C-------------------------------------- 256C Then the electrostatic contribution 257C-------------------------------------- 258 CALL CCMM_ELSTAT(WORK(KTAO),WORK(KWRK1),LWRK1) 259C 260C 261C 262 IF (IPQMMM .GT. 14 ) THEN 263 WRITE(LUPRI,*) 'NORM of G after electrostatic cont: ', 264 * DDOT(NNBASX,WORK(KTAO),1,WORK(KTAO),1) 265 WRITE (LUPRI,*) ' G matrix: ' 266 CALL OUTPAK(WORK(KTAO),NBAST,1,LUPRI) 267 ENDIF 268C-------------------------------------------------- 269C Add total contribution to effective AO fock matrix: 270C-------------------------------------------------- 271 272 LABEL= 'GIVE INT' 273 CALL CCMMINT(LABEL,-1,WORK(KINT),WORK(KTAO), 274 & IRREP,ISYM,IERR,WORK(KWRK1),LWRK1) 275C 276C Print G matrix to file 277 CALL PUT_TO_FILE('G_MATRIXHF',N2BST(ISYMOP),WORK(KINT)) 278C Put LHFFIL to TRUE to indicate that we read it from file in all future 279C inner iterations. Unlike LGFILE this keyword is never reset. 280 LHFFIL=.TRUE. 281 282 ELSE IF (LHFFIL) THEN 283 284 CALL GET_FROM_FILE('G_MATRIXHF',N2BST(ISYMOP),WORK(KINT)) 285 286 END IF 287 288 IF (IPQMMM .GT.14 ) THEN 289 CALL AROUND( 'CCMM cont. to AO matrix' ) 290 CALL OUTPUT(WORK(KINT),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI) 291 ENDIF 292C 293 IF (IPQMMM .GT.14) THEN 294 CALL AROUND( 'Usual Fock AO matrix' ) 295 CALL CC_PRFCKAO(FOCK,ISYMOP) 296 ENDIF 297C 298 CALL DAXPY(N2BST(ISYMOP),ONE,WORK(KINT),1,FOCK,1) 299C 300 IF (IPQMMM .GT.14) THEN 301 CALL AROUND( 'Modified Fock AO matrix when leaving ADDGHF') 302 CALL CC_PRFCKAO(FOCK,ISYMOP) 303 ENDIF 304C 305 306 CALL QEXIT('CCMM_ADDGHF') 307 END 308C 309C*********************************************************************** 310C 311C /* Deck ccmm_transformer */ 312 SUBROUTINE CCMM_TRANSFORMER(RHO1,RHO2,CTR1,CTR2,MODEL, 313 * ISYMTR,LR,WORK,LWORK) 314C 315C----------------------------------------------------------------------------- 316C 317C Construct effective operators from contracted densities and do 318C response transformation. The type of density and transformation 319C is controlled by the LR flag. 320C 321C IF (LR NE '0') 322C 1) Construct auxiliary density 323C 4) Construct effective operator,e.g G^C=-sum_a\∆mu^{DC}_a eps_a 324C ELSE IF (LR EQ '0') 325C Simply read GCC matrix from file constructed in the preceeding ADDG routine. 326C END IF 327C 4) Do xi/eta-transformation depending on LR 328C 5) Add to rho vector 329C 330C Sneskov, oct 09 331C----------------------------------------------------------------------------- 332C 333#include "implicit.h" 334#include "maxorb.h" 335#include "mxcent.h" 336#include "qmmm.h" 337#include "dummy.h" 338#include "iratdef.h" 339#include "priunit.h" 340#include "ccorb.h" 341#include "ccsdsym.h" 342#include "ccsdio.h" 343#include "ccsdinp.h" 344#include "ccfield.h" 345#include "exeinf.h" 346#include "ccfdgeo.h" 347#include "ccslvinf.h" 348#include "qm3.h" 349#include "ccinftap.h" 350#include "nuclei.h" 351#include "orgcom.h" 352 353 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, IZERO = 0 , TWO = 2.0D0) 354 PARAMETER (HALF = 0.5D0) 355 INTEGER NDIM 356C 357 DIMENSION WORK(LWORK),RHO1(*),RHO2(*),CTR1(*),CTR2(*) 358 CHARACTER LISDUM*(2) 359 INTEGER IDLDUM, ISYDUM 360C 361 CHARACTER MODEL*10 362 LOGICAL LOCDEB 363 PARAMETER (LOCDEB=.FALSE.) 364C 365 CHARACTER*8 LABEL,LR*(1), LIST*(2) 366C 367C 368 CALL QENTER('CCMM_TRANSFORMER') 369 370 IF (IPQMMM.GT.10) THEN 371 WRITE(LUPRI,*)'CCMM_TRANS : ISYMTR:', ISYMTR, 372 & ' and LR: ', LR 373 ENDIF 374C 375C--------------------- 376C Init parameters. 377C--------------------- 378C 379 NTAMP1 = NT1AM(ISYMTR) 380 NTAMP2 = NT2AM(ISYMTR) 381 NTAMP = NTAMP1 + NTAMP2 382 NDIM = 3*NNZAL 383C 384 IF (CCS) CALL QUIT('CCMM_TRANSFORMER: Not implemented for CCS') 385C 386 IF (DISCEX) CALL QUIT('CCMM_TRANSFORM:Not implemented for DISCEX') 387C fixed field approximation - no derivative response terms for fixed 388C ccfld 389 IF (CCFIXF .AND. (LR.NE.'0') ) THEN 390 CALL QEXIT('CCMM_TRANSFORMER') 391 RETURN 392 END IF 393C 394C------------------------------------ 395C Dynamical allocation for CCMM : 396C------------------------------------ 397C 398 KDENS = 1 399 KGMAT = KDENS + N2BST(ISYMTR) 400 KETA = KGMAT + N2BST(ISYMTR) 401 KWRK = KETA + NTAMP 402 LWRK = LWORK - KWRK 403 404 IF (LWRK .LT. 0) THEN 405 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KWRK 406 CALL QUIT( 'Too little work in CCMM_TRANSFORMER') 407 END IF 408C 409C Initialize the working matrices 410 CALL DZERO(WORK(KDENS),N2BST(ISYMTR)) 411 CALL DZERO(WORK(KETA),NTAMP) 412 CALL DZERO(WORK(KGMAT),N2BST(ISYMTR)) 413C 414C Print section 415 IF ( IPQMMM .GT. 10 .OR. LOCDEB) THEN 416 RHO1N = DDOT(NT1AM(ISYMTR),RHO1,1,RHO1,1) 417 RHO2N = DDOT(NT2AM(ISYMTR),RHO2,1,RHO2,1) 418 WRITE(LUPRI,*) ' Norm of RHO1 in CCMM_TRANS on input:', RHO1N 419 WRITE(LUPRI,*) ' Norm of RHO2 in CCMM_TRANS on input:', RHO2N 420 RHO1N = DDOT(NT1AM(ISYMTR),CTR1,1,CTR1,1) 421 RHO2N = DDOT(NT2AM(ISYMTR),CTR2,1,CTR2,1) 422 WRITE(LUPRI,*) ' Norm af C1AM in CCMM_TRANS on input:', RHO1N 423 WRITE(LUPRI,*) ' Norm af C2AM in CCMM_TRANS on input:', RHO2N 424 ENDIF 425 426 IF ( (LR.NE.'0' .AND. IPOLTP .GT. 0) .AND. (.NOT. HFFLD) ) THEN 427 428C---------------------------- 429C Construct auxiliary densities 430C---------------------------- 431 CALL CCMM_D1AO(WORK(KDENS),CTR1,CTR2,MODEL,LR, 432 * LISDUM,IDLDUM,ISYDUM, 433 * WORK(KWRK),LWRK) 434C---------------------------- 435C Construct effective G operator 436C---------------------------- 437 CALL CCMM_GEFF(WORK(KDENS),WORK(KGMAT),WORK(KWRK),LWRK) 438 439C 440 441 ELSE IF (LR .EQ. '0' ) THEN 442 IF (HFFLD ) THEN 443 CALL GET_FROM_FILE('G_MATRIXHF',N2BST(ISYMOP),WORK(KGMAT)) 444 ELSE 445 CALL GET_FROM_FILE('G_MATRIX',N2BST(ISYMOP),WORK(KGMAT)) 446 END IF 447 END IF 448 449C Transform the effective operator according to LR 450C--------------------------------------------------------------------- 451C (LR.EQ.'L').OR.(LR.EQ.'F') Calculate eta^Geff 452C (LR.EQ.'R').OR.(LR.EQ.'P').OR.(LR.EQ.'0') Calculate xi^Geff 453C--------------------------------------------------------------------- 454C 455 KETA1 = KETA 456 KETA2 = KETA1 + NTAMP1 457C 458 459 IF ( (LR.EQ.'L') .OR. (LR.EQ.'F') ) THEN 460 LABEL = 'GIVE INT' 461 LIST = 'L0' 462 IDLINO = 1 463C 464 CALL CC_ETAC(ISYMTR,LABEL,WORK(KETA), 465 * LIST,IDLINO,0,WORK(KGMAT),WORK(KWRK),LWRK) 466C 467 ELSE IF ( (LR.EQ.'R') .OR. (LR.EQ.'P') .OR. (LR.EQ.'0')) THEN 468 LABEL = 'GIVE INT' 469C 470 CALL CC_XKSI(WORK(KETA),LABEL,ISYMTR,0,WORK(KGMAT), 471 * WORK(KWRK),LWRK) 472C 473 IF (LR .EQ. 'R') THEN 474 CALL CCLR_DIASCL(WORK(KETA2),TWO,ISYMTR) 475 END IF 476 END IF 477C 478 479 IF(LOCDEB .OR. IPQMMM .GT. 14) THEN 480 TAL1= DDOT(NT1AM(ISYMTR),WORK(KETA1),1,WORK(KETA1),1) 481 TAL2= DDOT(NT2AM(ISYMTR),WORK(KETA2),1,WORK(KETA2),1) 482 WRITE(LUPRI,*) 'Printing TRANS contribution. 483 & Norm2 of singles: ', TAL1, 484 & 'Norm2 of doubles: ', TAL2 485 END IF 486 487 CALL DAXPY(NT1AM(ISYMTR),ONE,WORK(KETA1),1,RHO1,1) 488 CALL DAXPY(NT2AM(ISYMTR),ONE,WORK(KETA2),1,RHO2,1) 489C 490 IF (LOCDEB .OR. IPQMMM .GT. 14) THEN 491 TAL1= DDOT(NT1AM(ISYMTR),RHO1,1,RHO1,1) 492 TAL2= DDOT(NT2AM(ISYMTR),RHO2,1,RHO2,1) 493 WRITE(LUPRI,*) 'Printing RHO: 494 & Norm2 of singles: ', TAL1, 495 & 'Norm2 of doubles: ', TAL2 496 END IF 497C 498 CALL QEXIT('CCMM_TRANSFORMER') 499 500 END 501C 502C*********************************************************************** 503C /* Deck ccmm_d1ao */ 504 SUBROUTINE CCMM_D1AO(AODEN,CTR1,CTR2,MODEL,LR, 505 * LISTB,IDLSTB,ISYMTB, 506 * WORK,LWORK) 507C 508C Kristian Sneskov July 2009 inspired by CCX_D1AO 509C 510C Purpose: To calculate contributions to the one electron 511C auxiliary density matrix appearing in TINE CCMM_TRANSFORMER 512C and return it backtransformed to AO-basis in AODEN. 513C 514C IF(LR=R.OR.LR=F) -> Calc D^C 515C IF(LR=L.OR.LR=P) -> Calc D^B 516C Current models: CCS, CC2 and CCSD 517C KS sep 10: Modified to calculate QR density when LR=Q 518! Note that the list arguments are only dummy when LR=R,L 519C 520 USE PELIB_INTERFACE, ONLY: USE_PELIB 521#include "implicit.h" 522#include "priunit.h" 523#include "dummy.h" 524#include "maxash.h" 525#include "maxorb.h" 526#include "mxcent.h" 527#include "aovec.h" 528#include "iratdef.h" 529#include "ccorb.h" 530#include "infind.h" 531#include "blocks.h" 532#include "ccsdinp.h" 533#include "ccsdsym.h" 534#include "ccexci.h" 535#include "ccsdio.h" 536#include "distcl.h" 537#include "cbieri.h" 538#include "eribuf.h" 539#include "cclr.h" 540#include "qm3.h" 541C 542 PARAMETER (ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, TWO = 2.0D0) 543 DIMENSION AODEN(*), WORK(LWORK), CTR1(*), CTR2(*) 544 545 CHARACTER*(*) LISTB 546 INTEGER IDLSTB, ISYMTB 547 LOGICAL LOCDEB 548 PARAMETER (LOCDEB=.FALSE.) 549 CHARACTER MODEL*10,MODDUM*10 550 CHARACTER LR*1 551C 552 CALL QENTER('CCMM_D1AO') 553C 554C--------------------------------------------- 555C Find symmetry of left and right vectors. 556C--------------------------------------------- 557C 558 ISYML = ILSTSYM('L0',1) ! Symmetry of multipliers and b-trial 559 ISYMR = ILSTSYM('R0',1) ! symmetry of (t2)amplitudes 560 IF (LR .EQ. 'Q' ) ISYML=ISYMTB 561 562 IF (LR .NE. 'Q') THEN 563 IF (ISYMR .NE. ISYML) 564 & CALL QUIT('CCMM_D1AO: Density not total sym.') 565 END IF 566C 567 568C----------------------------------- 569C Initial work space allocation. 570C----------------------------------- 571C 572 KONEAI = 1 573 KONEAB = KONEAI + NT1AMX 574 KONEIJ = KONEAB + NMATAB(1) 575 KONEIA = KONEIJ + NMATIJ(1) 576 KL1AM = KONEIA + NT1AMX 577 KL2AM = KL1AM + NT1AM(ISYML) ! stores t2 when D^B and tbar2 when D^C 578 KT1AM = KL2AM + NT2SQ(ISYML) 579 KR2AM = KT1AM + NT1AM(ISYMR) 580 KEND0 = KR2AM + NT2AM(ISYMTR) 581C 582 KR2AMT = KEND0 583 KEND1 = KR2AMT + NT2AM(ISYMTR) 584 LWRK1 = LWORK - KEND1 585C 586 IF (LWRK1 .LT. 0) THEN 587 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 588 CALL QUIT('CCMM_D1AO: Not enough memory') 589 ENDIF 590C 591C-------------------------------------------- 592C Initialize one electron density arrays. 593C-------------------------------------------- 594C 595 CALL DZERO(WORK(KONEAB),NMATAB(1)) 596 CALL DZERO(WORK(KONEIJ),NMATIJ(1)) 597 CALL DZERO(WORK(KONEAI),NT1AMX) 598 CALL DZERO(WORK(KONEIA),NT1AMX) 599C-------------------------------------------- 600C Initialize working arrays. 601C NOTE if D^B -> KL1AM is just a dummy array 602C-------------------------------------------- 603 CALL DZERO(WORK(KL1AM),NT1AM(ISYML)) 604 CALL DZERO(WORK(KL2AM),NT2SQ(ISYML)) 605 CALL DZERO(WORK(KT1AM),NT1AM(ISYMR)) 606 CALL DZERO(WORK(KR2AM),NT2AM(ISYMTR)) 607 CALL DZERO(WORK(KR2AMT),NT2AM(ISYMTR)) 608C 609C---------------------- 610C Read left vector. 611C---------------------- 612C 613C 614! DH: bugfix, .OR.(LR.EQ.'Q') was missing 615! IF((LR.EQ.'R').OR.(LR.EQ.'F')) THEN ! old IF-THEN 616 IF((LR.EQ.'R').OR.(LR.EQ.'F').OR.(LR.EQ.'Q')) THEN 617 IOPT = 3 618 IF (LR.EQ.'Q') THEN 619 CALL CC_RDRSP(LISTB,IDLSTB,ISYMTB,IOPT,MODDUM, 620 & WORK(KL1AM),WORK(KR2AMT)) 621 ELSE 622 CALL CC_RDRSP('L0',1,ISYML,IOPT,MODDUM,WORK(KL1AM), 623 & WORK(KR2AMT)) 624 END IF 625C Read in t1 amplitudes 626 IOPT = 1 627 CALL CC_RDRSP('R0',0,1,IOPT,MODEL,WORK(KT1AM),DUMMY) 628C Copy CTR2 to R2AM since we will scale diagonal 629 CALL DCOPY(NT2AM(ISYMR),CTR2,1,WORK(KR2AM),1) 630 631 IF (ISYMR.EQ.1) CALL CCLR_DIASCL(WORK(KR2AM),TWO,ISYMR) 632 IF (LOCDEB) THEN 633 TAL1 = DDOT(NT1AMX,WORK(KL1AM),1,WORK(KL1AM),1) 634 TAL2 = DDOT(NT2AMX,WORK(KR2AMT),1,WORK(KR2AMT),1) 635 WRITE(LUPRI,*) 'CCMM_TGDEN: Norm of L1: ', TAL1, 636 & 'Norm of L2: ', TAL2 637 END IF 638 ELSE IF ( (LR .EQ. 'L') .OR. (LR .EQ. 'P') ) THEN 639 CALL DCOPY(NT2AM(ISYMTR),CTR2,1,WORK(KR2AMT),1) 640 IOPT = 3 641 CALL CC_RDRSP('R0',1,ISYMR,IOPT,MODDUM, 642 * WORK(KT1AM),WORK(KR2AM)) 643 IF (LOCDEB) THEN 644 TAL1 = DDOT(NT1AMX,WORK(KT1AM),1,WORK(KT1AM),1) 645 TAL2 = DDOT(NT2AMX,WORK(KR2AM),1,WORK(KR2AM),1) 646 WRITE(LUPRI,*) 'CCMM_TGDEN: Norm of T1: ', TAL1, 647 & 'Norm of T2: ', TAL2 648 END IF 649 END IF 650C 651C---------------------- 652C Square up tbar_2 multipliers or b_2 elements 653C---------------------- 654C 655 CALL CC_T2SQ(WORK(KR2AMT),WORK(KL2AM),ISYML) 656C 657C--------------------------------------- 658C Set up 2C-E of cluster amplitudes. 659C--------------------------------------- 660C 661C Note that diascl amplitudes are expected 662 IF (.NOT. MP2) THEN 663 CALL DCOPY(NT2AM(ISYMR),WORK(KR2AM),1,WORK(KR2AMT),1) 664 IOPTTCME = 1 665 CALL CCSD_TCMEPK(WORK(KR2AMT),1.0D0,ISYMR,IOPTTCME) 666 ENDIF 667C 668C------------------------------- 669C Work space allocation one. 670C Note that D(ia) is stored transposed 671C------------------------------- 672C 673 KXMAT = KEND1 674 KYMAT = KXMAT + NMATIJ(1) 675 KEND2 = KYMAT + NMATAB(1) 676 LWRK2 = LWORK - KEND1 677C 678 IF (LWRK2 .LT. 0) THEN 679 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 680 CALL QUIT('Insufficient memory for 2. alloc. in CCMM_D1AO') 681 ENDIF 682 CALL DZERO(WORK(KXMAT),NMATIJ(1)) 683 CALL DZERO(WORK(KYMAT),NMATAB(1)) 684C 685C----------------------------------------------------------- 686C Calculate X-intermediate of L2AM- and R2AM-amplitudes. 687C----------------------------------------------------------- 688C Note that KL2AM must be squared when input in these routines 689C 690 CALL CC_XI(WORK(KXMAT),WORK(KL2AM),ISYML,WORK(KR2AM),ISYMR, 691 * WORK(KEND2),LWRK2) 692C 693C----------------------------------------------------------- 694C Calculate Y-intermediate of L2AM- and R2AM-amplitudes. 695C----------------------------------------------------------- 696C 697 CALL CC_YI(WORK(KYMAT),WORK(KL2AM),ISYML,WORK(KR2AM),ISYMR, 698 * WORK(KEND2),LWRK2) 699C 700C-------------------------------------------------------------- 701C Construct <L2|[Emn,R2]|HF> contribution to DXaa and DXii. 702C-------------------------------------------------------------- 703C 704 CALL DCOPY(NMATAB(1),WORK(KYMAT),1,WORK(KONEAB),1) 705 CALL CC_EITR(WORK(KONEAB),WORK(KONEIJ),WORK(KEND2),LWRK2,1) 706 CALL DAXPY(NMATIJ(1),-ONE,WORK(KXMAT),1,WORK(KONEIJ),1) 707C 708 IF ( (LR .EQ. 'R') .OR. (LR .EQ. 'F') .OR. (LR .EQ. 'Q')) THEN 709C 710C-------------------------------------------------------------- 711C Construct <HF|[Emn,R1]|HF> contribution. 712C-------------------------------------------------------------- 713C 714 IF ( (LR .EQ. 'R') .OR. (LR .EQ. 'F') ) THEN 715 CALL DAXPY(NT1AMX,TWO,CTR1,1,WORK(KONEIA),1) 716 END IF 717C 718C 719C-------------------------------------------------------------- 720C Construct <L1|[Emn,R1]|HF> contribution to DXaa and DXii. 721C-------------------------------------------------------------- 722 CALL CC_DXIJ(WORK(KONEIJ),WORK(KL1AM),ISYML,CTR1,ISYMR) 723 CALL CC_DXAB(WORK(KONEAB),WORK(KL1AM),ISYML,CTR1,ISYMR) 724C 725C-------------------------------------------------------------- 726C Construct <L1|[Eia,R2]|HF> contribution to DXia 727C-------------------------------------------------------------- 728C 729 CALL CC_DXIA12(WORK(KONEIA),WORK(KL1AM),ISYML, 730 * WORK(KR2AMT),ISYMR) 731C 732C---------------------------------------------------------- 733C Construct <L2|[[Eia,R1],T2]|HF> contribution to DXia. 734C---------------------------------------------------------- 735 KT2AM = KEND0 736 KEND3 = KT2AM + NT2AM(ISYMOP) 737 LWRK3 = LWORK - KEND3 738C read t2 amplitudes here rather than in CC_DXIA21 739 IOPT = 2 740 CALL CC_RDRSP('R0',0,1,IOPT,MODEL,DUMMY,WORK(KT2AM)) 741C 742 743 IF (LWRK3 .LT. 0) THEN 744 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 745 CALL QUIT('Insufficient memory for 3. alloc. in CCMM_D1AO') 746 ENDIF 747 CALL CC_DXIA21(WORK(KONEIA),WORK(KL2AM),ISYML, 748 * CTR1,ISYMR,WORK(KT2AM),ISYMR, 749 * WORK(KEND3),LWRK3) 750 751 ELSE IF ((LR .EQ. 'L') .OR. (LR .EQ. 'P')) THEN 752C 753C---------------------------------------------------------- 754C Construct <B1|[Eai|HF> contribution to DXai. 755C---------------------------------------------------------- 756 CALL DAXPY(NT1AMX,ONE,CTR1,1,WORK(KONEAI),1) 757 758C-------------------------------------------------------------- 759C Construct <B1|[Eia,T2]|HF> contribution to DXia 760C-------------------------------------------------------------- 761 CALL CC_DXIA12(WORK(KONEIA),CTR1,ISYMTR,WORK(KR2AMT),ISYMR) 762 END IF 763C 764C-------------------------- 765C Write out test norms. 766C-------------------------- 767C 768 IF (LOCDEB) THEN 769 XIJ = DDOT(NMATIJ(1),WORK(KONEIJ),1,WORK(KONEIJ),1) 770 XAB = DDOT(NMATAB(1),WORK(KONEAB),1,WORK(KONEAB),1) 771 XAI = DDOT(NT1AMX,WORK(KONEAI),1,WORK(KONEAI),1) 772 XIA = DDOT(NT1AMX,WORK(KONEIA),1,WORK(KONEIA),1) 773 WRITE(LUPRI,*) 'Norms: DXIJ',XIJ 774 WRITE(LUPRI,*) 'Norms: DXAB',XAB 775 WRITE(LUPRI,*) 'Norms: DXAI',XAI 776 WRITE(LUPRI,*) 'Norms: DXIA',XIA 777 ENDIF 778C 779C---------------------------------- 780C Calculate the lambda matrices. 781C---------------------------------- 782C 783 KLAMDP = KEND0 784 KLAMDH = KLAMDP + NLAMDT 785 KEND4 = KLAMDH + NLAMDT 786 LWRK4 = LWORK - KEND4 787 IF (LWRK4 .LT. 0) THEN 788 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 789 CALL QUIT('Insufficient memory for 4. alloc. in CCMM_D1AO') 790 ENDIF 791 CALL LAMMAT(WORK(KLAMDP),WORK(KLAMDH),WORK(KT1AM),WORK(KEND4), 792 * LWRK4) 793C 794C-------------------------------------------------------- 795C Add the one electron density in the AO-basis. 796C-------------------------------------------------------- 797CC 798 ISDEN = 1 799 CALL CC_DENAO(AODEN,ISDEN,WORK(KONEAI),WORK(KONEAB), 800 * WORK(KONEIJ),WORK(KONEIA),ISDEN,WORK(KLAMDP),1, 801 * WORK(KLAMDH),1,WORK(KEND4),LWRK4) 802C 803 CALL QEXIT('CCMM_D1AO') 804C 805 1 FORMAT(1x,A35,1X,E20.10) 806 RETURN 807 END 808C 809C*********************************************************************** 810C /* Deck ccmm_f2dip */ 811 SUBROUTINE CCMM_F2DIP(ELF,INDMOM,WORK,LWORK) 812C 813C April 2010, KS 814C 815C Purpose: Transform electric field to induced dipole 816C by either direct inversion (MMMAT) or by iterative means (MMITER). 817C 818#include "implicit.h" 819#include "maxorb.h" 820#include "mxcent.h" 821#include "dummy.h" 822#include "iratdef.h" 823#include "priunit.h" 824#include "ccorb.h" 825#include "ccsdsym.h" 826#include "ccsdio.h" 827#include "ccsdinp.h" 828#include "ccfield.h" 829#include "exeinf.h" 830#include "ccfdgeo.h" 831#include "ccslvinf.h" 832#include "ccinftap.h" 833#include "nuclei.h" 834#include "orgcom.h" 835#include "qm3.h" 836#include "qmmm.h" 837C 838 DIMENSION WORK(LWORK) 839 DOUBLE PRECISION ELF(*),INDMOM(*) 840 LOGICAL LOCDEB,FNDLAB 841 PARAMETER (LOCDEB=.FALSE.) 842 PARAMETER (ZERO=0.0D0, ONE=1.0D0) 843 INTEGER NDIM 844 845 CALL QENTER('CCMM_F2DIP') 846 847 NDIM=3*NNZAL 848 NDIM2P = NDIM * (NDIM+1)/2!packed response matrix 849 850 KRELAY=1 851 KWRK1=KRELAY+NDIM2P 852 LWRK1 =LWORK-KWRK1 853 IF (LWRK1.LT.0) CALL QUIT('Too little work in CCMM_F2DIP') 854 855 CALL DZERO(WORK(KRELAY),NDIM2P) 856C 857C 858 IF (MMMAT) THEN 859 LUQMMM = -1 860 IF(LUQMMM .LT. 0) THEN 861 CALL GPOPEN(LUQMMM,'QMMMIM','UNKNOWN','SEQUENTIAL', 862 & 'UNFORMATTED',IDUMMY,.FALSE.) 863 ENDIF 864 REWIND(LUQMMM) 865C 866 IF (FNDLAB('QQMMMMAT',LUQMMM)) THEN 867 CALL READT(LUQMMM,NDIM2P,WORK(KRELAY)) 868 ELSE 869 CALL QUIT('Problem reading the relay 870 & matrix from QMMMIM file') 871 ENDIF 872 873 CALL GPCLOSE(LUQMMM,'KEEP') 874C 875 IF (IPQMMM .GE. 5) THEN 876 WRITE(LUPRI,*) 'The relay matrix is read from file' 877 CALL OUTPAK(WORK(KRELAY),NDIM,1,LUPRI) 878 ENDIF 879C 880 IF (IPRINT .GT. 1) THEN 881 WRITE(LUPRI,*) 882 WRITE(LUPRI,1051) 883 WRITE(LUPRI,1050) 884 WRITE(LUPRI,1051) 885 WRITE(LUPRI,*) 886 ENDIF 887C 888 IF(LOCDEB) THEN 889 WRITE(LUPRI,*) 'F vector to be transformed to dipole: ' 890 DO I=1,NDIM 891 WRITE(LUPRI,*) ELF(I) 892 END DO 893 END IF 894 895C 896C calculate the corresponding induced moments 897 CALL DSPMV('L', NDIM, ONE, WORK(KRELAY), ELF, 1,ZERO, INDMOM,1) 898 899C 900 ELSE IF (MMITER) THEN 901C convert the F vector into induced moments in an iterative fashion 902 IOPT=1 903 CALL F2QMMM(ELF,NNZAL,INDMOM,WORK(KRELAY), 904 * LWORK,IOPT,IPQMMM) 905C 906 IF (IPQMMM .GT. 1) THEN 907 WRITE(LUPRI,*) 908 WRITE(LUPRI,1030) 909 WRITE(LUPRI,*) 910 WRITE(LUPRI,1000) 911 WRITE(LUPRI,1010) 912 WRITE(LUPRI,1000) 913 ENDIF 914C 915 IINIM = 0 916 917 DO I=1,MMCENT 918 IF (ZEROAL(I) .EQ. -1) THEN 919 DIPX = 0.0D0 920 DIPY = 0.0D0 921 DIPZ = 0.0D0 922 ELSE 923 DIPX = INDMOM(IINIM) 924 DIPY = INDMOM(IINIM +1) 925 DIPZ = INDMOM(IINIM +2) 926 IINIM = IINIM +3 927 ENDIF 928 IF (IPQMMM .GT. 1) THEN 929 WRITE(LUPRI,1020) I, DIPX, DIPY, DIPZ 930 END IF 931 END DO 932 933 IF(IPQMMM .GT. 1) THEN 934 WRITE(LUPRI,1000) 935 WRITE(LUPRI,*) 936 ENDIF 937 END IF 938C 939C 940 1050 FORMAT(' Induced dipole moments ') 941 1051 FORMAT(2X,'=',22('-'),'=',2X) 942 1000 FORMAT(1X,51('=')) 943 1010 FORMAT(' | Site | X | Y | Z |') 944 1020 FORMAT(1X,I6,3(4X,F10.6)) 945 1030 FORMAT(' Total induced dipole moments: ') 946 947 CALL QEXIT('CCMM_F2DIP') 948 END 949CC 950C*********************************************************************** 951C /* Deck ccmm_elstat */ 952 SUBROUTINE CCMM_ELSTAT(NSAO,WORK,LWORK) 953C 954C April 2010, KS 955C 956C Purpose: Calculate electrostatic contribution to G operator. In 957C old qm3 formulation we refered to this contribution 958C as Ns and it is included in the 959C effective operator (G) by setting OLDTG=TRUE 960C We get the multipole integrals from the HF/MM 961C calculation available from file. 962C We do all the multipole corrections simultaneosly. We will analyze 963C the HF/MM results if we want to consider each multipole contribution 964C separately. 965C 966C NSAO is a vector containing the sum (of all orders) of 967C the multipole integrals. This is the one we need at output. Note that 968C it is returned packed i.e. LENGTH=NNBASX 969C Note that all signs are carried in the integrals themselves (cf. 970C sirqmmm.F) and daxpy is always with a + 971C 972#include "implicit.h" 973#include "maxorb.h" 974#include "mxcent.h" 975#include "dummy.h" 976#include "iratdef.h" 977#include "priunit.h" 978#include "ccorb.h" 979#include "ccsdsym.h" 980#include "ccsdio.h" 981#include "ccsdinp.h" 982#include "ccfield.h" 983#include "exeinf.h" 984#include "ccfdgeo.h" 985#include "ccslvinf.h" 986#include "ccinftap.h" 987#include "nuclei.h" 988#include "orgcom.h" 989#include "qm3.h" 990#include "qmmm.h" 991C 992 DOUBLE PRECISION NSAO(*) 993 DIMENSION WORK(LWORK) 994 PARAMETER (ONE=1.0D0) 995 PARAMETER (MONE=-1.0D0) 996 CALL QENTER('CCMM_ELSTAT') 997 998 KMULTI = 1 999 KWRK = KMULTI + NNBASX 1000 LWRK = LWORK - KWRK 1001 IF (LWRK.LT.0) CALL QUIT('Too little work in CCMM_ELSTAT') 1002 CALL DZERO(WORK(KMULTI),NNBASX) 1003 1004 IF (NMULT .GE. 0) THEN 1005C 1006 LUQMMM = -1 1007 CALL GPOPEN(LUQMMM,'MU0INT','OLD',' ', 1008 & 'UNFORMATTED',IDUMMY,.FALSE.) 1009 REWIND (LUQMMM) 1010 CALL READT(LUQMMM,NNBASX,WORK(KMULTI)) 1011 CALL GPCLOSE(LUQMMM,'KEEP') 1012 IF (IPQMMM .GT. 14) THEN 1013 WRITE(LUPRI,*) 'CCMM_ELSTAT: Print the stored 1014 & charge integrals:' 1015 CALL OUTPAK(WORK(KMULTI),NBAST,1,LUPRI) 1016 END IF 1017 1018 CALL DAXPY(NNBASX,ONE,WORK(KMULTI),1,NSAO,1) 1019 1020 END IF 1021C 2) The dipole contribution 1022 IF (NMULT .GE. 1) THEN 1023 CALL DZERO(WORK(KMULTI),NNBASX) 1024 CALL GPOPEN(LUQMMM,'MU1INT','OLD',' ', 1025 & 'UNFORMATTED',IDUMMY,.FALSE.) 1026 REWIND (LUQMMM) 1027 CALL READT(LUQMMM,NNBASX,WORK(KMULTI)) 1028 CALL GPCLOSE(LUQMMM,'KEEP') 1029 IF (IPQMMM .GT. 14) THEN 1030 WRITE(LUPRI,*) 'CCMM_ELSTAT: Print the stored 1031 & dipole integrals:' 1032 CALL OUTPAK(WORK(KMULTI),NBAST,1,LUPRI) 1033 END IF 1034 1035 CALL DAXPY(NNBASX,ONE,WORK(KMULTI),1,NSAO,1) 1036 1037 END IF 1038C 3) The quadrupole contribution 1039 IF (NMULT .GE. 2) THEN 1040 CALL DZERO(WORK(KMULTI),NNBASX) 1041 CALL GPOPEN(LUQMMM,'MU2INT','OLD',' ', 1042 & 'UNFORMATTED',IDUMMY,.FALSE.) 1043 REWIND (LUQMMM) 1044 CALL READT(LUQMMM,NNBASX,WORK(KMULTI)) 1045 CALL GPCLOSE(LUQMMM,'KEEP') 1046 IF (IPQMMM .GT. 14) THEN 1047 WRITE(LUPRI,*) 'CCMM_ELSTAT: Print the stored 1048 & quadrupole integrals:' 1049 CALL OUTPAK(WORK(KMULTI),NBAST,1,LUPRI) 1050 END IF 1051 1052 CALL DAXPY(NNBASX,ONE,WORK(KMULTI),1,NSAO,1) 1053 1054 END IF 1055C 1056 1057 CALL QEXIT('CCMM_ELSTAT') 1058 END 1059C*********************************************************************** 1060C /* Deck ccmm_dens2f */ 1061 SUBROUTINE CCMM_DENS2F(DENS,ELF,WORK,LWORK) 1062C 1063C April 2010, KS 1064C 1065C Purpose: Transform aux density to electric field. 1066C As of now (april-10) it only handles the aux density 1067C and not the construction of the normal induced dipoles. 1068C These are instead handled in FSTATIC 1069C KEPSAO contains the packed electric field integrals. 1070C KEPSx/y/z contains the unpacked field integrals. 1071C 1072#include "implicit.h" 1073#include "maxorb.h" 1074#include "mxcent.h" 1075#include "dummy.h" 1076#include "iratdef.h" 1077#include "priunit.h" 1078#include "ccorb.h" 1079#include "ccsdsym.h" 1080#include "ccsdio.h" 1081#include "ccsdinp.h" 1082#include "ccfield.h" 1083#include "exeinf.h" 1084#include "ccfdgeo.h" 1085#include "ccslvinf.h" 1086#include "ccinftap.h" 1087#include "nuclei.h" 1088#include "orgcom.h" 1089#include "qm3.h" 1090#include "qmmm.h" 1091#include "cclr.h" 1092C 1093 DIMENSION WORK(LWORK), DENS(*), ELF(*) 1094 INTEGER NDIM 1095 LOGICAL LOCDEB, LSKIP 1096 CHARACTER*8 LABEL 1097C 1098 CALL QENTER('CCMM_DENS2F') 1099C 1100 NDIM=3*NNZAL 1101 1102C Memory allocation 1103 KEPSAO = 1 1104 KEPSx = KEPSAO + 3*NNBASX 1105 KEPSy = KEPSx + N2BST(ISYMTR) 1106 KEPSz = KEPSy + N2BST(ISYMTR) 1107 KWRK = KEPSz + N2BST(ISYMTR) 1108 LWRK = LWORK - KWRK 1109 IF (LWRK.LT.0) CALL QUIT('Too little work in CCMM_DENS2F') 1110C 1111 1112 LRI=1 1113C 1114 DO 200 I = 1,MMCENT 1115 1116 LSKIP = .FALSE. 1117 1118 CALL CCMM_EPSAO(WORK(KEPSAO),I,LSKIP,WORK(KWRK),LWRK) 1119 1120 IF (LSKIP) THEN 1121 IF (ZEROAL(I) .EQ. -1) THEN 1122 GOTO 200 1123 ELSE 1124 LRI = LRI + 3 1125 GOTO 200 1126 END IF 1127 END IF 1128 1129C Unpack the integrals 1130 LABEL = 'GIVE INT' 1131 1132 CALL DZERO(WORK(KEPSx),3*N2BST(ISYMTR)) 1133 1134 CALL CCMMINT(LABEL,LUQMMM,WORK(KEPSx),WORK(KEPSAO), 1135 & IRREP,ISYM,IERR,WORK(KWRK),LWRK) 1136 CALL CCMMINT(LABEL,LUQMMM,WORK(KEPSy),WORK(KEPSAO+NNBASX), 1137 & IRREP,ISYM,IERR,WORK(KWRK),LWRK) 1138 CALL CCMMINT(LABEL,LUQMMM,WORK(KEPSz),WORK(KEPSAO+2*NNBASX), 1139 & IRREP,ISYM,IERR,WORK(KWRK),LWRK) 1140C 1141C Do the trace operation of the product matrices 1142C 1143 PRODx = DDOT(N2BST(ISYMTR),WORK(KEPSx),1,DENS,1) 1144 PRODy = DDOT(N2BST(ISYMTR),WORK(KEPSy),1,DENS,1) 1145 PRODz = DDOT(N2BST(ISYMTR),WORK(KEPSz),1,DENS,1) 1146 1147 ELF(LRI + 0) = PRODx 1148 ELF(LRI + 1) = PRODy 1149 ELF(LRI + 2) = PRODz 1150 1151 LRI = LRI + 3 1152 1153 200 CONTINUE 1154 1155 1156 IF (IPQMMM .GE. 5) THEN 1157 WRITE(LUPRI,*) 'Transformed electric fields and NDIM', NDIM 1158 CALL OUTPUT(ELF,0,NDIM,1,1,NDIM,1,1,LUPRI) 1159 WRITE(LUPRI,*) 1160 END IF 1161 CALL QEXIT('CCMM_DENS2F') 1162 1163 END 1164C*********************************************************************** 1165C /* Deck ccmm_pol */ 1166 SUBROUTINE CCMM_POL(GPOL,WORK,LWORK) 1167C 1168C April 2010, KS 1169C 1170C Purpose: Polarization contribution to effective operator G 1171C The first time we enter this routine we read the converged 1172C induced dipoles from the preceeding HF/MM calculation. 1173C In later iterations we read the induced dipoles written in 1174C CC_QMMME 1175C G=-sum_a mu^ind_a eps_a 1176C 1177C 1178#include "implicit.h" 1179#include "maxorb.h" 1180#include "mxcent.h" 1181#include "dummy.h" 1182#include "iratdef.h" 1183#include "priunit.h" 1184#include "ccorb.h" 1185#include "ccsdsym.h" 1186#include "ccsdio.h" 1187#include "ccsdinp.h" 1188#include "ccfield.h" 1189#include "exeinf.h" 1190#include "ccfdgeo.h" 1191#include "ccslvinf.h" 1192#include "ccinftap.h" 1193#include "nuclei.h" 1194#include "orgcom.h" 1195#include "qm3.h" 1196#include "qmmm.h" 1197#include "cclr.h" 1198C 1199 DIMENSION WORK(LWORK) 1200 DOUBLE PRECISION GPOL(*) 1201 INTEGER NDIM 1202 LOGICAL LSKIP 1203 1204 CALL QENTER('CCMM_POL') 1205 1206 NDIM=3*NNZAL 1207 1208 KEPSAO = 1 1209 KINDIP = KEPSAO + 3*NNBASX 1210 KWRK = KINDIP + NDIM 1211 LWRK = LWORK - KWRK 1212 1213 IF (LWRK.LT.0) CALL QUIT('Too little work in CCMM_POL') 1214 CALL DZERO(WORK(KINDIP),NDIM) 1215 1216C Read the induced moments from file 1217 CALL GET_FROM_FILE_1('INDUCED_DIPOLES',NNZAL,WORK(KINDIP)) 1218C 1219C------------------- 1220C Print section: 1221C------------------- 1222C 1223 1224 IF (IPQMMM .GT. 10) THEN 1225 WRITE(LUPRI,'(//,A)') 1226 * ' +============================================' 1227 * //'==============+' 1228 WRITE(LUPRI,'(A)') 1229 * ' | MMcent | <mu^ind>_x | <mu^ind>_y |' 1230 * //' <mu^ind>_z |' 1231 WRITE(LUPRI,'(1X,A)') 1232 * '+============================================' 1233 * //'==============+' 1234 1235 LRI = 0 1236 1237 DO 100 I = 1, MMCENT 1238 1239 IF (ZEROAL(I) .EQ. -1) GOTO 100 1240 1241 1242 WRITE(LUPRI,'(1X,A,I3,A,F16.10,A, 1243 & F16.10,A,F16.10,A)') 1244 * '| ', I,'|', WORK(KINDIP+LRI+0),' |', 1245 * WORK(KINDIP+LRI+1),' |', 1246 * WORK(KINDIP+LRI+2),' |' 1247 WRITE(LUPRI,'(1X,A)') 1248 * '+---------------------------------------------' 1249 * //'-------------+' 1250 LRI = LRI + 3 1251 100 CONTINUE 1252 WRITE(LUPRI,'(//,A)') 1253 END IF 1254 1255C 1256C--------------------------------------- 1257C Calculate electric field due to electrons at s in ao basis 1258C--------------------------------------- 1259 1260 LRI = 0 1261C 1262 DO 200 I = 1, MMCENT 1263 1264 LSKIP =.FALSE. 1265 1266 CALL CCMM_EPSAO(WORK(KEPSAO),I,LSKIP,WORK(KWRK),LWRK) 1267 1268 IF (LSKIP) THEN 1269 IF (ZEROAL(I) .EQ. -1) THEN 1270 GOTO 200 1271 ELSE 1272 LRI = LRI + 3 1273 GOTO 200 1274 END IF 1275 END IF 1276 1277 FACx = - WORK(KINDIP + LRI +0 ) 1278 FACy = - WORK(KINDIP + LRI +1 ) 1279 FACz = - WORK(KINDIP + LRI +2 ) 1280C 1281 CALL DAXPY(NNBASX,FACx,WORK(KEPSAO), 1282 * 1,GPOL,1) 1283C 1284 CALL DAXPY(NNBASX,FACy,WORK(KEPSAO+NNBASX), 1285 * 1,GPOL,1) 1286C 1287 CALL DAXPY(NNBASX,FACz,WORK(KEPSAO+2*NNBASX), 1288 * 1,GPOL,1) 1289 1290 LRI = LRI + 3 1291 200 CONTINUE 1292C 1293 IF (IPQMMM.GT.14) THEN 1294 WRITE (LUPRI,*) ' NORM of G matrix after pol. contr.: ', 1295 * DDOT(NNBASX,GPOL,1,GPOL,1) 1296 WRITE (LUPRI,'(/A)') ' G matrix: ' 1297 CALL OUTPAK(GPOL,NBAST,1,LUPRI) 1298 ENDIF 1299 1300 CALL QEXIT('CCMM_POL') 1301 END 1302 1303C*********************************************************************** 1304C /* Deck ccmm_fstatic */ 1305 SUBROUTINE CCMM_FSTATIC(ELF,AODEN,WORK,LWORK) 1306C 1307C April 2010, KS 1308C 1309C Purpose: Calculate the total static field as 1310C F^sta=F^mul+F^el+F^nuc 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 "ccslvinf.h" 1327#include "ccinftap.h" 1328#include "nuclei.h" 1329#include "orgcom.h" 1330#include "qm3.h" 1331#include "qmmm.h" 1332#include "cclr.h" 1333C 1334 DIMENSION WORK(LWORK) 1335 DIMENSION ELF(*),AODEN(*) 1336 INTEGER NDIM 1337 LOGICAL LOCDEB, LSKIP 1338 PARAMETER (LOCDEB=.FALSE.) 1339 CHARACTER LABEL*8 1340 1341 CALL QENTER('CCMM_FSTATIC') 1342 1343 NDIM=3*NNZAL 1344 ISYMOP=1 1345 1346 KEPSAO=1 1347 KEPSx=KEPSAO+3*NNBASX 1348 KEPSy=KEPSx+N2BST(ISYMOP) 1349 KEPSz=KEPSy+N2BST(ISYMOP) 1350 KWRK=KEPSz+N2BST(ISYMOP) 1351 LWRK=LWORK-KWRK 1352 1353 LRI = 1 ! Row-index in supervector 1354 1355 DO 100 I = 1,MMCENT 1356 1357C First we calculate the field due to the QM electrons by calculating 1358C the electric field operator and dotting with the AO CC density to 1359C obtain the CC expectation value 1360 LSKIP = .FALSE. 1361 1362 CALL CCMM_EPSAO(WORK(KEPSAO),I,LSKIP,WORK(KWRK),LWRK) 1363 1364 IF (LSKIP) THEN 1365 IF (ZEROAL(I) .EQ. -1) THEN 1366 GOTO 100 1367 ELSE 1368 LRI = LRI + 3 1369 GOTO 100 1370 END IF 1371 END IF 1372C 1373 1374C Unpack the integrals 1375 CALL DZERO(WORK(KEPSx),3*N2BST(ISYMOP)) 1376 LABEL = 'GIVE INT' 1377 CALL CCMMINT(LABEL,LUQMMM,WORK(KEPSx),WORK(KEPSAO), 1378 * IRREP,ISYM,IERR,WORK(KWRK),LWRK) 1379 CALL CCMMINT(LABEL,LUQMMM,WORK(KEPSy), 1380 * WORK(KEPSAO+NNBASX),IRREP,ISYM,IERR, 1381 * WORK(KWRK),LWRK) 1382 CALL CCMMINT(LABEL,LUQMMM,WORK(KEPSz), 1383 * WORK(KEPSAO+2*NNBASX),IRREP,ISYM,IERR, 1384 * WORK(KWRK),LWRK) 1385C Dot with the AO density 1386 1387 EXELCO = DDOT(N2BST(ISYMOP),AODEN,1,WORK(KEPSx),1) 1388 EYELCO = DDOT(N2BST(ISYMOP),AODEN,1,WORK(KEPSy),1) 1389 EZELCO = DDOT(N2BST(ISYMOP),AODEN,1,WORK(KEPSz),1) 1390 1391 IF (LOCDEB) THEN 1392 WRITE(LUPRI,*) 'At MM site ', I, ' the field due to the QM', 1393 & ' electrons: ',EXELCO,EYELCO,EZELCO 1394 END IF 1395 1396C Add to the electric field from the multipoles 1397 ELF(LRI+0) = ELF(LRI+0) + EXELCO 1398 ELF(LRI+1) = ELF(LRI+1) + EYELCO 1399 ELF(LRI+2) = ELF(LRI+2) + EZELCO 1400 1401C Second we add the field due to the distributed multipoles 1402 CALL CCMM_FMUL(ELF,LRI,I) 1403 1404C Finally add the field from the QM nuclei 1405 CALL CCMM_FNUC(ELF,LRI,I) 1406 1407 LRI = LRI + 3 1408 1409 100 CONTINUE !Sum over all sites - done with the static field 1410 1411 CALL QEXIT('CCMM_FSTATIC') 1412 END 1413C************************************************************************************** 1414C /* Deck put_to_file */ 1415 SUBROUTINE PUT_TO_FILE(FLNAME,NULOOP,DDATA) 1416C************************************************************** 1417C 1418#include "implicit.h" 1419#include "dummy.h" 1420C 1421 CHARACTER*(*) FLNAME 1422 INTEGER NMBU,NULOOP 1423 DIMENSION DDATA(*) 1424C 1425 NMBU = -1 1426 CALL GPOPEN(NMBU,FLNAME,'UNKNOWN',' ','FORMATTED',IDUMMY,.FALSE.) 1427C 1428 REWIND (NMBU) 1429 DO 820 L = 1,NULOOP 1430 WRITE(NMBU,'(I5,3E25.15)') L,DDATA(L) 1431 820 CONTINUE 1432C 1433 CALL GPCLOSE(NMBU,'KEEP') 1434C 1435 END 1436C 1437C************************************************************** 1438C /* Deck get_from_file */ 1439 SUBROUTINE GET_FROM_FILE(FLNAME,NULOOP,DDATA) 1440C************************************************************** 1441C 1442#include "implicit.h" 1443#include "dummy.h" 1444C 1445 CHARACTER*(*) FLNAME 1446 INTEGER NMBU,NULOOP 1447 DIMENSION DDATA(*) 1448C 1449 NMBU = -1 1450 CALL GPOPEN(NMBU,FLNAME,'UNKNOWN',' ','FORMATTED',IDUMMY,.FALSE.) 1451C 1452 REWIND (NMBU) 1453 DO 820 L = 1,NULOOP 1454 READ(NMBU,'(I5,3E25.15)') LK,DDATA(L) 1455 820 CONTINUE 1456C 1457 IF (LK.NE.NULOOP) THEN 1458 CALL QUIT('Problem in dimension in GET_FROM_FILE') 1459 ENDIF 1460 1461 CALL GPCLOSE(NMBU,'KEEP') 1462C 1463 END 1464C 1465C*********************************************************************** 1466C /* Deck ccmm_epsao */ 1467 SUBROUTINE CCMM_EPSAO(EPSAO,SITE,LSKIP,WORK,LWORK) 1468C 1469C July 2010, KS 1470C 1471C Purpose: Calculate the electric field at a site SITE due to 1472C the electrons. LSKIP keeps track of whether a site is skipped 1473C either due to having no polarizability or because of distance 1474C cutoff. The vector is returned in the AO basis. 1475C Note we are doing everything integral-direct to avoid 1476C a possible storage-problem when working with large systems. 1477 1478#include "implicit.h" 1479#include "maxorb.h" 1480#include "mxcent.h" 1481#include "dummy.h" 1482#include "iratdef.h" 1483#include "priunit.h" 1484#include "ccorb.h" 1485#include "ccsdsym.h" 1486#include "ccsdio.h" 1487#include "ccsdinp.h" 1488#include "ccfield.h" 1489#include "exeinf.h" 1490#include "ccfdgeo.h" 1491#include "ccslvinf.h" 1492#include "ccinftap.h" 1493#include "nuclei.h" 1494#include "orgcom.h" 1495#include "qmmm.h" 1496#include "qm3.h" 1497#include "cclr.h" 1498 1499 PARAMETER (D3I = 1.0D0/3.0D0, D6I = 1.0D0/6.0D0) 1500 DIMENSION WORK(LWORK) 1501 DIMENSION EPSAO(*) 1502 INTEGER SITE 1503 LOGICAL LSKIP 1504 LOGICAL TOFILE,TRIMAT,EXP1VL,LOCDEB 1505 PARAMETER (LOCDEB=.FALSE.) 1506 DIMENSION INTREP(9*MXCENT), INTADR(9*MXCENT) 1507 CHARACTER*8 LABINT(9*MXCENT) 1508 1509 KWRK = 1 1510 LWRK = LWORK - KWRK 1511 1512 IF (ZEROAL(SITE) .EQ. -1) THEN 1513 LSKIP = .TRUE. 1514 RETURN 1515 END IF 1516 1517C Backup dipole origin 1518 OBKPX = DIPORG(1) 1519 OBKPY = DIPORG(2) 1520 OBKPZ = DIPORG(3) 1521 1522 DIPORG(1) = MMCORD(1,SITE) 1523 DIPORG(2) = MMCORD(2,SITE) 1524 DIPORG(3) = MMCORD(3,SITE) 1525 1526 DIST2 = (MMCORD(1,SITE)-QMCOM(1))**2 + 1527 & (MMCORD(2,SITE)-QMCOM(2))**2 + 1528 & (MMCORD(3,SITE)-QMCOM(3))**2 1529 DIST = SQRT(DIST2) 1530 1531 IF (DIST .GT. RCUTMM) THEN 1532 LSKIP =.TRUE. 1533 IF (LOCDEB) THEN 1534 WRITE(LUPRI,*) 'Skipping site in mu^ind vector 1535 & in CCMM_POL', SITE 1536 ENDIF 1537 RETURN 1538 ENDIF 1539 1540 CALL DZERO(EPSAO,3*NNBASX) 1541 KPATOM = 0 1542 NOCOMP = 3 1543 TOFILE = .FALSE. 1544 TRIMAT = .TRUE. 1545 EXP1VL = .FALSE. 1546 1547 RUNQM3 = .TRUE. 1548 CALL GET1IN(EPSAO,'NEFIELD',NOCOMP,WORK(KWRK), 1549 * LWRK,LABINT,INTREP,INTADR,SITE,TOFILE,KPATOM, 1550 * TRIMAT,DUMMY,EXP1VL,DUMMY,IPQMMM) 1551 RUNQM3 = .FALSE. 1552 1553 IF (QMDAMP) THEN 1554 IF ( (IDAMP .EQ. 3) .AND. (NQMNUC .NE. NUCIND) ) THEN 1555 CALL QUIT('ERROR in no. of assigned QM 1556 & polarizabilities') 1557 ENDIF 1558 IF ( (IDAMP .EQ. 1) .OR. (IDAMP .EQ. 3) ) THEN 1559 DIST = 9.99D+99 1560 MHIT = 0 1561 DO 205 M=1,NUCIND 1562 DISTC = (DIPORG(1)-CORD(1,M))**2 + 1563 & (DIPORG(2)-CORD(2,M))**2 + 1564 & (DIPORG(3)-CORD(3,M))**2 1565 IF (DISTC .LE. DIST) THEN 1566 DIST = DISTC 1567 MHIT = M 1568 ENDIF 1569 205 CONTINUE 1570 ELSE IF (IDAMP .EQ. 2) THEN 1571 DIST = (DIPORG(1)-QMCOM(1))**2 + 1572 & (DIPORG(2)-QMCOM(2))**2 + 1573 & (DIPORG(3)-QMCOM(3))**2 1574 ENDIF 1575 DIST = SQRT(DIST) 1576C 1577 IF (IDAMP .EQ. 3) THEN 1578 IF (IPOLTP .EQ. 2) THEN 1579 TEMPI = D3I*(POLMM(1,SITE)+POLMM(4,SITE)+POLMM(6,SITE)) 1580 ELSE IF (IPOLTP .EQ. 1) THEN 1581 TEMPI = POLIMM(SITE) 1582 ENDIF 1583 TEMP = (TEMPI*QMPOL(MHIT))**(D6I) 1584 SIJ = 2.1304D0*DIST/TEMP 1585 DFACT = 1.0D0 - (1.0D0+SIJ+0.50D0*SIJ*SIJ)*exp(-SIJ) 1586 ELSE 1587 DFACT = (1-exp(-ADAMP*DIST))**3 1588 ENDIF 1589 CALL DSCAL(3*NNBASX,DFACT,EPSAO,1) 1590 ENDIF 1591 1592C Print section 1593 IF ( (IPQMMM.GT.15) .OR. (LOCDEB) ) THEN 1594 WRITE(LUPRI,*) 'Packed AO electric field_x matrix due to QM 1595 & electrons at MM site ',SITE 1596 CALL OUTPAK(EPSAO(1),NBAST,1,LUPRI) 1597 1598 WRITE(LUPRI,*) 'Packed AO electric field_y matrix due to QM 1599 & electrons at MM site ',SITE 1600 CALL OUTPAK(EPSAO(1+NNBASX),NBAST,1,LUPRI) 1601 1602 WRITE(LUPRI,*) 'Packed AO electric field_z matrix due to QM 1603 & electrons at MM site ',SITE 1604 CALL OUTPAK(EPSAO(1+2*NNBASX),NBAST,1,LUPRI) 1605 ENDIF 1606C Restore dipole origin 1607 DIPORG(1) = OBKPX 1608 DIPORG(2) = OBKPY 1609 DIPORG(3) = OBKPZ 1610 1611 END 1612C************************************************************** 1613C /* Deck cc_qmmme */ 1614 SUBROUTINE CC_QMMME(ECCMM,EPOL,EEL,EELMUL,EREP,AODEN,WORK,LWORK) 1615C************************************************************** 1616C 1617C KS, 09 1618C Purpose: Updates the induced dipoles besides calculating E(QM/MM) 1619C 1620C In new QMMM implementation (NYQMMM) the strategy is as follows 1621C 1) Calculate F^sta=F^el+F^mul+F^nuc 1622C 2) Calculate induced dipole by transforming with relay 1623C matrix (MMMAT) or by iterative means (MMITER) 1624C 3) Put induced moments to file 1625C 4) Calculate Epol=-1/2sum_a\mu^ind_aF^sta_a 1626C 5) Calculate Eelesta=-sum_s<N_s> 1627C 1628C Also the LGFILE keyword is reset such that a new G matrix is 1629C written to file when TGT is called next time. 1630C 1631C************************************************************** 1632C 1633#include "implicit.h" 1634#include "priunit.h" 1635#include "dummy.h" 1636#include "mxcent.h" 1637#include "qmmm.h" 1638#include "qm3.h" 1639#include "iratdef.h" 1640#include "maxorb.h" 1641#include "orgcom.h" 1642#include "nuclei.h" 1643#include "ccinftap.h" 1644#include "ccorb.h" 1645#include "ccsdsym.h" 1646#include "ccsdio.h" 1647#include "ccsdinp.h" 1648#include "ccfield.h" 1649#include "exeinf.h" 1650#include "ccfdgeo.h" 1651#include "thrzer.h" 1652 1653C 1654 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, IZERO = 0 , TWO = 2.0D0) 1655 PARAMETER (HALF = 1.0D0/2.0D0) 1656 PARAMETER (D3I = 1.0D0/3.0D0, D6I = 1.0D0/6.0D0) 1657C 1658 DIMENSION INTREP(9*MXCENT), INTADR(9*MXCENT) 1659 CHARACTER*8 LABINT(9*MXCENT) 1660 LOGICAL TOFILE, TRIMAT, EXP1VL, LOCDEB, FNDLAB, EXCENT 1661 PARAMETER (LOCDEB=.FALSE.) 1662 DIMENSION WORK(LWORK),AODEN(*) 1663 CHARACTER LABEL*8 1664 INTEGER NTOTI 1665C 1666 CALL QENTER('CC_QMMME') 1667C------------------------------------------------------------------ 1668C Init parameters 1669C------------------------------------------------------------------ 1670 1671 ISYMOP = 1 1672C 1673C------------------------------------------------------------------ 1674C Dynamical allocation for CCMM 1675C------------------------------------------------------------------ 1676C 1677 NDIM = 3*NNZAL 1678 LGFILE=.FALSE. ! reset G matrix construction 1679 KELF = 1 1680 KINDIP = KELF + NDIM 1681 KNSAO = KINDIP + NDIM 1682 KINT = KNSAO + NNBASX 1683 KWRK1 = KINT + N2BST(ISYMOP) 1684 LWRK1 = LWORK - KWRK1 1685C 1686 CALL DZERO(WORK(KELF),NDIM) 1687 CALL DZERO(WORK(KINDIP),NDIM) 1688C 1689C 1690 IF (LWRK1.LT.0) CALL QUIT( 'Too little work in CC_QMMME') 1691C 1692C SPLDIP option not implemented for CC 1693 IF (SPLDIP) WRITE(LUPRI,*) 'WARNING: SPLDIP option 1694 & not implemented for CC/MM. Keyword ignored!' 1695 IF (FIXDIP) WRITE(LUPRI,*) 'WARNING: FIXDIP option 1696 & not implemented for CC/MM. Keyword ignored!' 1697 IF (LGSPOL) CALL QUIT( 'LGSPOL not implemented for CC/MM') 1698 1699C 1700 IF ( ( IPOLTP .GT. 0 ) .AND. (.NOT. HFFLD) ) THEN 1701 CALL CCMM_FSTATIC(WORK(KELF),AODEN,WORK(KWRK1),LWRK1) 1702C 1703CC Transform the field to induced dipoles 1704 CALL CCMM_F2DIP(WORK(KELF),WORK(KINDIP),WORK(KWRK1),LWRK1) 1705C 1706CC Write the updated induced moments to file 1707 CALL PUT_TO_FILE_1('INDUCED_DIPOLES',NNZAL,WORK(KINDIP)) 1708 ELSE IF ( HFFLD .AND. (IPOLTP .GT. 0 ) ) THEN 1709 CALL CCMM_FSTATIC(WORK(KELF),AODEN,WORK(KWRK1),LWRK1) 1710 CALL GET_FROM_FILE_1('INDUCED_DIPOLES',NNZAL,WORK(KINDIP)) 1711CC 1712 END IF 1713 1714 ECCMM = ZERO 1715 EEL = ZERO 1716 EELMUL = ZERO 1717 EPOL = ZERO 1718 ETEMP = ZERO 1719C------------------------------------- 1720C Add up the energy contributions: 1721C 1) Epol: 1722C------------------------------------- 1723 1724 ETEMP = DDOT(NDIM,WORK(KINDIP),1,WORK(KELF),1) 1725 1726 EPOL = -HALF*ETEMP 1727 1728C Now add the electrostatic energy contribution 1729C we do this by getting the multipole integrals from file 1730C and taking the expectation value. The integrals are at our disposal 1731C from the preceeding HF/MM calculation handled by sirius 1732C-------- 1733C 2) The electrostatic contribution: 1734C-------- 1735C 1736 CALL DZERO(WORK(KNSAO),NNBASX) 1737 CALL CCMM_ELSTAT(WORK(KNSAO),WORK(KWRK1),LWRK1) 1738 1739C Unpack the multipole integrals and dot with the density to 1740C obtain the expectation value 1741 1742 LABEL = 'GIVE INT' 1743 CALL CCMMINT(LABEL,LUQMMM,WORK(KINT),WORK(KNSAO), 1744 * IRREP,ISYM,IERR,WORK(KWRK1),LWRK1) 1745 1746 EELMUL = DDOT(N2BST(ISYMOP),WORK(KINT),1,AODEN,1) 1747 1748C Finally we add the interaction energy between QM nuclei and MM multipoles. This 1749C was found in the previous SCF/MM calculation and is readily available 1750C in the COMMON block. 1751 EEL = EELMUL + ENUMUL 1752 1753 ECCMM = EPOL + EEL 1754 IF (LOCDEB .OR. IPQMMM .GT. 4) THEN 1755 WRITE(LUPRI,*) 'EELMUL: ', EELMUL 1756 WRITE(LUPRI,*) 'EELSTAT: ', EEL 1757 WRITE(LUPRI,*) 'EPOL: ', EPOL 1758 WRITE(LUPRI,*) 'ECCMM: ', ECCMM 1759 WRITE(LUPRI,*) 'The nuclear-multipole contribution: ', ENUMUL 1760 END IF 1761C 1762C 1763C------------- 1764C 3) Edisp 1765C------------- 1766C 1767C sne - New qmmm not implemented for dispersion yet 1768C sne - Anyhow it is just a static energy contribution (for now) so 1769C straightforward to do - copy how it is done in dft 1770C ECCMM = ECCMM + ECLVDW 1771C 1772C----------------------------------------------- 1773C sne New qmmm not implemented for explicit repulsion yet 1774C 4) E_repulsion 1775C 1776C CALL CCMM_REP2(EREP,AODEN,WORK(KWRK1),LWRK1) 1777C ECCMM = ECCMM + EREP 1778C 1779C----------------------------------------------- 1780C Writing out the final energy contributions one by one: 1781C----------------------------------------------- 1782C 1783 1050 FORMAT(' Induced dipole moments ') 1784 1051 FORMAT(2X,'=',22('-'),'=',2X) 1785 1000 FORMAT(1X,51('=')) 1786 1020 FORMAT(1X,I6,3(4X,F10.6)) 1787 1010 FORMAT(' | Site | X | Y | Z |') 1788 1030 FORMAT(' Total induced dipole moments: ') 1789 1790 CALL QEXIT('CC_QMMME') 1791C 1792 END 1793C 1794C*********************************************************************** 1795C /* Deck ccmm_fmul */ 1796 SUBROUTINE CCMM_FMUL(ELF,LRI,SITEI) 1797C 1798C Aug 2010, KS 1799C 1800C Purpose: Calculate the electric field due to the MM multipoles 1801C at SITEI 1802C 1803 1804#include "implicit.h" 1805#include "maxorb.h" 1806#include "mxcent.h" 1807#include "dummy.h" 1808#include "iratdef.h" 1809#include "priunit.h" 1810#include "ccorb.h" 1811#include "ccsdsym.h" 1812#include "ccsdio.h" 1813#include "ccsdinp.h" 1814#include "ccfield.h" 1815#include "exeinf.h" 1816#include "ccfdgeo.h" 1817#include "ccslvinf.h" 1818#include "ccinftap.h" 1819#include "nuclei.h" 1820#include "orgcom.h" 1821#include "qm3.h" 1822#include "qmmm.h" 1823#include "cclr.h" 1824 1825 DIMENSION ELF(*) 1826 INTEGER SITEI, SITEJ, LRI 1827 LOGICAL LOCDEB, EXCENT 1828 PARAMETER (LOCDEB=.FALSE.) 1829 1830 CALL QENTER('CCMM_FMUL') 1831 1832 DO 200 SITEJ=1,MMCENT 1833 1834 IF (SITEJ .NE. SITEI) THEN 1835 DIST2 = (MMCORD(1,SITEI)-MMCORD(1,SITEJ))**2 + 1836 & (MMCORD(2,SITEI)-MMCORD(2,SITEJ))**2 + 1837 & (MMCORD(3,SITEI)-MMCORD(3,SITEJ))**2 1838 DIST = SQRT(DIST2) 1839 1840 IF(DIST .GT. RCUTMM) THEN 1841 IF (LOCDEB) THEN 1842 WRITE(LUPRI,*) 1843 & 'Skipping element in T2 CC', SITEI,SITEJ 1844 ENDIF 1845 GOTO 200 1846 ENDIF 1847 1848 EXCENT=.FALSE. 1849 1850C Check if any centers were in the same LoProp calculation 1851 IF (NEWEXC) THEN 1852 DO L = 1, NEXLST 1853 IF (EXLIST(1,SITEI) .EQ. EXLIST(L,SITEJ)) EXCENT = .TRUE. 1854 ENDDO 1855 ELSE 1856 DO L = 1, NEXLST 1857 IF (EXLIST(L,SITEI) .EQ. EXLIST(1,SITEJ)) EXCENT = .TRUE. 1858 ENDDO 1859 ENDIF 1860 1861C Form F vector due to the permament moments 1862 1863 IF (.NOT. (EXCENT) ) THEN 1864 1865C A) The point-charge contribution 1866 IF((ABS(MUL0MM(SITEJ)) .GT. THRMM) 1867 & .AND. (NMULT.GE.0) ) THEN 1868 1869 ELFLDX = 0.0D0 1870 ELFLDY = 0.0D0 1871 ELFLDZ = 0.0D0 1872 1873 CALL GET_CHARGE_ELFLD(MUL0MM(SITEJ), 1874 * MMCORD(1,SITEJ),MMCORD(2,SITEJ),MMCORD(3,SITEJ), 1875 * MMCORD(1,SITEI),MMCORD(2,SITEI),MMCORD(3,SITEI), 1876 * ELFLDX,ELFLDY,ELFLDZ) 1877 ELF(LRI+0) = ELF(LRI+0)+ELFLDX 1878 ELF(LRI+1) = ELF(LRI+1)+ELFLDY 1879 ELF(LRI+2) = ELF(LRI+2)+ELFLDZ 1880 1881 ENDIF 1882 1883C B) The dipole contribution 1884 IF (NMULT .GE. 1) THEN 1885 ELFLDX = 0.0D0 1886 ELFLDY = 0.0D0 1887 ELFLDZ = 0.0D0 1888 CALL GET_DIPOLE_ELFLD(MUL1MM(1,SITEJ), 1889 * MUL1MM(2,SITEJ),MUL1MM(3,SITEJ), 1890 * MMCORD(1,SITEJ),MMCORD(2,SITEJ),MMCORD(3,SITEJ), 1891 * MMCORD(1,SITEI),MMCORD(2,SITEI),MMCORD(3,SITEI), 1892 * ELFLDX,ELFLDY,ELFLDZ) 1893 ELF(LRI+0) = ELF(LRI+0)+ELFLDX 1894 ELF(LRI+1) = ELF(LRI+1)+ELFLDY 1895 ELF(LRI+2) = ELF(LRI+2)+ELFLDZ 1896 ENDIF 1897 1898C C) The quadrupole contribution 1899 IF (NMULT .GE. 2) THEN 1900 ELFLDX = 0.0D0 1901 ELFLDY = 0.0D0 1902 ELFLDZ = 0.0D0 1903 CALL GET_QUADRUPOLE_ELFLD( 1904 * MUL2MM(1,SITEJ),MUL2MM(2,SITEJ),MUL2MM(3,SITEJ), 1905 * MUL2MM(4,SITEJ),MUL2MM(5,SITEJ),MUL2MM(6,SITEJ), 1906 * MMCORD(1,SITEJ),MMCORD(2,SITEJ),MMCORD(3,SITEJ), 1907 * MMCORD(1,SITEI),MMCORD(2,SITEI),MMCORD(3,SITEI), 1908 * ELFLDX,ELFLDY,ELFLDZ) 1909 ELF(LRI+0) = ELF(LRI+0)+ELFLDX 1910 ELF(LRI+1) = ELF(LRI+1)+ELFLDY 1911 ELF(LRI+2) = ELF(LRI+2)+ELFLDZ 1912 ENDIF 1913 ENDIF 1914 END IF 1915 200 CONTINUE 1916 1917 CALL QEXIT('CCMM_FMUL') 1918 END 1919C*********************************************************************** 1920C /* Deck ccmm_fnuc */ 1921 SUBROUTINE CCMM_FNUC(ELF,LRI,SITEI) 1922C 1923C Aug 2010, KS 1924C 1925C Purpose: Calculate the electric field due to the QM nuclei 1926C at SITEI 1927C 1928 1929#include "implicit.h" 1930#include "maxorb.h" 1931#include "mxcent.h" 1932#include "dummy.h" 1933#include "iratdef.h" 1934#include "priunit.h" 1935#include "ccorb.h" 1936#include "ccsdsym.h" 1937#include "ccsdio.h" 1938#include "ccsdinp.h" 1939#include "ccfield.h" 1940#include "exeinf.h" 1941#include "ccfdgeo.h" 1942#include "ccslvinf.h" 1943#include "ccinftap.h" 1944#include "nuclei.h" 1945#include "orgcom.h" 1946#include "qm3.h" 1947#include "qmmm.h" 1948#include "cclr.h" 1949 1950 DIMENSION ELF(*) 1951 INTEGER SITEI, LRI, SITEJ 1952 LOGICAL LOCDEB 1953 PARAMETER (LOCDEB=.FALSE.) 1954 PARAMETER (D3I = 1.0D0/3.0D0, D6I = 1.0D0/6.0D0) 1955 1956 TJKX = 0.0D0 1957 TJKY = 0.0D0 1958 TJKZ = 0.0D0 1959 1960 DO 145 SITEJ=1,NUCIND ! Dalton keyword for nr of QM nuclei 1961 ELFLDX = 0.0D0 1962 ELFLDY = 0.0D0 1963 ELFLDZ = 0.0D0 1964 CALL GET_CHARGE_ELFLD(CHARGE(SITEJ), 1965 * CORD(1,SITEJ),CORD(2,SITEJ),CORD(3,SITEJ), 1966 * MMCORD(1,SITEI),MMCORD(2,SITEI),MMCORD(3,SITEI), 1967 * ELFLDX,ELFLDY,ELFLDZ) 1968 1969 IF (QMDAMP) THEN 1970 IF ((IDAMP .EQ. 3).AND.(NQMNUC .NE. NUCIND)) THEN 1971 CALL QUIT('ERROR in no. of assigned QM 1972 & polarizabilities') 1973 ENDIF 1974 IF ( (IDAMP .EQ. 1) .OR. (IDAMP .EQ. 3) ) THEN 1975 DIQM = 9.99D+99 1976 MHIT = 0 1977 DO 155 M=1,NUCIND 1978 DIQMC = (MMCORD(1,SITEI)-CORD(1,M))**2 + 1979 & (MMCORD(2,SITEI)-CORD(2,M))**2 + 1980 & (MMCORD(3,SITEI)-CORD(3,M))**2 1981 IF (DIQMC .LE. DIQM) THEN 1982 DIQM = DIQMC 1983 MHIT = M 1984 ENDIF 1985 155 CONTINUE 1986 ELSE IF (IDAMP .EQ. 2) THEN 1987 DIQM = (MMCORD(1,SITEI)-QMCOM(1))**2 + 1988 & (MMCORD(2,SITEI)-QMCOM(2))**2 + 1989 & (MMCORD(3,SITEI)-QMCOM(3))**2 1990 ENDIF 1991 DIQM = SQRT(DIQM) 1992 1993 IF (IDAMP .EQ. 3) THEN 1994 IF (IPOLTP .EQ. 2) THEN 1995 TEMPI =D3I*(POLMM(1,SITEI)+POLMM(4,SITEI)+POLMM(6,SITEI)) 1996 ELSE IF (IPOLTP .EQ. 1) THEN 1997 TEMPI = POLIMM(SITEI) 1998 ENDIF 1999 TEMP = (TEMPI*QMPOL(MHIT))**(D6I) 2000 SIJ = 2.1304D0*DIQM/TEMP 2001 DFACT=1.0D0-(1.0D0+SIJ+0.50D0*SIJ*SIJ)*exp(-SIJ) 2002 ELSE 2003 DFACT = (1-exp(-ADAMP*DIQM))**3 2004 ENDIF 2005 2006 ELFLDX = ELFLDX*DFACT 2007 ELFLDY = ELFLDY*DFACT 2008 ELFLDZ = ELFLDZ*DFACT 2009 END IF 2010 2011 TJKX = TJKX + ELFLDX 2012 TJKY = TJKY + ELFLDY 2013 TJKZ = TJKZ + ELFLDZ 2014 145 CONTINUE 2015 ELF(LRI+0)=ELF(LRI+0)+TJKX 2016 ELF(LRI+1)=ELF(LRI+1)+TJKY 2017 ELF(LRI+2)=ELF(LRI+2)+TJKZ 2018 2019 2020 IF (LOCDEB) THEN 2021 WRITE(LUPRI,*) 'Nuclear field: CC ', TJKX,TJKY,TJKZ 2022 ENDIF 2023 2024 END 2025C*********************************************************************** 2026C 2027C /* Deck ccmm_qrtransformer */ 2028 SUBROUTINE CCMM_QRTRANSFORMER(RHO1,RHO2,ISYRES, 2029 * LISTB,IDLSTB,ISYMTB, 2030 * LISTC,IDLSTC,ISYMTC, 2031 * MODEL,RSPTYP,WORK,LWORK) 2032C 2033C----------------------------------------------------------------------------- 2034C 2035C Construct effective operators from contracted densities and do 2036C response transformation. This is the QR analog to the TINE 2037C CCMM_LRTRANSFORMER. 2038C For B and G we have 3 contributions 2039C For K transformation 2 contributions 2040 2041C A: Lagrangian multipliers to be used for F transformation 2042C B: B trial vector 2043C C: C trial vector 2044C 2045C The strategy is as follows: 2046C 1) Construct auxiliary density 2047C 2) Transform density to electric field 2048C 3) Transform electric field to induced dipoles 2049C 4) Construct effective operator. G^C=-sum_a\∆mu^C_a eps_a 2050C 5) Do response-transformation the type of which determined by QR flag 2051C 6) Add to rho vector 2052C 2053C KS, aug 10 2054C KS, Sep 10: Modified to be called from FMATNEW to generate TPA pseudo 2055C B matrix contributions 2056C----------------------------------------------------------------------------- 2057C 2058#include "implicit.h" 2059#include "maxorb.h" 2060#include "mxcent.h" 2061#include "qmmm.h" 2062#include "dummy.h" 2063#include "iratdef.h" 2064#include "priunit.h" 2065#include "ccorb.h" 2066#include "ccsdsym.h" 2067#include "ccsdio.h" 2068#include "ccsdinp.h" 2069#include "ccfield.h" 2070#include "exeinf.h" 2071#include "ccfdgeo.h" 2072#include "ccslvinf.h" 2073#include "qm3.h" 2074#include "ccinftap.h" 2075#include "nuclei.h" 2076#include "orgcom.h" 2077 2078 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, IZERO = 0 , TWO = 2.0D0) 2079 PARAMETER (HALF = 0.5D0) 2080C 2081 DIMENSION WORK(LWORK),RHO1(*),RHO2(*) 2082C 2083 CHARACTER*(*) LISTB, LISTC 2084 INTEGER IDLSTB, IDLSTC 2085 INTEGER ISYCB, ISYRES, ISYMTB, ISYMTC 2086 CHARACTER MODEL*10 2087 LOGICAL LOCDEB, LSAME 2088 SAVE LOCDEB 2089 DATA LOCDEB /.FALSE./ 2090 CHARACTER*(2) LISTL 2091C 2092 CHARACTER*8 LABEL,DENTYP*(1),RSPTYP*(1) 2093C 2094C 2095 CALL QENTER('CCMM_QRTRANSFORMER') 2096 LISTL='L0' 2097 2098 IF (CCS) CALL QUIT('CCMM_QRTRANS: Not implemented for CCS') 2099C 2100 IF (DISCEX) CALL QUIT('CCMM_QRTRANS:Not implemented for DISCEX') 2101 2102 IF (IPQMMM.GT.10 .OR. LOCDEB) THEN 2103 WRITE(LUPRI,*) 'CCMM_QRTRANSFORMER: RSPTYP: ',RSPTYP 2104 WRITE(LUPRI,*)'CCMM_QRTRANSFORMER: ISYRES= ',ISYRES 2105 WRITE(LUPRI,*)'CCMM_QRTRANSFORMER: ISYMTB= ',ISYMTB 2106 WRITE(LUPRI,*)'CCMM_QRTRANSFORMER: ISYMTC= ',ISYMTC 2107 WRITE(LUPRI,*)'CCMM_QRTRANSFORMER: LISTB= ',LISTB 2108 WRITE(LUPRI,*)'CCMM_QRTRANSFORMER: LISTC= ',LISTC 2109 WRITE(LUPRI,*)'CCMM_QRTRANSFORMER: IDLSTB= ',IDLSTB 2110 WRITE(LUPRI,*)'CCMM_QRTRANSFORMER: IDLSTC= ',IDLSTC 2111 CALL FLSHFO(LUPRI) 2112 ENDIF 2113 2114 IF (IPOLTP .EQ. 0) THEN 2115 WRITE(LUPRI,*) 'No QM/MM polarization so no contribution 2116 & from quadratic response transformer.' 2117 CALL QEXIT('CCMM_QRTRANSFORMER') 2118 RETURN 2119 END IF 2120C 2121C Check the symmetry 2122 ISYCB=MULD2H(ISYMTB,ISYMTC) 2123 IF (ISYCB .NE. ISYRES) THEN 2124 CALL QUIT('Symmetry problem in CCMM_QRTRANSFORMER') 2125 END IF 2126C 2127C Check what kind of response. 2128C Find out which density and hence which operator to construct 2129C For F start with left density and change flag when done with this cont. 2130 IF ( (RSPTYP.EQ.'G') .OR. (RSPTYP.EQ.'B') ) THEN 2131 DENTYP='R' 2132 ELSE IF (RSPTYP.EQ.'K' .OR. (RSPTYP .EQ. 'F')) THEN 2133 DENTYP='L' 2134 END IF 2135C 2136C--------------------- 2137C Init parameters. 2138C--------------------- 2139C 2140 NB1 = NT1AM(ISYMTB) 2141 NB2 = NT2AM(ISYMTB) 2142 NDIM = 3*NNZAL 2143C 2144C------------------------------------ 2145C Dynamical allocation for CCMM : 2146C------------------------------------ 2147 KB1=1 2148 KB2=KB1+NB1 2149 KDENS=KB2+NB2 2150 KGBMAT=KDENS+N2BST(ISYMTB) 2151 KWRK1=KGBMAT+N2BST(ISYMTB) 2152 LWRK1=LWORK-KWRK1 2153 2154 IF (LWRK1 .LT. 0) THEN 2155 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KWRK1 2156 CALL QUIT( 'Too little work in CCMM_QRTRANSFORMER, 1') 2157 END IF 2158C 2159C Initialize the working matrices 2160 CALL DZERO(WORK(KB1),NB1) 2161 CALL DZERO(WORK(KB2),NB2) 2162 CALL DZERO(WORK(KDENS),N2BST(ISYMTB)) 2163 CALL DZERO(WORK(KGBMAT),N2BST(ISYMTB)) 2164C 2165C Read C vector from file 2166 IOPT = 3 2167 CALL CC_RDRSP(LISTB,IDLSTB,ISYMTB,IOPT,MODEL, 2168 * WORK(KB1),WORK(KB2)) 2169C 2170C Print section 2171 IF ( IPQMMM .GT. 10 .OR. LOCDEB) THEN 2172 RHO1N = DDOT(NT1AM(ISYCB),RHO1,1,RHO1,1) 2173 RHO2N = DDOT(NT2AM(ISYCB),RHO2,1,RHO2,1) 2174 WRITE(LUPRI,*)'Norm of RHO1 in CCMM_QRTRANS on input:,1', RHO1N 2175 WRITE(LUPRI,*)'Norm of RHO2 in CCMM_QRTRANS on input:,1', RHO2N 2176 RHO1N = DDOT(NB1,WORK(KB1),1,WORK(KB1),1) 2177 RHO2N = DDOT(NB2,WORK(KB2),1,WORK(KB2),1) 2178 WRITE(LUPRI,*)'Norm af B1AM in CCMM_QRTRANS on input:,1', RHO1N 2179 WRITE(LUPRI,*)'Norm af B2AM in CCMM_QRTRANS on input:,1', RHO2N 2180 ENDIF 2181C---------------------------- 2182C Construct auxiliary densities 2183C---------------------------- 2184 CALL CCMM_D1AO(WORK(KDENS),WORK(KB1),WORK(KB2),MODEL,DENTYP, 2185 * LISTC,IDLSTC,ISYMTC, 2186 * WORK(KWRK1),LWRK1) 2187 2188C---------------------------- 2189C Construct effective operator 2190C---------------------------- 2191 CALL CCMM_GEFF(WORK(KDENS),WORK(KGBMAT),WORK(KWRK1),LWRK1) 2192C 2193C Transform the effective operator according to RESPTYP flag 2194 IF(LOCDEB .OR. IPQMMM .GT. 14) THEN 2195 TAL1= DDOT(N2BST(ISYCB),WORK(KGBMAT),1,WORK(KGBMAT),1) 2196 WRITE(LUPRI,*) 'Print Norm2 GBMAT: ', TAL1 2197 END IF 2198C 2199C Dynamical memory allocation to store the response transformed operator 2200 KETA1=KWRK1 2201 KETA2=KETA1+NT1AM(ISYCB) 2202 KWRK2=KETA2+NT2AM(ISYCB) 2203 LWRK2=LWORK-KWRK2 2204 2205 LABEL='GIVE INT' 2206 IF ( (RSPTYP.EQ.'G') .OR. (RSPTYP.EQ.'F') ) THEN 2207 CALL CCLR_FA(LABEL,ISYMTB,LISTC,IDLSTC,LISTL,0, 2208 * WORK(KGBMAT),WORK(KETA1),LWRK2) 2209 ELSE IF (RSPTYP.EQ.'B') THEN 2210 CALL CCCR_AA(LABEL,ISYMTB,LISTC,IDLSTC,WORK(KGBMAT), 2211 * WORK(KETA1),LWRK2) 2212 CALL CCLR_DIASCL(WORK(KETA2),TWO,ISYCB) 2213 ELSE IF (RSPTYP.EQ.'K') THEN 2214 CALL CC_ETAC(ISYMTB,LABEL,WORK(KETA1),LISTC,IDLSTC,0, 2215 * WORK(KGBMAT),WORK(KWRK2),LWRK2) 2216 END IF 2217 2218 2219C Print section 2220 IF (LOCDEB .OR. IPQMMM .GT. 14) THEN 2221 TAL1= DDOT(NT1AM(ISYCB),WORK(KETA1),1,WORK(KETA1),1) 2222 TAL2= DDOT(NT2AM(ISYCB),WORK(KETA2),1,WORK(KETA2),1) 2223 WRITE(LUPRI,*) 'Printing transformed G^B*C contribution. 2224 & Norm2 of singles: ', TAL1, 2225 & 'Norm2 of doubles: ', TAL2 2226 END IF 2227 2228C Note that if B=C then we just multiply the already found contribution by 2 2229C If we are dealing with the pseudo B matrix contribution no contributions 2230C will be the same and we prepare to make a right density. 2231 IF (RSPTYP .EQ. 'F') THEN 2232 FACTKS=ONE 2233 LSAME=.FALSE. 2234 DENTYP='R' 2235 ELSE 2236 LSAME= (LISTC.EQ.LISTB .AND. IDLSTC.EQ.IDLSTB) 2237 IF (LSAME) THEN 2238 FACTKS=TWO 2239 ELSE 2240 FACTKS=ONE 2241 END IF 2242 END IF 2243 2244 CALL DAXPY(NT1AM(ISYCB),FACTKS,WORK(KETA1),1,RHO1,1) 2245 CALL DAXPY(NT2AM(ISYCB),FACTKS,WORK(KETA2),1,RHO2,1) 2246C 2247 IF (LOCDEB .OR. IPQMMM .GT. 14) THEN 2248 TAL1= DDOT(NT1AM(ISYCB),RHO1,1,RHO1,1) 2249 TAL2= DDOT(NT2AM(ISYCB),RHO2,1,RHO2,1) 2250 WRITE(LUPRI,*) 'Printing RHO: 2251 & Norm2 of singles: ', TAL1, 2252 & 'Norm2 of doubles: ', TAL2 2253 END IF 2254 2255C Then the second contribution ^BG*C. 2256 IF ( .NOT. LSAME) THEN 2257 2258C Initial parameters 2259 NC1 = NT1AM(ISYMTC) 2260 NC2 = NT2AM(ISYMTC) 2261 2262C Renaming (for readability) and resetting 2263 KC1=KB1 2264 KC2=KC1+NC1 2265 KDENS=KC2+NC2 2266 KGCMAT=KDENS+N2BST(ISYMTC) 2267 KWRK1=KGCMAT+N2BST(ISYMTC) 2268 LWRK1=LWORK-KWRK1 2269 2270C Zero the workspace 2271 CALL DZERO(WORK(KC1),NC1) 2272 CALL DZERO(WORK(KC2),NC2) 2273 CALL DZERO(WORK(KDENS),N2BST(ISYMTC)) 2274 CALL DZERO(WORK(KGCMAT),N2BST(ISYMTC)) 2275C 2276C Read C vector from file 2277 IOPT = 3 2278 CALL CC_RDRSP(LISTC,IDLSTC,ISYMTC,IOPT,MODEL, 2279 * WORK(KC1),WORK(KC2)) 2280C 2281C Print section 2282 IF ( IPQMMM .GT. 10 .OR. LOCDEB) THEN 2283 RHO1N = DDOT(NT1AM(ISYCB),RHO1,1,RHO1,1) 2284 RHO2N = DDOT(NT2AM(ISYCB),RHO2,1,RHO2,1) 2285 WRITE(LUPRI,*)'Norm of RHO1 in CCMM_QRTRANS on input,2:',RHO1N 2286 WRITE(LUPRI,*)'Norm of RHO2 in CCMM_QRTRANS on input,2:',RHO2N 2287 RHO1N = DDOT(NC1,WORK(KC1),1,WORK(KC1),1) 2288 RHO2N = DDOT(NC2,WORK(KC2),1,WORK(KC2),1) 2289 WRITE(LUPRI,*)'Norm af C1AM in CCMM_QRTRANS on input,2:',RHO1N 2290 WRITE(LUPRI,*)'Norm af C2AM in CCMM_QRTRANS on input,2:',RHO2N 2291 ENDIF 2292C---------------------------- 2293C Construct auxiliary densities 2294C---------------------------- 2295 CALL CCMM_D1AO(WORK(KDENS),WORK(KC1),WORK(KC2),MODEL,DENTYP, 2296 * LISTB,IDLSTB,ISYMTB, 2297 * WORK(KWRK1),LWRK1) 2298C 2299C---------------------------- 2300C Construct effective operator 2301C---------------------------- 2302 CALL CCMM_GEFF(WORK(KDENS),WORK(KGCMAT),WORK(KWRK1),LWRK1) 2303 2304C Transform the effective operator according to RSPTYP flag 2305C 2306C Dynamical memory allocation to store the response transformed operator 2307 KETA1=KWRK1 2308 KETA2=KETA1+NT1AM(ISYCB) 2309 KWRK2=KETA2+NT2AM(ISYCB) 2310 LWRK2=LWORK-KWRK2 2311 2312 IF (LWRK2 .LT. 0) THEN 2313 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KWRK2 2314 CALL QUIT( 'Too little work in CCMM_QRTRANSFORMER, 2') 2315 END IF 2316 LABEL='GIVE INT' 2317 IF (RSPTYP.EQ.'G') THEN 2318 CALL CCLR_FA(LABEL,ISYMTC,LISTB,IDLSTB,LISTL,0, 2319 * WORK(KGCMAT),WORK(KETA1),LWRK2) 2320 ELSE IF (RSPTYP.EQ.'B') THEN 2321 CALL CCCR_AA(LABEL,ISYMTC,LISTB,IDLSTB,WORK(KGCMAT), 2322 * WORK(KETA1),LWRK2) 2323 CALL CCLR_DIASCL(WORK(KETA2),TWO,ISYCB) 2324 ELSE IF ( (RSPTYP.EQ.'K') .OR. (RSPTYP.EQ.'F')) THEN 2325 CALL CC_ETAC(ISYMTC,LABEL,WORK(KETA1),LISTB,IDLSTB,0, 2326 * WORK(KGCMAT),WORK(KWRK2),LWRK2) 2327 DENTYP='Q' ! only effect for F term since K contains no more conts 2328 END IF 2329 2330C Print section 2331 IF(LOCDEB .OR. IPQMMM .GT. 14) THEN 2332 TAL1= DDOT(NT1AM(ISYCB),WORK(KETA1),1,WORK(KETA1),1) 2333 TAL2= DDOT(NT2AM(ISYCB),WORK(KETA2),1,WORK(KETA2),1) 2334 WRITE(LUPRI,*) 'Printing transformed G^C*B contribution. 2335 & Norm2 of singles: ', TAL1, 2336 & 'Norm2 of doubles: ', TAL2 2337 END IF 2338 2339 CALL DAXPY(NT1AM(ISYCB),ONE,WORK(KETA1),1,RHO1,1) 2340 CALL DAXPY(NT2AM(ISYCB),ONE,WORK(KETA2),1,RHO2,1) 2341C 2342 IF (LOCDEB .OR. IPQMMM .GT. 14) THEN 2343 TAL1= DDOT(NT1AM(ISYCB),RHO1,1,RHO1,1) 2344 TAL2= DDOT(NT2AM(ISYCB),RHO2,1,RHO2,1) 2345 WRITE(LUPRI,*) 'Printing RHO: 2346 & Norm2 of singles: ', TAL1, 2347 & 'Norm2 of doubles: ', TAL2 2348 END IF 2349 END IF 2350 2351C Finally the third contribution GBC 2352 IF ( (RSPTYP.EQ.'G') .OR. (RSPTYP.EQ.'B') 2353 & .OR. (RSPTYP.EQ.'F')) THEN 2354C Note we do not delete C vector from work 2355 KGBC=KDENS+N2BST(ISYCB) 2356 KWRK1=KGBC+N2BST(ISYCB) 2357 LWRK1=LWORK-KWRK1 2358 2359C Reset the primary workspace 2360 CALL DZERO(WORK(KDENS),N2BST(ISYCB)) 2361 CALL DZERO(WORK(KGBC),N2BST(ISYCB)) 2362 2363C Construct the quadratic CCMM auxiliary density 2364 IF ((RSPTYP .EQ. 'G') .OR. (RSPTYP .EQ. 'B') ) THEN 2365 CALL CCMM_D2AO(WORK(KDENS),ISYRES, 2366 * LISTB,IDLSTB,ISYMTB, 2367 * LISTC,IDLSTC,ISYMTC, 2368 * MODEL,WORK(KWRK1),LWRK1) 2369 ELSE IF (RSPTYP .EQ. 'F') THEN 2370 CALL CCMM_D1AO(WORK(KDENS),WORK(KC1),WORK(KC2),MODEL,DENTYP, 2371 * LISTB,IDLSTB,ISYMTB, 2372 * WORK(KWRK1),LWRK1) 2373 END IF 2374 2375C Construct effective operator 2376 IF(LOCDEB .OR. IPQMMM .GT. 14) THEN 2377 TAL1= DDOT(N2BST(ISYCB),WORK(KGBC),1,WORK(KGBC),1) 2378 WRITE(LUPRI,*) 'Printing {GBC} before build. 2379 & Norm2 of operator: ', TAL1 2380 END IF 2381 CALL CCMM_GEFF(WORK(KDENS),WORK(KGBC),WORK(KWRK1),LWRK1) 2382C 2383 IF(LOCDEB .OR. IPQMMM .GT. 14) THEN 2384 TAL1= DDOT(N2BST(ISYCB),WORK(KGBC),1,WORK(KGBC),1) 2385 WRITE(LUPRI,*) 'Printing {GBC} after build. 2386 & Norm2 of operator: ', TAL1 2387 END IF 2388 2389C Then do response transformation 2390 KETA1=KWRK1 2391 KETA2=KETA1+NT1AM(ISYCB) 2392 KWRK2=KETA2+NT2AM(ISYCB) 2393 LWRK2=LWORK-KWRK2 2394 IF (LWRK2 .LT. 0) THEN 2395 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KWRK2 2396 CALL QUIT( 'Too little work in CCMM_QRTRANSFORMER, 2') 2397 END IF 2398 CALL DZERO(WORK(KETA1),NT1AM(ISYCB)+NT2AM(ISYCB)) 2399 2400 LABEL='GIVE INT' 2401 IF ( (RSPTYP .EQ. 'G') .OR. (RSPTYP .EQ. 'F')) THEN 2402 CALL CC_ETAC(ISYCB,LABEL,WORK(KETA1),LISTL,0,0, 2403 * WORK(KGBC),WORK(KWRK2),LWRK2) 2404 ELSE IF (RSPTYP .EQ. 'B') THEN 2405 CALL CC_XKSI(WORK(KETA1),LABEL,ISYCB,0,WORK(KGBC), 2406 * WORK(KWRK2),LWRK2) 2407 CALL CCLR_DIASCL(WORK(KETA2),TWO,ISYCB) 2408 END IF 2409 2410 LOCDEB=.TRUE. 2411 IF(LOCDEB .OR. IPQMMM .GT. 14) THEN 2412 TAL1= DDOT(NT1AM(ISYCB),WORK(KETA1),1,WORK(KETA1),1) 2413 TAL2= DDOT(NT2AM(ISYCB),WORK(KETA2),1,WORK(KETA2),1) 2414 WRITE(LUPRI,*) 'Printing transformed G^{BC} contribution. 2415 & Norm2 of singles: ', TAL1, 2416 & 'Norm2 of doubles: ', TAL2 2417 END IF 2418 LOCDEB=.FALSE. 2419 2420 CALL DAXPY(NT1AM(ISYCB),ONE,WORK(KETA1),1,RHO1,1) 2421 CALL DAXPY(NT2AM(ISYCB),ONE,WORK(KETA2),1,RHO2,1) 2422 END IF 2423C 2424 CALL QEXIT('CCMM_QRTRANSFORMER') 2425 2426 END 2427C 2428C*********************************************************************** 2429C 2430C /* Deck ccmm_geff */ 2431 SUBROUTINE CCMM_GEFF(DENS,GEFF,WORK,LWORK) 2432C 2433C----------------------------------------------------------------------------- 2434C 2435C Construct effective operators from contracted densities. 2436C 2437C KS, aug 10 2438C----------------------------------------------------------------------------- 2439C 2440#include "implicit.h" 2441#include "maxorb.h" 2442#include "mxcent.h" 2443#include "qmmm.h" 2444#include "dummy.h" 2445#include "iratdef.h" 2446#include "priunit.h" 2447#include "ccorb.h" 2448#include "ccsdsym.h" 2449#include "ccsdio.h" 2450#include "ccsdinp.h" 2451#include "ccfield.h" 2452#include "exeinf.h" 2453#include "ccfdgeo.h" 2454#include "ccslvinf.h" 2455#include "qm3.h" 2456#include "ccinftap.h" 2457#include "nuclei.h" 2458#include "orgcom.h" 2459 2460 DIMENSION WORK(LWORK) 2461 DIMENSION DENS(*), GEFF(*) 2462 INTEGER NDIM 2463 LOGICAL LSKIP 2464 CHARACTER*8 LABEL 2465C 2466C Initial parameters 2467 NDIM=3*NNZAL 2468 2469C Workspace initialization 2470 KEPSAO = 1 2471 KEFIEL = KEPSAO + 3*NNBASX 2472 KINDIP = KEFIEL + NDIM 2473 KGAO = KINDIP + NDIM 2474 KWRK = KGAO + NNBASX 2475 LWRK = LWORK - KWRK 2476 IF (LWRK .LT. 0) THEN 2477 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KWRK 2478 CALL QUIT( 'Too little work in CCMM_GEFF') 2479 END IF 2480 2481C Convert auxiliary density to electric field 2482 CALL CCMM_DENS2F(DENS,WORK(KEFIEL),WORK(KWRK),LWRK) 2483 2484C Transform the field to induced dipoles 2485 CALL CCMM_F2DIP(WORK(KEFIEL),WORK(KINDIP),WORK(KWRK),LWRK) 2486 2487 CALL DZERO(WORK(KGAO),NNBASX) 2488 2489C Construct response G from the induced moments - eg. G^C=-sum_a ∆mu^C eps_a 2490 LRI = 0 2491 2492 DO 100 I=1,MMCENT 2493 2494 LSKIP = .FALSE. 2495 2496 CALL CCMM_EPSAO(WORK(KEPSAO),I,LSKIP,WORK(KWRK),LWRK) 2497 2498 IF (LSKIP) THEN 2499 IF (ZEROAL(I) .EQ. -1) THEN 2500 GOTO 100 2501 ELSE 2502 LRI = LRI + 3 2503 GOTO 100 2504 END IF 2505 END IF 2506C 2507 FACx = -WORK(KINDIP + LRI + 0) 2508 2509 FACy = -WORK(KINDIP + LRI + 1) 2510 2511 FACz = -WORK(KINDIP + LRI + 2) 2512 2513 2514 CALL DAXPY(NNBASX,FACx,WORK(KEPSAO), 2515 * 1,WORK(KGAO),1) 2516 2517 CALL DAXPY(NNBASX,FACy,WORK(KEPSAO+NNBASX), 2518 * 1,WORK(KGAO),1) 2519 2520 CALL DAXPY(NNBASX,FACz,WORK(KEPSAO+2*NNBASX), 2521 * 1,WORK(KGAO),1) 2522 2523 LRI = LRI + 3 2524 100 CONTINUE 2525 2526C Unpack the integrals and store in GEFF 2527 LABEL = 'GIVE INT' 2528 CALL CCMMINT(LABEL,LUQMMM,GEFF,WORK(KGAO), 2529 & IRREP,ISYM,IERR,WORK(KWRK),LWRK) 2530 2531 END 2532C*********************************************************************** 2533C /* Deck ccmm_d2ao */ 2534 SUBROUTINE CCMM_D2AO(AODEN,ISYRES, 2535 * LISTB,IDLSTB,ISYMB, 2536 * LISTC,IDLSTC,ISYMC, 2537 * MODEL,WORK,LWORK) 2538C 2539C Purpose: To calculate contributions to the one electron 2540C quadratic auxiliary density matrix to be used 2541C in CCMM quadratic response theory. 2542C 2543C 2544C Current models: CCSD 2545C 2546#include "implicit.h" 2547#include "priunit.h" 2548#include "dummy.h" 2549#include "maxash.h" 2550#include "maxorb.h" 2551#include "mxcent.h" 2552#include "aovec.h" 2553#include "iratdef.h" 2554#include "ccorb.h" 2555#include "infind.h" 2556#include "blocks.h" 2557#include "ccsdinp.h" 2558#include "ccsdsym.h" 2559#include "ccexci.h" 2560#include "ccsdio.h" 2561#include "distcl.h" 2562#include "cbieri.h" 2563#include "eribuf.h" 2564#include "cclr.h" 2565C 2566 PARAMETER (ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, TWO = 2.0D0) 2567 DIMENSION AODEN(*), WORK(LWORK) 2568 CHARACTER*(10) MODDUM, MODEL 2569 CHARACTER*(2) LISTB, LISTC 2570 CHARACTER DENTYP*(1) 2571 INTEGER ISYRES, IDLSTB, IDLSTC, ISYMB, ISYMC 2572 LOGICAL LOCDEB 2573 PARAMETER (LOCDEB=.FALSE.) 2574 2575 CALL QENTER('CCMM_D2AO') 2576 2577 ISYML = ILSTSYM('L0',1) ! Symmetry of multipliers 2578 ISYMR = ILSTSYM('R0',1) ! Symmetry of amplitudes 2579C----------------------------------- 2580C Initial work space allocation. 2581C----------------------------------- 2582C 2583 KONEIA=1 2584 KL1AM =KONEIA+NT1AMX 2585 KL2AM =KL1AM +NT1AMX 2586 KB1 =KL2AM +NT2SQ(ISYML) 2587 KB2 =KB1 +NT1AM(ISYMB) 2588 KC1 =KB2 +NT2AM(ISYMB) 2589 KC2 =KC1 +NT1AM(ISYMC) 2590 KWRK1 =KC2 +NT2AM(ISYMC) 2591 LWRK1 =LWORK -KWRK1 2592 IF (LWRK1 .LT. 0) THEN 2593 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KWRK1 2594 CALL QUIT('CCMM_D2AO: Not enough memory') 2595 ENDIF 2596 2597C----------------------------------- 2598C Reset the working matrices 2599C----------------------------------- 2600 CALL DZERO(WORK(KONEIA),NT1AMX) 2601 CALL DZERO(WORK(KL1AM),NT1AM(ISYML)) 2602 CALL DZERO(WORK(KL2AM),NT2SQ(ISYML)) 2603 CALL DZERO(WORK(KB1),NT1AM(ISYMB)) 2604 CALL DZERO(WORK(KB2),NT2AM(ISYMB)) 2605 CALL DZERO(WORK(KC1),NT1AM(ISYMC)) 2606 CALL DZERO(WORK(KC2),NT2AM(ISYMC)) 2607C----------------------------------- 2608C Read in multipliers and singles part of trial vectors 2609C----------------------------------- 2610C WRITE(LUPRI,*) 'NT1AMX: ', NT1AMX 2611C WRITE(LUPRI,*) 'NT2AMX: ', NT2AMX 2612C WRITE(LUPRI,*) 'NT2AM(ISYML): ', NT2AM(ISYML) 2613C WRITE(LUPRI,*) 'NB1: ', NT1AM(ISYMB) 2614C WRITE(LUPRI,*) 'NB2: ', NT2AM(ISYMB) 2615C WRITE(LUPRI,*) 'NC1: ', NT1AM(ISYMC) 2616C WRITE(LUPRI,*) 'NC2: ', NT2AM(ISYMC) 2617 2618 IOPT = 3 2619 CALL CC_RDRSP('L0',1,ISYML,IOPT,MODDUM, 2620 * WORK(KL1AM),WORK(KWRK1)) 2621C CALL DZERO(WORK(KL1AM),NT1AM(ISYMB)) 2622C CALL DZERO(WORK(KWRK1),NT2AM(ISYMB)) 2623 2624 CALL CC_RDRSP(LISTB,IDLSTB,ISYMB,IOPT,MODDUM, 2625 * WORK(KB1),WORK(KB2)) 2626 CALL CC_RDRSP(LISTC,IDLSTC,ISYMC,IOPT,MODDUM, 2627 * WORK(KC1),WORK(KC2)) 2628 2629C scale doubles diagonal 2630 IF (ISYMB .EQ. 1) CALL CCLR_DIASCL(WORK(KB2),TWO,ISYMB) 2631 IF (ISYMC .EQ. 1) CALL CCLR_DIASCL(WORK(KC2),TWO,ISYMC) 2632 2633C Square up the L0_2 block. 2634 CALL CC_T2SQ(WORK(KWRK1),WORK(KL2AM),ISYML) 2635 2636 IF (LOCDEB) THEN 2637 TAL1 = DDOT(NT1AMX,WORK(KL1AM),1,WORK(KL1AM),1) 2638 TAL2 = DDOT(NT2SQ(ISYML),WORK(KL2AM),1,WORK(KL2AM),1) 2639 WRITE(LUPRI,*) 'CCMM_D2AO: Norm of L0_1: ', TAL1, 2640 & 'Norm of L0_2: ', TAL2 2641 TAL1 = DDOT(NT1AMX,WORK(KB1),1,WORK(KB1),1) 2642 TAL2 = DDOT(NT2AMX,WORK(KB2),1,WORK(KB2),1) 2643 WRITE(LUPRI,*) 'CCMM_D2AO: Norm of B_1: ', TAL1, 2644 & 'Norm of B_2: ', TAL2 2645 TAL1 = DDOT(NT1AMX,WORK(KC1),1,WORK(KC1),1) 2646 TAL2 = DDOT(NT2AMX,WORK(KC2),1,WORK(KC2),1) 2647 WRITE(LUPRI,*) 'CCMM_D2AO: Norm of C_1: ', TAL1, 2648 & 'Norm of C_2: ', TAL2 2649 END IF 2650 2651C---------------------------------------------------------- 2652C Construct <L2|[[Eia,B1],C2]|HF> contribution to DXia. 2653C---------------------------------------------------------- 2654 IF (LOCDEB) THEN 2655 TAL1 = DDOT(NT1AMX,WORK(KONEIA),1,WORK(KONEIA),1) 2656 WRITE(LUPRI,*) 'CCMM_D2AO: Norm of DIA: ', TAL1 2657 END IF 2658C 2659 CALL CC_DXIA21(WORK(KONEIA),WORK(KL2AM),ISYML, 2660 * WORK(KB1),ISYMB,WORK(KC2),ISYMC, 2661 * WORK(KWRK1),LWRK1) 2662 IF (LOCDEB) THEN 2663 TAL1 = DDOT(NT1AMX,WORK(KONEIA),1,WORK(KONEIA),1) 2664 WRITE(LUPRI,*) 'CCMM_D2AO: Norm of DIA after 1: ', TAL1 2665 END IF 2666 2667C 2668C---------------------------------------------------------- 2669C Construct <L2|[[Eia,C1],B2]|HF> contribution to DXia. 2670C---------------------------------------------------------- 2671 CALL CC_DXIA21(WORK(KONEIA),WORK(KL2AM),ISYML, 2672 * WORK(KC1),ISYMC,WORK(KB2),ISYMB, 2673 * WORK(KWRK1),LWRK1) 2674 IF (LOCDEB) THEN 2675 TAL1 = DDOT(NT1AMX,WORK(KONEIA),1,WORK(KONEIA),1) 2676 WRITE(LUPRI,*) 'CCMM_D2AO: Norm of DIA after 2: ', TAL1 2677 END IF 2678 2679 KONEAI=KWRK1 2680 KONEAB=KONEAI+NT1AMX 2681 KONEIJ=KONEAB+NMATAB(1) 2682 KWRK2=KONEIJ+NMATIJ(1) 2683 LWRK2=LWORK-KWRK2 2684 IF (LWRK2 .LT. 0) THEN 2685 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KWRK2 2686 CALL QUIT('Insufficient memory for 3. alloc. in CCMM_D1AO') 2687 ENDIF 2688 CALL DZERO(WORK(KONEAI),NT1AMX) 2689 CALL DZERO(WORK(KONEAB),NMATAB(1)) 2690 CALL DZERO(WORK(KONEIJ),NMATIJ(1)) 2691C 2692C---------------------------------------------------------- 2693C Construct <L1|[[Eia,B1],C1]|HF> contribution to DXia. 2694C---------------------------------------------------------- 2695C We do this by doing two contractions of the type 2696C sum_ia t_ia b_iq c_pa = sum_i zeta_ip b_iq = D_pq (=D_ia) 2697C For efficiency we start by performing the sum over virtual orbitals 2698 CALL CC_DXIJ(WORK(KONEIJ),WORK(KL1AM),ISYML,WORK(KB1),ISYMB) 2699 CALL CCMM_DIA(WORK(KONEIA),WORK(KONEIJ),ISYRES,WORK(KC1),ISYMC) 2700 2701 IF (LOCDEB) THEN 2702 TAL1 = DDOT(NT1AMX,WORK(KONEIA),1,WORK(KONEIA),1) 2703 WRITE(LUPRI,*) 'CCMM_D2AO: Norm of DIA after 3: ', TAL1 2704 END IF 2705C and the related contribution 2706 CALL DZERO(WORK(KONEIJ),NMATIJ(1)) 2707 2708 CALL CC_DXIJ(WORK(KONEIJ),WORK(KL1AM),ISYML,WORK(KC1),ISYMC) 2709 CALL CCMM_DIA(WORK(KONEIA),WORK(KONEIJ),ISYRES,WORK(KB1),ISYMB) 2710 2711 IF (LOCDEB) THEN 2712 TAL1 = DDOT(NT1AMX,WORK(KONEIA),1,WORK(KONEIA),1) 2713 WRITE(LUPRI,*) 'CCMM_D2AO: Norm of DIA after 4: ', TAL1 2714 END IF 2715 2716 KT1AM = KWRK2 2717 KLAMDP = KT1AM + NT1AMX 2718 KLAMDH = KLAMDP + NLAMDT 2719 KWRK3 = KLAMDH + NLAMDT 2720 LWRK3 = LWORK - KWRK3 2721 IF (LWRK3 .LT. 0) THEN 2722 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KWRK3 2723 CALL QUIT('Insufficient memory for 2. alloc. in CCMM_D1AO') 2724 ENDIF 2725 2726 CALL DZERO(WORK(KT1AM),NT1AMX) 2727 IOPT = 1 2728 2729C Read singles to be used for construction of Lambda matrices 2730 KDUM = KWRK3 2731 CALL CC_RDRSP('R0',1,ISYMR,IOPT,MODDUM, 2732 * WORK(KT1AM),WORK(KDUM)) 2733 2734 CALL LAMMAT(WORK(KLAMDP),WORK(KLAMDH),WORK(KT1AM),WORK(KWRK3), 2735 * LWRK3) 2736C 2737C-------------------------------------------------------- 2738C Add the one electron density in the AO-basis. 2739C-------------------------------------------------------- 2740CC 2741 CALL DZERO(WORK(KONEIJ),NMATIJ(1)) 2742 2743 ISDEN = 1 2744 CALL CC_DENAO(AODEN,ISDEN,WORK(KONEAI),WORK(KONEAB), 2745 * WORK(KONEIJ),WORK(KONEIA),ISDEN,WORK(KLAMDP),1, 2746 * WORK(KLAMDH),1,WORK(KWRK3),LWRK3) 2747C 2748 CALL QEXIT('CCMM_D2AO') 2749 2750 END 2751 2752C*********************************************************************** 2753C /* Deck ccmm_addgdiff */ 2754 SUBROUTINE CCMM_ADDGDIFF(FOCK,WORK,LWORK) 2755C************************************************************** 2756C 2757C Purpose: For CC2 fullfld Construct the EFFective Diff G solvent operator FROM 2758C the diff between cc induced dipoles and the HF INDUCED DIPOLES and add to Fock 2759C (in AO basis) 2760C 2761C 2762C Get G operators from file dumped 2763C Subtract G^CC-G^HF 2764C Add effective operator G to FOCK 2765C 2766C KS 11 2767C************************************************************** 2768C 2769#include "implicit.h" 2770#include "maxorb.h" 2771#include "mxcent.h" 2772#include "qmmm.h" 2773#include "dummy.h" 2774#include "iratdef.h" 2775#include "priunit.h" 2776#include "ccorb.h" 2777#include "ccsdsym.h" 2778#include "ccsdio.h" 2779#include "ccsdinp.h" 2780#include "ccfield.h" 2781#include "exeinf.h" 2782#include "ccfdgeo.h" 2783#include "qm3.h" 2784#include "ccinftap.h" 2785#include "nuclei.h" 2786#include "orgcom.h" 2787C 2788 PARAMETER (ZERO = 0.0D0, MONE=-1.0D0 ,ONE = 1.0D0, 2789 & IZERO = 0 , TWO = 2.0D0) 2790 PARAMETER (D3I = 1.0D0/3.0D0, D6I = 1.0D0/6.0D0) 2791 DIMENSION WORK(LWORK),FOCK(*) 2792 DIMENSION INTREP(9*MXCENT), INTADR(9*MXCENT) 2793 CHARACTER*8 LABINT(9*MXCENT) 2794 LOGICAL TOFILE,TRIMAT,EXP1VL,LOCDEB 2795 PARAMETER (LOCDEB=.FALSE.) 2796C 2797 CHARACTER*8 LABEL 2798C 2799C 2800 CALL QENTER('CCMM_ADDGDIFF') 2801C------------------------- 2802C Init parameters 2803C------------------------- 2804 NDIM=3*NNZAL 2805C-------------------------- 2806C Dynamical allocation 2807C-------------------------- 2808C 2809 KCC = 1 2810 KHF = KCC + N2BST(ISYMOP) 2811 KWRK = KHF + N2BST(ISYMOP) 2812 LWRK1 = LWORK - KWRK 2813 2814 IF (LWRK1.LT.0) CALL QUIT( 'Too little work in CCMM_ADDGDIFF' ) 2815 2816 CALL DZERO(WORK(KCC),N2BST(ISYMOP)) 2817 CALL DZERO(WORK(KHF),N2BST(ISYMOP)) 2818C 2819 2820 WRITE(LUPRI,*) 'Reading G from file' 2821 CALL GET_FROM_FILE('G_MATRIX',N2BST(ISYMOP),WORK(KCC)) 2822 2823 WRITE(LUPRI,*) 'Reading GHF from file' 2824 CALL GET_FROM_FILE('G_MATRIXHF',N2BST(ISYMOP),WORK(KHF)) 2825 2826 2827 IF (IPQMMM .GT.14 ) THEN 2828 CALL AROUND( 'CCMM cont. to AO matrix' ) 2829 CALL OUTPUT(WORK(KCC),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI) 2830 ENDIF 2831C 2832 IF (IPQMMM .GT.14) THEN 2833 CALL AROUND( 'Usual Fock AO matrix' ) 2834 CALL CC_PRFCKAO(FOCK,ISYMOP) 2835 ENDIF 2836 2837 CALL DAXPY(N2BST(ISYMOP),-ONE,WORK(KHF),1,WORK(KCC),1) 2838 2839 TAL=DDOT(N2BST(ISYMOP),WORK(KCC),1,WORK(KCC),1) 2840 WRITE(LUPRI,*) 'Difference density cont norm ', TAL 2841C 2842 CALL DAXPY(N2BST(ISYMOP),ONE,WORK(KCC),1,FOCK,1) 2843 2844 IF (IPQMMM .GT.14) THEN 2845 CALL AROUND( 'Modified Fock AO matrix when leaving ADDGDIFF') 2846 CALL CC_PRFCKAO(FOCK,ISYMOP) 2847 ENDIF 2848C 2849 CALL QEXIT('CCMM_ADDGDIFF') 2850 END 2851