1! 2! Dalton, a molecular electronic structure program 3! Copyright (C) by the authors of Dalton. 4! 5! This program is free software; you can redistribute it and/or 6! modify it under the terms of the GNU Lesser General Public 7! License version 2.1 as published by the Free Software Foundation. 8! 9! This program is distributed in the hope that it will be useful, 10! but WITHOUT ANY WARRANTY; without even the implied warranty of 11! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 12! Lesser General Public License for more details. 13! 14! If a copy of the GNU LGPL v2.1 was not distributed with this 15! code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html. 16! 17! 18C 19c /* deck cclr_setup */ 20*=====================================================================* 21 SUBROUTINE CCLR_SETUP(MXTRAN, MXVEC, 22 & IFTRAN, IFDOTS, FCONS, NFTRAN, 23 & IJTRAN, IJDOTS, CJCON, NJTRAN, 24 & IXITRAN, IXIDOTS, XICONS, NXITRAN, 25 & IRTRAN, IRDOTS, RCONS, NRTRAN, 26 & IXETRAN,IXDOTS,IEDOTS,XCONS,ECONS,NXETRAN, 27 & RESULT, MXSOP, LADD, WORK, LWORK ) 28*---------------------------------------------------------------------* 29* 30* Purpose: set up for CC linear response section: 31* - list of F matrix transformations with Cauchy vectors 32* - list of XKSI and ETA vector calculations 33* - list of X intermediate contributions 34* - list of second-order reortho./relax. contributions 35* 36* Written by Christof Haettig, may 1999 based on CCCM_SETUP 37* 38*=====================================================================* 39 USE PELIB_INTERFACE, ONLY: USE_PELIB 40#if defined (IMPLICIT_NONE) 41 IMPLICIT NONE 42#else 43# include "implicit.h" 44#endif 45#include "priunit.h" 46#include "ccorb.h" 47#include "cclrinf.h" 48#include "ccroper.h" 49#include "ccr1rsp.h" 50#include "ccsdinp.h" 51#include "ccexpfck.h" 52#include "cclists.h" 53#include "ccsections.h" 54#include "ccslvinf.h" 55 56* local parameters: 57 CHARACTER*(20) MSGDBG 58 PARAMETER (MSGDBG = '[debug] CCLR_SETUP> ') 59 LOGICAL LOCDBG 60 PARAMETER (LOCDBG = .FALSE.) 61 62 LOGICAL LADD 63 INTEGER MXVEC, MXTRAN, MXSOP 64 65 INTEGER IFTRAN(MXDIM_FTRAN,MXTRAN) 66 INTEGER IFDOTS(MXVEC,MXTRAN) 67 INTEGER IJTRAN(MXDIM_JTRAN,MXTRAN) 68 INTEGER IJDOTS(MXVEC,MXTRAN) 69 INTEGER IXITRAN(1,MXTRAN) 70 INTEGER IXIDOTS(MXVEC,MXTRAN) 71 INTEGER IRTRAN(1,MXTRAN) 72 INTEGER IRDOTS(MXVEC,MXTRAN) 73 INTEGER IXETRAN(MXDIM_XEVEC,MXTRAN) 74 INTEGER IXDOTS(MXVEC,MXTRAN), IEDOTS(MXVEC,MXTRAN) 75 76 INTEGER NFTRAN, NXITRAN, NRTRAN, NXETRAN, LWORK 77 INTEGER NJTRAN 78 79#if defined (SYS_CRAY) 80 REAL RESULT(MXSOP) 81 REAL FCONS(MXVEC,MXTRAN) 82 REAL CJCON(MXVEC,MXTRAN) 83 REAL XICONS(MXVEC,MXTRAN) 84 REAL RCONS(MXVEC,MXTRAN) 85 REAL XCONS(MXVEC,MXTRAN), ECONS(MXVEC,MXTRAN) 86 REAL WORK(LWORK) 87 REAL ZERO, SIGN 88 REAL WSTAT, WNUCL, WREO, WXE1, WXE2, WXI1, WXI2, WF, WJ 89#else 90 DOUBLE PRECISION RESULT(MXSOP) 91 DOUBLE PRECISION FCONS(MXVEC,MXTRAN) 92 DOUBLE PRECISION CJCON(MXVEC,MXTRAN) 93 DOUBLE PRECISION XICONS(MXVEC,MXTRAN) 94 DOUBLE PRECISION RCONS(MXVEC,MXTRAN) 95 DOUBLE PRECISION XCONS(MXVEC,MXTRAN), ECONS(MXVEC,MXTRAN) 96 DOUBLE PRECISION WORK(LWORK) 97 DOUBLE PRECISION ZERO, SIGN 98 DOUBLE PRECISION WSTAT, WNUCL, WREO, WXE1, WXE2, WXI1, WXI2, WF,WJ 99#endif 100 PARAMETER (ZERO = 0.0D0) 101 102 LOGICAL LORXA, LORXB, LPDBSA, LPDBSB 103 INTEGER ISYMA, ISYMB, ITRAN, IVEC, ISYML, IDUM, I, N, IFREQ 104 INTEGER IR1VECA,IR1VECB,IOPERA,IOPERB,IEATA1A,IEATA1B, ISGNSOP 105 INTEGER IL1VECB, INUM, IOPER, NBSOP, ISIGN, ISYOP, IKAPA, IKAPB 106 INTEGER MEAVEC,MFVEC,MXAVEC,MXIVEC,MXRVEC,MXEVEC, IDX, IEXPV 107 INTEGER MJVEC, IL1VECA 108 109 CHARACTER LABELA*(8), LABELB*(8), LABSOP*(8) 110 111* external functions: 112 INTEGER IR1TAMP 113 INTEGER IR1KAPPA 114 INTEGER IL1ZETA 115 INTEGER IETA1 116 INTEGER IEXPECT 117#if defined (SYS_CRAY) 118 REAL CC_NUCCON 119#else 120 DOUBLE PRECISION CC_NUCCON 121#endif 122 123*---------------------------------------------------------------------* 124* initializations: 125*---------------------------------------------------------------------* 126 DO ITRAN = 1, MXTRAN 127 DO I = 1, MXDIM_XEVEC 128 IXETRAN(I,ITRAN) = 0 129 END DO 130 DO I = 1, MXDIM_FTRAN 131 IFTRAN(I,ITRAN) = 0 132 END DO 133 DO I = 1, MXDIM_JTRAN 134 IJTRAN(I,ITRAN) = 0 135 END DO 136 DO I = 1, 1 137 IXITRAN(I,ITRAN) = 0 138 IRTRAN(I,ITRAN) = 0 139 END DO 140 141 DO IVEC = 1, MXVEC 142 IFDOTS(IVEC,ITRAN) = 0 143 IJDOTS(IVEC,ITRAN) = 0 144 IXIDOTS(IVEC,ITRAN) = 0 145 IRDOTS(IVEC,ITRAN) = 0 146 IXDOTS(IVEC,ITRAN) = 0 147 IEDOTS(IVEC,ITRAN) = 0 148 END DO 149 END DO 150 151 NFTRAN = 0 152 NJTRAN = 0 153 NXITRAN = 0 154 NRTRAN = 0 155 NXETRAN = 0 156 157 MFVEC = 0 158 MJVEC = 0 159 MXIVEC = 0 160 MXRVEC = 0 161 MXEVEC = 0 162 163 NBSOP = 0 164 165*---------------------------------------------------------------------* 166* start loop over all requested second-order properties: 167*---------------------------------------------------------------------* 168 DO IOPER = 1, NLROP 169 IOPERA = IALROP(IOPER) 170 IOPERB = IBLROP(IOPER) 171 LORXA = LALORX(IOPER) 172 LORXB = LBLORX(IOPER) 173 ISYMA = ISYOPR(IOPERA) 174 ISYMB = ISYOPR(IOPERB) 175 LABELA = LBLOPR(IOPERA) 176 LABELB = LBLOPR(IOPERB) 177 LPDBSA = LPDBSOP(IOPERA) 178 LPDBSB = LPDBSOP(IOPERB) 179 180 181 IF (ISYMA.EQ.ISYMB) THEN 182 183 DO IFREQ = 1, NBLRFR 184 NBSOP = NBSOP + 1 185 186 IF (NBSOP.GT.MXSOP) THEN 187 CALL QUIT('NBSOP out of range in CCLR_SETUP.') 188 END IF 189 190 DO ISIGN = 1, -1, -2 191 SIGN = DBLE(ISIGN) 192 193*---------------------------------------------------------------------* 194* in all cases we need Eta{A} x R1^B 195*---------------------------------------------------------------------* 196 IR1VECB = IR1TAMP(LABELB,LORXB,SIGN*BLRFR(IFREQ),IDUM) 197C IEATA1A = IETA1(LABELA,LORXA,SIGN*ALRFR(IFREQ),IDUM) 198 IF (LORXA) THEN 199 IKAPA = IR1KAPPA(LABELA,SIGN*ALRFR(IFREQ),IDUM) 200 ELSE 201 IKAPA = 0 202 END IF 203 204 CALL CC_SETXE('Eta',IXETRAN,IEDOTS,MXTRAN,MXVEC, 205 & 0,IOPERA,IKAPA,0,0,0,IR1VECB,ITRAN,IVEC) 206 NXETRAN = MAX(NXETRAN,ITRAN) 207 MXEVEC = MAX(MXEVEC, IVEC) 208 WXE1 = ECONS(IVEC,ITRAN) 209 210 IF (.NOT. ASYMSD) THEN 211*---------------------------------------------------------------------* 212* symmetric approach: add F * R1^A * R1^B + Eta{B} x R1^A 213*---------------------------------------------------------------------* 214 IR1VECA = IR1TAMP(LABELA,LORXA,SIGN*ALRFR(IFREQ),IDUM) 215C IEATA1B = IETA1(LABELB,LORXB,SIGN*BLRFR(IFREQ),IDUM) 216 IF (LORXB) THEN 217 IKAPB = IR1KAPPA(LABELB,SIGN*BLRFR(IFREQ),IDUM) 218 ELSE 219 IKAPB = 0 220 END IF 221 222 IF (CIS) THEN 223 WF = ZERO 224 ELSE 225 CALL CCQR_SETF(IFTRAN,IFDOTS,MXTRAN,MXVEC, 226 & 0,IR1VECA,IR1VECB,ITRAN,IVEC) 227 NFTRAN = MAX(NFTRAN,ITRAN) 228 MFVEC = MAX(MFVEC, IVEC) 229 WF = FCONS(IVEC,ITRAN) 230C 231 IF (CCSLV.OR.USE_PELIB()) THEN 232 IL1VECA = IL1ZETA(LABELA,LORXA,SIGN*ALRFR(IFREQ),IDUM) 233 IL1VECB = IL1ZETA(LABELB,LORXB,SIGN*BLRFR(IFREQ),IDUM) 234C Here, we substitute amplitudes <-> multipilers 235 CALL CCQR_SETF(IJTRAN,IJDOTS,MXTRAN,MXVEC, 236 & 0,IL1VECA,IL1VECB,ITRAN,IVEC) 237 NJTRAN = MAX(NJTRAN,ITRAN) 238 MJVEC = MAX(MJVEC, IVEC) 239 WJ = CJCON(IVEC,ITRAN) 240 ELSE 241 WJ = ZERO 242 END IF 243 END IF 244 245 CALL CC_SETXE('Eta',IXETRAN,IEDOTS,MXTRAN,MXVEC, 246 & 0,IOPERB,IKAPB,0,0,0,IR1VECA,ITRAN,IVEC) 247 NXETRAN = MAX(NXETRAN,ITRAN) 248 MXEVEC = MAX(MXEVEC, IVEC) 249 WXE2 = ECONS(IVEC,ITRAN) 250 251 ELSE 252*---------------------------------------------------------------------* 253* asymmetric approach: add L1^B x Xksi{A} 254*---------------------------------------------------------------------* 255 WF = ZERO 256 WJ = ZERO 257 258 IL1VECB = IL1ZETA(LABELB,LORXB,SIGN*BLRFR(IFREQ),IDUM) 259 260 CALL CC_SETXE('Xi ',IXETRAN,IXDOTS,MXTRAN,MXVEC, 261 & 0,IOPERA,IKAPA,0,0,0,IL1VECB,ITRAN,IVEC) 262 NXETRAN = MAX(NXETRAN,ITRAN) 263 MXEVEC = MAX(MXEVEC, IVEC) 264 WXE2 = XCONS(IVEC,ITRAN) 265 266 END IF 267 268*---------------------------------------------------------------------* 269* for orbital relaxed second-order properties or if we have 270* perturbation-dependent basis sets involved add the contrib. 271* from the first-order effective Fock matrix times the 272* Q matrix (kappa + R) : 273*---------------------------------------------------------------------* 274 WXI1 = ZERO 275 WXI2 = ZERO 276 277 IF (LORXB.OR.LPDBSB) THEN 278 279 IKAPA = IR1KAPPA(LABELA,SIGN*ALRFR(IFREQ),IDUM) 280 IKAPB = IR1KAPPA(LABELB,SIGN*BLRFR(IFREQ),IDUM) 281 282 CALL CC_SETDOT(IXITRAN,IXIDOTS,MXTRAN,MXVEC, 283 & IKAPA,IKAPB,ITRAN,IVEC) 284 NXITRAN = MAX(NXITRAN,ITRAN) 285 MXIVEC = MAX(MXIVEC, IVEC) 286 WXI1 = WXI1 + XICONS(IVEC,ITRAN) 287 288 289 END IF 290 291 292 IF (LORXA.OR.LPDBSA) THEN 293 294 IKAPA = IR1KAPPA(LABELA,SIGN*ALRFR(IFREQ),IDUM) 295 IKAPB = IR1KAPPA(LABELB,SIGN*BLRFR(IFREQ),IDUM) 296 297 CALL CC_SETDOT(IXITRAN,IXIDOTS,MXTRAN,MXVEC, 298 & IKAPB,IKAPA,ITRAN,IVEC) 299 NXITRAN = MAX(NXITRAN,ITRAN) 300 MXIVEC = MAX(MXIVEC, IVEC) 301 WXI2 = WXI2 + XICONS(IVEC,ITRAN) 302 303 END IF 304 305*---------------------------------------------------------------------* 306* for derivatives we might need to include coupling between 307* relaxation and reorthonormalization: 308*---------------------------------------------------------------------* 309 IF ( (LPDBSA .AND. LORXB) .OR. (LPDBSB .AND. LORXA) ) THEN 310 311 IKAPA = IR1KAPPA(LABELA,SIGN*ALRFR(IFREQ),IDUM) 312 IKAPB = IR1KAPPA(LABELB,SIGN*BLRFR(IFREQ),IDUM) 313 314 CALL CC_SETDOT(IRTRAN,IRDOTS,MXTRAN,MXVEC, 315 & IKAPB,IKAPA,ITRAN,IVEC) 316 NRTRAN = MAX(NRTRAN,ITRAN) 317 MXRVEC = MAX(MXRVEC,IVEC) 318 WREO = RCONS(IVEC,ITRAN) 319 320 ELSE 321 WREO = ZERO 322 END IF 323 324*---------------------------------------------------------------------* 325* get "static" and nuclear contribution: 326*---------------------------------------------------------------------* 327 IF (LPDBSA .OR. LPDBSB) THEN 328 329 CALL CC_FIND_SO_OP(LABELA,LABELB,LABSOP,ISYOP,ISGNSOP, 330 & INUM,WORK,LWORK) 331 332 IF (INUM.LT.0) CALL QUIT('Operator error in CCLR_SETUP.') 333 334 IEXPV = IEXPECT(LABSOP,ISYOP,1) 335 WSTAT = EXPVALUE(1,IEXPV) + EXPVALUE(2,IEXPV) 336 337 WNUCL = CC_NUCCON(LABSOP,ISYOP) 338 339 ELSE 340 WSTAT = ZERO 341 WNUCL = ZERO 342 END IF 343 344*---------------------------------------------------------------------* 345* add contributions together: 346*---------------------------------------------------------------------* 347 IF (LADD) THEN 348 349 IDX = NBLRFR*(IOPER-1) + IFREQ 350 IF (ISIGN.EQ.-1) IDX = IDX + NBLRFR*NLROP 351 352 RESULT(IDX) = WXE1+WXE2+WF-WJ+WXI1+WXI2+WREO-WSTAT-WNUCL 353 354 IF (LOCDBG) THEN 355 WRITE (LUPRI,*) 'LABELA,LABELB:',LABELA,LABELB 356 WRITE (LUPRI,*) 'IFREQ:',IFREQ 357 WRITE (LUPRI,*) 'IDX = ',IDX 358 WRITE (LUPRI,*) 'WSTAT,WNUCL:',WSTAT,WNUCL 359 WRITE (LUPRI,*) 'WXI1,WXI2,WREO:', WXI1, WXI2, WREO 360 WRITE (LUPRI,*) 'WXE1, WXE2, WF, WJ :', 361 & WXE1, WXE2, WF, WJ 362 WRITE (LUPRI,*) 'RESULT:',RESULT(IDX) 363 END IF 364 365 END IF 366 367*---------------------------------------------------------------------* 368* end loop over second-order properties 369*---------------------------------------------------------------------* 370 END DO 371 END DO 372 373 END IF 374 END DO 375 376 IF (MFVEC.GT.MXVEC) THEN 377 CALL QUIT('MFVEC has been out of bounds in CCLR_SETUP.') 378 ELSE IF (MXEVEC.GT.MXVEC) THEN 379 CALL QUIT('MXEVEC has been out of bounds in CCLR_SETUP.') 380 ELSE IF (MXIVEC.GT.MXVEC) THEN 381 CALL QUIT('MXIVEC has been out of bounds in CCLR_SETUP.') 382 ELSE IF (MXRVEC.GT.MXVEC) THEN 383 CALL QUIT('MXRVEC has been out of bounds in CCLR_SETUP.') 384 ELSE IF (MJVEC.GT.MXVEC) THEN 385 CALL QUIT('MJVEC has been out of bounds in CCLR_SETUP.') 386 END IF 387 388*---------------------------------------------------------------------* 389* print the lists: 390*---------------------------------------------------------------------* 391* general statistics: 392 IF ((.NOT.LADD) .OR. LOCDBG) THEN 393 WRITE(LUPRI,'(/,/3X,A,I3,A)') 'For the requested',NBSOP, 394 & ' second-order properties' 395 WRITE(LUPRI,'((8X,A,I3,A))') 396 & ' - ',NFTRAN,' F matrix transformations with R1 vectors', 397 & ' - ',NJTRAN,' J matrix transformations with L1 vectors', 398 & ' - ',NXETRAN,' ETA and XKSI vector calculations ', 399 & ' - ',NXITRAN,' X intermediate calculations ', 400 & ' - ',NRTRAN, ' 2. order reortho./relax. contributions' 401 WRITE(LUPRI,'(3X,A,/,/)') 'will be performed.' 402 END IF 403 404 IF (LOCDBG) THEN 405 406 ! F matrix transformations: 407 WRITE(LUPRI,*)'List of F matrix transformations:' 408 DO ITRAN = 1, NFTRAN 409 WRITE(LUPRI,'(A,2I5,5X,(25I3,20X))') MSGDBG, 410 & (IFTRAN(I,ITRAN),I=1,2),(IFDOTS(I,ITRAN),I=1,MFVEC) 411 END DO 412 WRITE(LUPRI,*) 413 414 ! J matrix transformations: 415 WRITE(LUPRI,*)'List of J matrix transformations:' 416 DO ITRAN = 1, NJTRAN 417 WRITE(LUPRI,'(A,2I5,5X,(25I3,20X))') MSGDBG, 418 & (IJTRAN(I,ITRAN),I=1,2),(IJDOTS(I,ITRAN),I=1,MJVEC) 419 END DO 420 WRITE(LUPRI,*) 421 422 ! Xi{O} and ETA{O} vector calculations: 423 WRITE(LUPRI,*) 'List of Xi{O} and ETA{O} vector calculations:' 424 DO ITRAN = 1, NXETRAN 425 WRITE(LUPRI,'(A,5I5,5X,(25I3,20X))') MSGDBG, 426 & (IXETRAN(I,ITRAN),I=1,5),(IXDOTS(I,ITRAN),I=1,MXEVEC) 427 WRITE(LUPRI,'(A,25X,5X,(25I3,20X))') MSGDBG, 428 & (IEDOTS(I,ITRAN),I=1,MXEVEC) 429 END DO 430 WRITE(LUPRI,*) 431 432 ! X{O} intermediate calculations: 433 WRITE(LUPRI,*) 'List of X{O} intermediate calculations:' 434 DO ITRAN = 1, NXITRAN 435 WRITE(LUPRI,'(A,2I5,5X,(25I3,20X))') MSGDBG, 436 & (IXITRAN(I,ITRAN),I=1,1),(IXIDOTS(I,ITRAN),I=1,MXIVEC) 437 END DO 438 WRITE(LUPRI,*) 439 440 END IF 441 442 RETURN 443 END 444 445*---------------------------------------------------------------------* 446* END OF SUBROUTINE CCLR_SETUP * 447*---------------------------------------------------------------------* 448