1*=====================================================================* 2c /* deck ccxopa_eomsetup */ 3*=====================================================================* 4 SUBROUTINE CCXOPA_EOMSETUP(IFTRAN, IFDOTS, FCONS, 5 & NFTRAN, MXFTRAN, MXFVEC, 6 & IEATRAN, IEADOTS, EACONS, 7 & NXE1TRAN,MXATRAN, MXAVEC, 8 & IXE2TRAN,IX2DOTS, X2CONS, 9 & NXE2TRAN,MXXTRAN, MXXVEC, 10 & IEOMXE2TRAN,IEOMX2DOTS, EOMX2CONS, 11 & NEOMXE2TRAN,MXEOMXTRAN, MXEOMXVEC, 12 & IEOML0TRAN,IEOML0DOTS,EOML0CONS, 13 & NEOML0TRAN,MXEOML0TRAN,MXEOML0VEC, 14 & RESULT, MXOPA, LADD, WORK, LWORK ) 15*---------------------------------------------------------------------* 16* 17* Purpose: set up for CC first-order transition moments 18* - list of B matrix transformations with eigenvectors 19* - list of A{X} matrix transformations with eigenvectors 20* - list of XKSI vector contractions with Nbar multipliers 21* 22* Written by Christof Haettig, Oct 2003 23* 24* For EOM XOPA 25* - list of XKSI vector contractions with LE vectors 26* and of dot products tbar*RE) 27* Sonia Coriani, Nov 2015 28* 29*=====================================================================* 30 IMPLICIT NONE 31#include "priunit.h" 32#include "cclists.h" 33#include "ccxopainf.h" 34#include "ccroper.h" 35#include "ccexci.h" 36#include "ccsdinp.h" 37#include "ccorb.h" 38 39* local parameters: 40 CHARACTER*(22) MSGDBG 41 PARAMETER (MSGDBG = '[debug] CCXOPA_EOMSETUP> ') 42 LOGICAL LOCDBG 43 PARAMETER (LOCDBG = .false.) 44 45 LOGICAL LADD 46 INTEGER MXOPA,MXFTRAN,MXFVEC,MXATRAN,MXAVEC,MXXTRAN,MXXVEC 47 INTEGER MXEOMXTRAN,MXEOMXVEC 48 INTEGER MXEOML0TRAN,MXEOML0VEC 49 50 INTEGER IFTRAN(MXDIM_FTRAN,MXFTRAN) 51 INTEGER IFDOTS(MXFVEC,MXFTRAN) 52 INTEGER IEATRAN(MXDIM_XEVEC,MXATRAN) 53 INTEGER IEADOTS(MXAVEC,MXATRAN) 54 INTEGER IXE2TRAN(MXDIM_XEVEC,MXXTRAN) 55 INTEGER IX2DOTS(MXXVEC,MXXTRAN) 56! 57! EOM is one of these 58! 59 INTEGER IEOMXE2TRAN(MXDIM_XEVEC,MXEOMXTRAN) 60 INTEGER IEOMX2DOTS(MXEOMXVEC,MXEOMXTRAN) 61 INTEGER IEOML0TRAN(MXEOML0TRAN) 62 INTEGER IEOML0DOTS(MXEOML0VEC,MXEOML0VEC) 63 64 INTEGER NFTRAN, NXE1TRAN, NXE2TRAN, LWORK 65 INTEGER NEOMXE2TRAN, NEOML0TRAN 66 67#if defined (SYS_CRAY) 68 REAL RESULT(MXOPA) 69 REAL FCONS(MXFVEC,MXFTRAN) 70 REAL EACONS(MXAVEC,MXATRAN) 71 REAL X2CONS(MXXVEC,MXXTRAN) 72 REAL EOMX2CONS(MXEOMXVEC,MXEOMXTRAN) !(csi*LE) 73 REAL EOML0CONS(MXEOML0VEC,MXEOML0VEC) !(Tbar*RE) 74 !REAL EOMCONS(MXXVEC,MXEOMXTRAN) !(csi*LE).(Tbar*RE) 75 REAL WORK(LWORK) 76 REAL ZERO, SIGN, EIGVI, EIGVF 77 REAL WIAF, WXINIF, WIBF, WLiXi, WL0Rf, WLiXiRf 78#else 79 DOUBLE PRECISION RESULT(MXOPA) 80 DOUBLE PRECISION FCONS(MXFVEC,MXFTRAN) 81 DOUBLE PRECISION EACONS(MXAVEC,MXATRAN) 82 DOUBLE PRECISION X2CONS(MXXVEC,MXXTRAN) 83 DOUBLE PRECISION EOMX2CONS(MXEOMXVEC,MXEOMXTRAN)!(csi*LE) 84 DOUBLE PRECISION EOML0CONS(MXEOML0VEC,MXEOML0TRAN) !(Tbar*RE) 85 !DOUBLE PRECISION EOMCONS(MXEOMXVEC,MXEOMXTRAN) !(csi*LE).(Tbar*RE) 86 DOUBLE PRECISION WORK(LWORK) 87 DOUBLE PRECISION ZERO, SIGN, EIGVI, EIGVF 88 DOUBLE PRECISION WIAF, WXINIF, WIBF, WLiXi, WL0Rf, WLiXiRf 89#endif 90 PARAMETER (ZERO = 0.0D0) 91 92 CHARACTER LABEL*(8) 93 LOGICAL LORX, LPDBS 94 INTEGER ITRAN, I, IRSD, IRSDX, ISTATEI, ISTATEF, ISYMI, ISYMF, 95 & ISTISY, ISTFSY, IOP, IOPER, ISYMO, ISYME, ITURN, 96 & IKAP, MXEAVEC, MXE2VEC, IN2VEC, IR1VEC, MFVEC, 97 & ITMIF, IVEC, NBOPA, IDUM, 98 & IMULI, IMULF 99 100 INTEGER MXEOME2VEC, MXEOM2 101 INTEGER NL0TRAN 102 103* external functions: 104 INTEGER IR1TAMP 105 INTEGER IN2AMP 106 107*---------------------------------------------------------------------* 108* initializations: 109* initialize for EOM as well .... 110*---------------------------------------------------------------------* 111 DO ITRAN = 1, MXATRAN 112 IEATRAN(1,ITRAN) = 0 113 IEATRAN(2,ITRAN) = 0 114 IEATRAN(3,ITRAN) = -1 115 IEATRAN(4,ITRAN) = -1 116 IEATRAN(5,ITRAN) = 0 117 DO IVEC = 1, MXAVEC 118 IEADOTS(IVEC,ITRAN) = 0 119 END DO 120 END DO 121 122 DO ITRAN = 1, MXXTRAN 123 IXE2TRAN(1,ITRAN) = 0 124 IXE2TRAN(2,ITRAN) = 0 125 IXE2TRAN(3,ITRAN) = -1 126 IXE2TRAN(4,ITRAN) = -1 127 IXE2TRAN(5,ITRAN) = 0 128 DO IVEC = 1, MXXVEC 129 IX2DOTS(IVEC,ITRAN) = 0 130 END DO 131 END DO 132 133 DO ITRAN = 1, MXEOMXTRAN 134 IEOMXE2TRAN(1,ITRAN) = 0 135 IEOMXE2TRAN(2,ITRAN) = 0 136 IEOMXE2TRAN(3,ITRAN) = -1 137 IEOMXE2TRAN(4,ITRAN) = -1 138 IEOMXE2TRAN(5,ITRAN) = 0 139 DO IVEC = 1, MXEOMXVEC 140 IEOMX2DOTS(IVEC,ITRAN) = 0 141 END DO 142 END DO 143 144 !megaredundant but I am too tired to make it smarter... 145 DO ITRAN = 1, MXEOML0TRAN 146 IEOML0TRAN(ITRAN) = 0 147 DO IVEC = 1, MXEOML0VEC 148 IEOML0DOTS(IVEC,ITRAN) = 0 149 END DO 150 END DO 151 152 DO ITRAN = 1, MXFTRAN 153 DO I = 1, 3 154 IFTRAN(I,ITRAN) = 0 155 END DO 156 DO IVEC = 1, MXFVEC 157 IFDOTS(IVEC,ITRAN) = 0 158 END DO 159 END DO 160 161 NFTRAN = 0 162 NXE1TRAN = 0 163 NXE2TRAN = 0 164 NEOMXE2TRAN = 0 165 NEOML0TRAN = 0 166 167 NBOPA = 0 168 MFVEC = 0 169 MXE2VEC = 0 170 MXEAVEC = 0 171 172 MXEOME2VEC = 0 173 MXEOM2 = 0 174 NL0TRAN = 0 175 176! mi manca qualcosa qua sopra... 177 178*---------------------------------------------------------------------* 179* start loop over all requested transition moments: 180*---------------------------------------------------------------------* 181 DO IRSDX = 1, 2*NXQR2ST 182 ITURN = 1 + (IRSDX-1)/NXQR2ST 183 IRSD = IRSDX - (ITURN-1)*NXQR2ST 184 185 IF (ITURN.EQ.1) THEN 186 ISTATEI = IQR2ST(IRSD,1) 187 ISTATEF = IQR2ST(IRSD,2) 188 ELSE IF (ITURN.EQ.2) THEN 189 ! switch state indices (and thereby also the sign of the freqs) 190 ! to get the conjugated transition moments 191 ISTATEI = IQR2ST(IRSD,2) 192 ISTATEF = IQR2ST(IRSD,1) 193 ELSE 194 CALL QUIT('Error in CCXOPA_EOMSETUP') 195 END IF 196 197 ISYMI = ISYEXC(ISTATEI) 198 ISYMF = ISYEXC(ISTATEF) 199 ISYME = MULD2H(ISYMI,ISYMF) 200 ISTISY = ISTATEI - ISYOFE(ISYMI) 201 ISTFSY = ISTATEF - ISYOFE(ISYMF) 202 EIGVI = EIGVAL(ISTATEI) 203 EIGVF = EIGVAL(ISTATEF) 204CRF We need to know multiplicities as well. 205 IMULI = IMULTE(ISTATEI) 206 IMULF = IMULTE(ISTATEF) 207 208 IF (LOCDBG) THEN 209 WRITE(LUPRI,*) 'CCXOPA_EOMSETUP:' 210 WRITE(LUPRI,*) 'ITURN,IRSD:',ITURN,IRSD 211 WRITE(LUPRI,*) 'ISTATEI,ISTATEF:',ISTATEI,ISTATEF 212 WRITE(LUPRI,*) 'ISYMI,ISYMF:',ISYMI,ISYMF 213 WRITE(LUPRI,*) 'ISTISY,ISTFSY:',ISTISY,ISTFSY 214 WRITE(LUPRI,*) 'EIGVI,EIGVF:',EIGVI,EIGVF 215 WRITE(LUPRI,*) 'IMULI,IMULF:',IMULI,IMULF 216 END IF 217 218 IF (IMULI.NE.IMULF) THEN 219 WRITE(LUPRI,*) 'Only singlet operators are currently' 220 & //' implemented!' 221 CYCLE 222 END IF 223 224 DO IOP = 1, NQR2OP 225 IOPER = IQR2OP(IOP) 226 LORX = .FALSE. 227 ISYMO = ISYOPR(IOPER) 228 LABEL = LBLOPR(IOPER) 229 LPDBS = LPDBSOP(IOPER) 230 IKAP = 0 231 232 IF (LPDBS) CALL QUIT('perturbation-dependent basis sets not '// 233 & 'implemented in CCXOPA_EOMSETUP.') 234 235 IF (ISYMO.EQ.ISYME) THEN 236 237 NBOPA = NBOPA + 1 238 239 IF (NBOPA.GT.MXOPA) THEN 240 CALL QUIT('NBOPA out of range in CCXOPA_EOMSETUP.') 241 END IF 242 243*---------------------------------------------------------------------* 244* in all cases we need LE x A{X} x RE 245*---------------------------------------------------------------------* 246 !write(lupri,*) "Call CC_SETXE('Eta')" 247 !call flshfo(lupri) 248 249 CALL CC_SETXE('Eta',IEATRAN,IEADOTS,MXATRAN,MXAVEC, 250 & ISTATEI,IOPER,IKAP,0,0,0,ISTATEF,ITRAN,IVEC) 251 NXE1TRAN = MAX(NXE1TRAN,ITRAN) 252 MXEAVEC = MAX(MXEAVEC, IVEC) 253 WIAF = EACONS(IVEC,ITRAN) 254 255*---------------------------------------------------------------------* 256* add N2 * Xksi{X} or LE * B * RE * R1, depending on QR22N1 257*---------------------------------------------------------------------* 258 WXINIF = ZERO 259 WIBF = ZERO 260 WLiXi = ZERO 261 WL0Rf = ZERO 262 WLiXiRf = ZERO 263 264 IF (.NOT.CIS) THEN 265 IF (QR22N1) THEN 266 IN2VEC=IN2AMP(ISTATEI,-EIGVI,ISYMI,ISTATEF,+EIGVF,ISYMF) 267 CALL CC_SETXE('Xi ',IXE2TRAN,IX2DOTS,MXXTRAN,MXXVEC, 268 & 0,IOPER,IKAP,0,0,0,IN2VEC,ITRAN,IVEC) 269 NXE2TRAN = MAX(NXE2TRAN,ITRAN) 270 MXE2VEC = MAX(MXE2VEC, IVEC) 271 WXINIF = X2CONS(IVEC,ITRAN) 272 ELSE IF (LEOMXOPA.AND.(IMULF.EQ.1).AND.(ISYMF.EQ.1)) then 273 ! Only do this term if the F state is totally symmetric in 274 ! both spin and space! 275 ! 276 ! recover index of LE and do dot prod with Xi(X) 277 ! 278 CALL CC_SETXE('Xi ',IEOMXE2TRAN,IEOMX2DOTS, 279 & MXEOMXTRAN,MXEOMXVEC, 280 & 0,IOPER,IKAP,0,0,0,ISTATEI,ITRAN,IVEC) 281 NEOMXE2TRAN = MAX(NEOMXE2TRAN,ITRAN) 282 MXEOME2VEC = MAX(MXEOME2VEC, IVEC) 283 WLiXi = EOMX2CONS(IVEC,ITRAN) 284 !write(lupri,*) "WLiXi:", WLiXi 285 ! 286 !now we fill in the list required for the dot products 287 ! 288 CALL CC_SETDOT(IEOML0TRAN,IEOML0DOTS, 289 & MXEOML0TRAN,MXEOML0VEC, 290 & 0,ISTATEF,ITRAN,IVEC) 291 NEOML0TRAN = MAX(NEOML0TRAN,ITRAN) 292 MXEOM2 = MAX(MXEOM2,IVEC) 293 294 WL0Rf = EOML0CONS(IVEC,ITRAN) 295 !write(lupri,*) "WL0Rf:", WL0Rf 296 WLiXiRf = WLiXi*WL0Rf 297 ELSE IF(.NOT.LEOMXOPA) THEN 298 IR1VEC = IR1TAMP(LABEL,LORX,EIGVI-EIGVF,IDUM) 299 CALL CC_SETF12(IFTRAN,IFDOTS,MXFTRAN,MXFVEC, 300 & ISTATEI,ISTATEF,IR1VEC,ITRAN,IVEC) 301 NFTRAN = MAX(NFTRAN,ITRAN) 302 MFVEC = MAX(MFVEC, IVEC) 303 WIBF = FCONS(IVEC,ITRAN) 304 END IF 305 !end if 306 END IF 307 308*---------------------------------------------------------------------* 309* add contributions together: 310*---------------------------------------------------------------------* 311 IF (LADD) THEN 312 313 ITMIF = (NQR2OP*(IRSD-1) + IOP-1)*2 + ITURN 314 315 RESULT(ITMIF) = WIAF + WXINIF + WIBF - WLIXiRF 316 317 IF (LOCDBG) THEN 318 WRITE (LUPRI,*) '----- Summary after add ------' 319 WRITE (LUPRI,*) 'IDX of result = ',ITMIF 320 WRITE (LUPRI,*) 'ISTATEI, EIGVI:',ISTATEI,EIGVI 321 WRITE (LUPRI,*) 'ISTATEF, EIGVF:',ISTATEF,EIGVF 322 WRITE (LUPRI,*) 'OPERATOR:',LABEL 323 WRITE (LUPRI,*) 'L^i A{X} x R^f :',WIAF 324 WRITE (LUPRI,*) 'N^if x Xksi{X}:',WXINIF 325 WRITE (LUPRI,*) 'L^i x B x R^f x R^X:',WIBF 326 WRITE (LUPRI,*) 'L^i x Xksi{X}:',WLiXi 327 WRITE (LUPRI,*) 'R^f x L0:',WL0Rf 328 WRITE (LUPRI,*) '(L^i x Xksi{X})(R^f x L0):',WLiXiRF 329 WRITE (LUPRI,*) 'Total result:',RESULT(ITMIF) 330 WRITE (LUPRI,*) '------------------------------' 331 END IF 332 333 END IF 334 335*---------------------------------------------------------------------* 336* end loop over transition moments 337*---------------------------------------------------------------------* 338 339 END IF 340 END DO 341 END DO 342 343 IF (MFVEC.GT.MXFVEC) THEN 344 CALL QUIT('MFVEC has been out of bounds CCXOPA_EOMSETUP.') 345 ELSE IF (MXEAVEC.GT.MXAVEC) THEN 346 CALL QUIT('MXEAVEC has been out of bounds CCXOPA_EOMSETUP.') 347 ELSE IF (MXE2VEC.GT.MXXVEC) THEN 348 CALL QUIT('MXE2VEC has been out of bounds CCXOPA_EOMSETUP.') 349 ELSE IF (MXEOME2VEC.GT.MXEOMXVEC) THEN 350 CALL QUIT('MXEOME2VEC has been out of bounds CCXOPA_EOMSETUP.') 351 ELSE IF (MXEOM2.GT.MXEOML0VEC) THEN 352 CALL QUIT('MXEOM2 has been out of bounds CCXOPA_EOMSETUP.') 353 ELSE IF (NFTRAN.GT.MXFTRAN) THEN 354 CALL QUIT('NFTRAN has been out of bounds CCXOPA_EOMSETUP.') 355 ELSE IF (NXE1TRAN.GT.MXATRAN) THEN 356 CALL QUIT('NXE1TRAN has been out of bounds CCXOPA_EOMSETUP.') 357 ELSE IF (NXE2TRAN.GT.MXXTRAN) THEN 358 CALL QUIT('NXE2TRAN has been out of bounds CCXOPA_EOMSETUP.') 359 ELSE IF (NEOMXE2TRAN.GT.MXEOMXTRAN) THEN 360 CALL QUIT('NEOMXE2TRAN has been out bounds CCXOPA_EOMSETUP.') 361 ELSE IF (NEOML0TRAN.GT.MXEOML0TRAN) THEN 362 CALL QUIT('NEOML0TRAN has been out bounds CCXOPA_EOMSETUP.') 363 END IF 364 365*---------------------------------------------------------------------* 366* print the lists: 367*---------------------------------------------------------------------* 368* general statistics: 369 IF ((.NOT.LADD) .OR. LOCDBG) THEN 370 WRITE(LUPRI,'(/,/3X,A,I3,A)') 'For the requested',NBOPA, 371 & ' transition moments' 372 WRITE(LUPRI,'((8X,A,I3,A))') 373 & ' - ',NFTRAN, ' F matrix transformations with RE vectors', 374 & ' - ',NXE1TRAN,' A{X} matrix transformations with LE vectors', 375 & ' - ',NXE2TRAN,' extra XKSI vector calculations ', 376 & ' - ',NEOMXE2TRAN,' extra XKSI vector calculations (EOM)', 377 & ' - ',NEOML0TRAN,' extra L0 vector (EOM)' 378 WRITE(LUPRI,'(3X,A,/,/)') 'will be performed.' 379 END IF 380 381 IF (LOCDBG) THEN 382 383 ! F matrix transformations: 384 WRITE(LUPRI,*)'List of F matrix transformations:' 385 DO ITRAN = 1, NFTRAN 386 WRITE(LUPRI,'(A,2I5,5X,(25I3,20X))') MSGDBG, 387 & (IFTRAN(I,ITRAN),I=1,2),(IFDOTS(I,ITRAN),I=1,MFVEC) 388 END DO 389 WRITE(LUPRI,*) 390 391 ! LE x A{X} vector calculations: 392 WRITE(LUPRI,*) 'List of A{O} matrix transformations:' 393 DO ITRAN = 1, NXE1TRAN 394 WRITE(LUPRI,'(A,5I5,5X,(25I3,20X))') MSGDBG, 395 & (IEATRAN(I,ITRAN),I=1,5),(IEADOTS(I,ITRAN),I=1,MXEAVEC) 396 END DO 397 WRITE(LUPRI,*) 398 399 ! extra Xi{O} vector calculations: 400 WRITE(LUPRI,*) 'List of extra Xi{O} vector calculations:' 401 DO ITRAN = 1, NXE2TRAN 402 WRITE(LUPRI,'(A,5I5,5X,(25I3,20X))') MSGDBG, 403 & (IXE2TRAN(I,ITRAN),I=1,5),(IX2DOTS(I,ITRAN),I=1,MXE2VEC) 404 END DO 405 WRITE(LUPRI,*) 406 407 ! extra Xi{O} vector calculations for EOM: 408 WRITE(LUPRI,*) 'List of EOM Xi{O} vector calculations:' 409 DO ITRAN = 1, NEOMXE2TRAN 410 WRITE(LUPRI,'(A,5I5,5X,(25I3,20X))') MSGDBG, 411 & (IEOMXE2TRAN(I,ITRAN),I=1,5), 412 & (IEOMX2DOTS(I,ITRAN),I=1,MXEOME2VEC) 413 END DO 414 WRITE(LUPRI,*) 415 416 417 !extra L0 x RE vector dot products: 418 WRITE (LUPRI,*) 'List of L0 x RE dot products:' 419 DO ITRAN = 1, NEOML0TRAN 420 WRITE(LUPRI,'(A,I5,5X,(12I5,20X))') MSGDBG, 421 & IEOML0TRAN(ITRAN),(IEOML0DOTS(I,ITRAN),I=1,MXEOML0VEC) 422 END DO 423 WRITE (LUPRI,*) 424 425 END IF 426 427 RETURN 428 END 429 430*---------------------------------------------------------------------* 431* END OF SUBROUTINE CCXOPA_EOMSETUP * 432*---------------------------------------------------------------------* 433