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 19*=====================================================================* 20 SUBROUTINE CCXIMCON(IXIM,NXIM,LISTQ,IXDOTS,XCONS, 21 & MXVEC,WORK,LWORK) 22*---------------------------------------------------------------------* 23* 24* Purpose: calculate contributions from the generalized relaxation 25* & reorthogonalization X intermediate to response 26* functions involving orbital relaxation or perturbation- 27* dependent basis sets. 28* 29* IXIM -- array with the perturbation indeces 30* NXIM -- length of IXIM 31* LISTQ -- type of vectors the X intermediates are to 32* be dotted on 33* IXDOTS -- matrix of indeces for the vectors with which 34* the dot products are to be calculated 35* XCONS -- matrix with the dot product results 36* MXVEC -- leading dimension of IXDOTS and XCONS 37* 38* Christof Haettig 1-6-1999 39* 40* N.B.: this routine is not yet adapted for QMATP diff. from QMATH 41* 42*---------------------------------------------------------------------* 43 IMPLICIT NONE 44#include "priunit.h" 45#include "dummy.h" 46#include "ccorb.h" 47#include "ccsdsym.h" 48#include "ccexpfck.h" 49#include "cc1dxfck.h" 50#include "ccr1rsp.h" 51#include "ccfro.h" 52#include "ccroper.h" 53#include "inftap.h" 54#include "iratdef.h" 55 56 LOGICAL LOCDBG 57 PARAMETER (LOCDBG = .FALSE.) 58 59 INTEGER ISYM0, LUFCK 60 PARAMETER( ISYM0 = 1) 61 CHARACTER LABEL0*(8) 62 PARAMETER( LABEL0 = 'HAM0 ' ) 63 64 LOGICAL LORX, LFOCK0, LORXQ, LPDBS 65 CHARACTER*(3) LISTQ 66 INTEGER ISYXIM, LWORK, NXIM, MXVEC 67 68 INTEGER IXIM(NXIM), IXDOTS(MXVEC,NXIM) 69 70#if defined (SYS_CRAY) 71 REAL FREQ, XCONS(MXVEC,NXIM), WORK(LWORK) 72 REAL HALF, ONE, TWO, ZERO, FTRACE 73#else 74 DOUBLE PRECISION FREQ, XCONS(MXVEC,NXIM), WORK(LWORK) 75 DOUBLE PRECISION HALF, ONE, TWO, ZERO, FTRACE 76#endif 77 PARAMETER( ONE = 1.0D0, TWO = 2.0D0, ZERO = 0.0D0, HALF = 0.5D0 ) 78 79 CHARACTER*(10) MODEL 80 CHARACTER*(8) LABEL, LABELQ 81 LOGICAL NOKAPPA 82 INTEGER NCMOT, NASHT, N2ASHX, LCINDX 83 INTEGER IFOCK, IEXPV, IADRF, KOVERLP, KFOCK0, KFOCK1, KCMO, KAPPA 84 INTEGER IKAPPA, IOPT, IOPER, ISYM, IDXXIM, KOFF1, IFCK1 85 INTEGER KEND1,LWRK1,KEND2,LWRK2,KXIM,KSCR1,KRMAT,KQMATP,KQMATH 86 INTEGER IVEC, IFILE, ISYMQ, IOPERQ, KEND3, LWRK3 87 INTEGER ISYM1, ISYM2, KOFF2, KQTRP, IREAL 88 89* external functions: 90 INTEGER IEFFFOCK 91 INTEGER IEXPECT 92 INTEGER IROPER 93 INTEGER I1DXFCK 94 INTEGER ILSTSYMRLX 95#if defined (SYS_CRAY) 96 REAL DDOT 97#else 98 DOUBLE PRECISION DDOT 99#endif 100 101*---------------------------------------------------------------------* 102* check, if there is anything at all to do: 103*---------------------------------------------------------------------* 104 105 IF (NXIM.LE.0) RETURN 106 107*---------------------------------------------------------------------* 108* get some constants from sirius common block: 109*---------------------------------------------------------------------* 110 111 CALL CC_SIRINF(NCMOT,NASHT,N2ASHX,LCINDX) 112 113*---------------------------------------------------------------------* 114* allocate memory for perturbation-independent stuff needed: 115*---------------------------------------------------------------------* 116 KCMO = 1 117 KFOCK0 = KCMO + NCMOT 118 KOVERLP = KFOCK0 + N2BST(ISYM0) 119 KEND1 = KOVERLP + N2BST(ISYM0) 120 LWRK1 = LWORK - KEND1 121 122 IF (LWRK1.LT.0) THEN 123 CALL QUIT('Insufficient work space in CCXIMCON.') 124 END IF 125 126*---------------------------------------------------------------------* 127* read MO coefficients from file: 128*---------------------------------------------------------------------* 129 130 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY, 131 & .FALSE.) 132 REWIND LUSIFC 133 CALL MOLLAB('SIR IPH ',LUSIFC,LUPRI) 134 READ (LUSIFC) 135 READ (LUSIFC) 136 CALL READI(LUSIFC,IRAT*NCMOT,WORK(KCMO)) 137 CALL GPCLOSE(LUSIFC,'KEEP') 138 139*---------------------------------------------------------------------* 140* read overlap matrix from file: 141*---------------------------------------------------------------------* 142 143 IF (LWRK1.LT.NBAST) THEN 144 CALL QUIT('Insufficient work space in CCXIMCON.') 145 END IF 146 147 CALL RDONEL('OVERLAP ',.TRUE.,WORK(KEND1),NBAST) 148 CALL CCSD_SYMSQ(WORK(KEND1),ISYM0,WORK(KOVERLP)) 149 150*---------------------------------------------------------------------* 151* read zeroth-order effective Fock matrix if available: 152*---------------------------------------------------------------------* 153 call quit('hjaaj: code error, ISYM not defined, please fix!') 154 IFOCK = IEFFFOCK(LABEL0,ISYM,1) 155 156 IF (LEXPFCK(2,IFOCK)) THEN 157 IADRF = IADRFCK(1,IFOCK) 158 159 LUFCK = -1 160 CALL WOPEN2(LUFCK,FILFCKEFF,64,0) 161 CALL GETWA2(LUFCK,FILFCKEFF,WORK(KFOCK0),IADRF,N2BST(ISYM0)) 162 CALL WCLOSE2(LUFCK,FILFCKEFF,'KEEP') 163 164 CALL CC_EFFCKMO(WORK(KFOCK0),ISYM0,WORK(KCMO),WORK(KOVERLP), 165 & WORK(KEND1),LWRK1) 166 167 LFOCK0 = .TRUE. 168 169 IF (LOCDBG) THEN 170 WRITE (LUPRI,*) 171 & 'CCXIMCON> zeroth-order effective Fock in ao/AO:' 172 CALL CC_PRONELAO(WORK(KFOCK0),ISYM0) 173 END IF 174 175 ELSE 176 LFOCK0 = .FALSE. 177 END IF 178 179*---------------------------------------------------------------------* 180* start loop over X intermediates: 181*---------------------------------------------------------------------* 182 DO IDXXIM = 1, NXIM 183 184 IKAPPA = IXIM(IDXXIM) 185 LABEL = LRTHFLBL(IKAPPA) 186 IOPER = IROPER(LABEL,ISYXIM) 187 LORX = .TRUE. 188 LPDBS = LPDBSOP(IOPER) 189 ISYXIM = ILSTSYMRLX('R1',IKAPPA) 190 191 IF (LOCDBG) THEN 192 WRITE (LUPRI,*) 'CCXIMCON> IDXXIM:',IDXXIM 193 WRITE (LUPRI,*) 'CCXIMCON> IKAPPA:',IKAPPA 194 WRITE (LUPRI,*) 'CCXIMCON> LABEL :',LABEL 195 WRITE (LUPRI,*) 'CCXIMCON> LORX :',LORX 196 END IF 197 198 KXIM = KEND1 199 KRMAT = KXIM + N2BST(ISYXIM) 200 KQMATP = KRMAT + N2BST(ISYXIM) 201 KQMATH = KQMATP + N2BST(ISYXIM) 202 KQTRP = KQMATH + N2BST(ISYXIM) 203 KAPPA = KQTRP + N2BST(ISYXIM) 204 KEND2 = KAPPA + 2*NALLAI(ISYXIM) 205 206 LWRK2 = LWORK - KEND2 207 208 IF (LWRK2 .LT. 0) THEN 209 CALL QUIT('Insufficient memory in CCXIMCON.') 210 END IF 211 212 ! read first-order effective Fock matrix from file 213 call quit('hjaaj: code error, ISYM not defined, please fix!') 214 IFOCK = IEFFFOCK(LABEL,ISYM,1) 215 IEXPV = IEXPECT(LABEL,ISYM) 216 IADRF = IADRFCK(1,IFOCK) 217 218 CALL WOPEN2(LUFCK,FILFCKEFF,64,0) 219 CALL GETWA2(LUFCK,FILFCKEFF,WORK(KXIM),IADRF,N2BST(ISYXIM)) 220 CALL WCLOSE2(LUFCK,FILFCKEFF,'KEEP') 221 222 IF (LOCDBG) THEN 223 FTRACE = ZERO 224 IF (ISYXIM.EQ.1) THEN 225 DO ISYM = 1, NSYM 226 KOFF1 = KXIM + IAODIS(ISYM,ISYM) 227 DO I = 1, NBAS(ISYM) 228 FTRACE = FTRACE + WORK(KOFF1+NBAS(ISYM)*(I-1)+I-1) 229 END DO 230 END DO 231 END IF 232 WRITE (LUPRI,*) 'LABEL:',LABEL 233 WRITE (LUPRI,*) 'ISYXIM,IFOCK,IEXPV:',ISYXIM,IFOCK,IEXPV 234 WRITE (LUPRI,*) 'FTRACE of read matrix:',FTRACE 235 WRITE (LUPRI,*) 'one-electron expect:',EXPVALUE(1,IEXPV) 236 WRITE (LUPRI,*) 'two-electron expect:',EXPVALUE(2,IEXPV) 237 END IF 238 239 ! transform effective Fock matrix to MO basis 240 CALL CC_EFFCKMO(WORK(KXIM),ISYXIM,WORK(KCMO),WORK(KOVERLP), 241 & WORK(KEND2),LWRK2) 242 243 IF (LOCDBG) THEN 244 FTRACE = ZERO 245 IF (ISYXIM.EQ.1) THEN 246 DO ISYM = 1, NSYM 247 KOFF1 = KXIM + IAODIS(ISYM,ISYM) 248 DO I = 1, NBAS(ISYM) 249 FTRACE = FTRACE + WORK(KOFF1+NBAS(ISYM)*(I-1)+I-1) 250 END DO 251 END DO 252 END IF 253 WRITE (LUPRI,*) 'LABEL:',LABEL 254 WRITE (LUPRI,*) 'ISYXIM:',ISYXIM 255 WRITE (LUPRI,*) 256 & 'FTRACE of matrix generated in CC_EFCKMO:',FTRACE 257 END IF 258 259 ! scale with 2 to get contribution to X intermediate 260 CALL DSCAL(N2BST(ISYXIM),TWO,WORK(KXIM),1) 261 262 IF (LOCDBG) THEN 263 WRITE (LUPRI,*) 264 & 'CCXIMCON> direct contribution to X intermediate:' 265 CALL CC_PRONELAO(WORK(KXIM),ISYXIM) 266 END IF 267 268 IF (LORX .OR. LPDBS) THEN 269 270 KSCR1 = KEND2 271 KFOCK1 = KSCR1 + N2BST(ISYXIM) 272 KEND3 = KFOCK1 + N2BST(ISYXIM) 273 LWRK3 = LWORK - KEND3 274 275 IF (LWRK3 .LT. 0) THEN 276 CALL QUIT('Insufficient memory in CCXIMCON.') 277 END IF 278 279 IF (LORX) THEN 280 CALL CC_RDHFRSP('R1 ',IKAPPA,ISYXIM,WORK(KAPPA)) 281 ELSE 282 CALL DZERO(WORK(KAPPA),2*NALLAI(ISYXIM)) 283 END IF 284 285 IOPER = IROPER(LABEL,ISYXIM) 286 CALL CC_GET_RMAT(WORK(KRMAT),IOPER,1,ISYXIM, 287 & WORK(KEND3),LWRK3) 288 289 NOKAPPA = .FALSE. 290 IREAL = ISYMAT(IOPER) 291 CALL CC_QMAT(WORK(KQMATP),WORK(KQMATH), 292 & WORK(KRMAT),WORK(KAPPA), 293 & IREAL,ISYXIM,NOKAPPA, 294 & WORK(KCMO),WORK(KEND3),LWRK3) 295 296 DO ISYM1 = 1, NSYM 297 ISYM2 = MULD2H(ISYM1,ISYXIM) 298 KOFF1 = KQMATP + IAODIS(ISYM1,ISYM2) 299 KOFF2 = KQTRP + IAODIS(ISYM2,ISYM1) 300 CALL TRSREC(NBAS(ISYM1),NBAS(ISYM2), 301 & WORK(KOFF1),WORK(KOFF2)) 302 END DO 303 CALL DSCAL(N2BST(ISYXIM),-ONE,WORK(KQTRP ),1) 304 305 CALL CC_MMOMMO('N','N',ONE,WORK(KFOCK0),ISYM0, 306 & WORK(KQTRP ),ISYXIM,ZERO,WORK(KSCR1),ISYXIM) 307 308 IF (LOCDBG) THEN 309 WRITE (LUPRI,*) 310 & 'CCXIMCON> relax. contrib. 1 to X intermediate:' 311 CALL CC_PRONELAO(WORK(KSCR1),ISYXIM) 312 FTRACE = ZERO 313 DO ISYM = 1, NSYM 314 KOFF1 = KSCR1 + IAODIS(ISYM,ISYM) 315 DO I = 1, NBAS(ISYM) 316 FTRACE = FTRACE + WORK(KOFF1+NBAS(ISYM)*(I-1)+I-1) 317 END DO 318 END DO 319 WRITE (LUPRI,*) 'trace:',FTRACE 320 END IF 321 322 ! add to result matrix: 323 CALL DAXPY(N2BST(ISYXIM),ONE,WORK(KSCR1),1,WORK(KXIM),1) 324 325 ! read contribution from 'one-index' transformed density 326 ! from file 327 FREQ = FRQLRTHF(IKAPPA) 328 IFCK1 = I1DXFCK('HAM0 ','R1 ',LABEL,FREQ,ISYXIM) 329 IADRF = IADR1DXF(1,IFCK1) 330 CALL WOPEN2(LUFCK,FIL1DXFCK,64,0) 331 CALL GETWA2(LUFCK,FIL1DXFCK,WORK(KFOCK1),IADRF,N2BST(ISYXIM)) 332 CALL WCLOSE2(LUFCK,FIL1DXFCK,'KEEP') 333 334 CALL CC_EFFCKMO(WORK(KFOCK1),ISYXIM,WORK(KCMO),WORK(KOVERLP), 335 & WORK(KEND3),LWRK3) 336 337 IF (LOCDBG) THEN 338 WRITE (LUPRI,*) 339 & 'CCXIMCON> relax. contrib. 2 to X intermediate:' 340 CALL CC_PRONELAO(WORK(KFOCK1),ISYXIM) 341 FTRACE = ZERO 342 DO ISYM = 1, NSYM 343 KOFF1 = KFOCK1 + IAODIS(ISYM,ISYM) 344 DO I = 1, NBAS(ISYM) 345 FTRACE = FTRACE + WORK(KOFF1+NBAS(ISYM)*(I-1)+I-1) 346 END DO 347 END DO 348 WRITE (LUPRI,*) 'trace:',FTRACE 349 END IF 350 351 ! add to result matrix: 352 CALL DAXPY(N2BST(ISYXIM),ONE,WORK(KFOCK1),1,WORK(KXIM),1) 353 354 END IF 355 356 IF (LOCDBG) THEN 357 WRITE (LUPRI,*) 'CCXIMCON> final result for X intermediate:' 358 CALL CC_PRONELAO(WORK(KXIM),ISYXIM) 359 FTRACE = ZERO 360 DO ISYM = 1, NSYM 361 KOFF1 = KXIM-1 + IAODIS(ISYM,ISYM) 362 DO I = 1, NBAS(ISYM) 363 FTRACE = FTRACE + WORK(KOFF1+NBAS(ISYM)*(I-1)+I) 364 END DO 365 END DO 366 WRITE (LUPRI,*) 'trace of X intermediate:',FTRACE 367 END IF 368 369*---------------------------------------------------------------------* 370* calculate all required dot products with this X intermediate: 371*---------------------------------------------------------------------* 372 IVEC = 1 373 374 DO WHILE (IXDOTS(IVEC,IDXXIM).NE.0 .AND. IVEC.LE.MXVEC) 375 376 377 IFILE = IXDOTS(IVEC,IDXXIM) 378 ISYMQ = ILSTSYMRLX(LISTQ,IFILE) 379 380 IF (ISYMQ.NE.ISYXIM) THEN 381 WRITE (LUPRI,*) IDXXIM, ISYXIM, LISTQ, IFILE, ISYMQ 382 CALL QUIT('symmetry mismatch in CCXIMCON.') 383 END IF 384 385 IF (LISTQ(1:3).EQ.'R1 ') THEN 386 LORXQ = .TRUE. 387 LABELQ = LRTHFLBL(IFILE) 388 ELSE 389 CALL QUIT('Unknown LISTQ in CCXIMCON.') 390 END IF 391 392 IOPERQ = IROPER(LABELQ,ISYMQ) 393 IREAL = ISYMAT(IOPERQ) 394 395 IF (LORXQ) THEN 396 CALL CC_RDHFRSP(LISTQ,IFILE,ISYMQ,WORK(KAPPA)) 397 ELSE 398 CALL DZERO(WORK(KAPPA),2*NALLAI(ISYMQ)) 399 END IF 400 401 CALL CC_GET_RMAT(WORK(KRMAT),IOPERQ,1,ISYMQ, 402 & WORK(KEND2),LWRK2) 403 404 NOKAPPA = .FALSE. 405 CALL CC_QMAT(WORK(KQMATP),WORK(KQMATH), 406 & WORK(KRMAT),WORK(KAPPA), 407 & IREAL,ISYMQ,NOKAPPA, 408 & WORK(KCMO),WORK(KEND2),LWRK2) 409 410 DO ISYM1 = 1, NSYM 411 ISYM2 = MULD2H(ISYM1,ISYMQ) 412 KOFF1 = KQMATH + IAODIS(ISYM1,ISYM2) 413 KOFF2 = KQTRP + IAODIS(ISYM2,ISYM1) 414 CALL TRSREC(NBAS(ISYM1),NBAS(ISYM2), 415 & WORK(KOFF1),WORK(KOFF2)) 416 END DO 417 CALL DCOPY(N2BST(ISYMQ),WORK(KQTRP),1,WORK(KQMATH),1) 418 CALL DSCAL(N2BST(ISYMQ),HALF,WORK(KQMATH),1) 419 420 DO ISYM1 = 1, NSYM 421 ISYM2 = MULD2H(ISYM1,ISYMQ) 422 KOFF1 = KQMATP + IAODIS(ISYM1,ISYM2) 423 KOFF2 = KQTRP + IAODIS(ISYM2,ISYM1) 424 CALL TRSREC(NBAS(ISYM1),NBAS(ISYM2), 425 & WORK(KOFF1),WORK(KOFF2)) 426 END DO 427 CALL DCOPY(N2BST(ISYMQ),WORK(KQTRP),1,WORK(KQMATP),1) 428 CALL DSCAL(N2BST(ISYMQ),HALF,WORK(KQMATP),1) 429 430 XCONS(IVEC,IDXXIM) = 431 & DDOT(N2BST(ISYMQ),WORK(KQMATH),1,WORK(KXIM),1)+ 432 & DBLE(IREAL)*DDOT(N2BST(ISYMQ),WORK(KQMATP),1,WORK(KXIM),1) 433 434 435 IVEC = IVEC + 1 436 END DO 437 438 END DO 439 440*---------------------------------------------------------------------* 441* print the results: 442*---------------------------------------------------------------------* 443 IF (LOCDBG) THEN 444 WRITE(LUPRI,*) 'CCXIMCON> results for '// 445 & 'X intermediate contribs.:' 446 IF (MXVEC.NE.0) THEN 447 DO IDXXIM = 1, NXIM 448 WRITE (LUPRI,*) 'IDXXIM = ',IDXXIM 449 IVEC = 1 450 DO WHILE(IXDOTS(IVEC,IDXXIM).NE.0 .AND. IVEC.LE.MXVEC) 451 WRITE(LUPRI,'(A,2I5,2X,E19.12)') 'CCXIMCON> ', 452 & IXIM(IDXXIM),IXDOTS(IVEC,IDXXIM),XCONS(IVEC,IDXXIM) 453 IVEC = IVEC + 1 454 END DO 455 END DO 456 ELSE 457 WRITE (LUPRI,*) 'MXVEC.EQ.0 --> nothing calculated.' 458 END IF 459 END IF 460 461 RETURN 462 END 463*======================================================================* 464