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*=====================================================================* 20c /* deck cc_mcdcal */ 21*=====================================================================* 22 SUBROUTINE CC_MCDCAL(WORK,LWORK) 23*---------------------------------------------------------------------* 24* 25* Purpose: drive final calculation of B terms for MCD 26* 27* Sonia Coriani Fall 1997 28* Restructured January 2000. Sonia 29* Calculation of LAO-B term introduced in Spring 2000. Sonia 30* New output style March 2001. Sonia 31*=====================================================================* 32#if defined (IMPLICIT_NONE) 33 IMPLICIT NONE 34#else 35# include "implicit.h" 36#endif 37#include "priunit.h" 38#include "ccsdinp.h" 39#include "ccorb.h" 40#include "ccsdsym.h" 41#include "ccroper.h" 42C#include "ccropr2.h" 43#include "ccmcdinf.h" 44#include "ccexci.h" 45#include "cclists.h" 46#include "ccr1rsp.h" 47#include "ccpl1rsp.h" 48#include "ccel1rsp.h" 49#include "ccer1rsp.h" 50#include "cclrmrsp.h" 51#include "cco1rsp.h" 52#include "ccx1rsp.h" 53#include "cco2rsp.h" 54#include "dummy.h" 55* arguments 56 INTEGER LWORK 57#if defined (SYS_CRAY) 58 REAL WORK(LWORK) 59#else 60 DOUBLE PRECISION WORK(LWORK) 61#endif 62* local parameters: 63 CHARACTER*(19) MSGDBG 64 PARAMETER (MSGDBG = '[debug] CC_MCDCAL> ') 65 LOGICAL LOCDBG 66 PARAMETER (LOCDBG = .FALSE.) 67#if defined (SYS_CRAY) 68 REAL DDOT,DNRM2 69 REAL ZERO, D05, TWO 70 REAL FREQEX, XLEFTC, XRIGHTC, XRIGHTAB,XLEFTAB 71 REAL BTERM_A,BTERM_S, XMAB_MC, XMC_MAB 72 REAL FORZADIP, SQRFORZA 73#else 74 DOUBLE PRECISION DDOT,DNRM2 75 DOUBLE PRECISION ZERO, D05, TWO, ONE 76 DOUBLE PRECISION FREQEX, XLEFTC, XRIGHTC, XRIGHTAB,XLEFTAB 77 DOUBLE PRECISION BTERM_A,BTERM_S,XMAB_MC,XMC_MAB 78 DOUBLE PRECISION FORZADIP, SQRFORZA 79#endif 80 PARAMETER (ZERO=0.0D00, D05=0.5D00, TWO=2.0D00, ONE=1.0D00) 81C PARAMETER (DUMMY = -12345.67890) 82 83 CHARACTER MODEL*10, MODPRI*5 84 CHARACTER*8 LABELA, LABELB, LABELC, LABSOP 85 CHARACTER*8 CRLXA, CRLXB, CRLXC 86 LOGICAL LORXA, LORXB, LORXC, LPDBSA, LPDBSB, LPDBSC, LPROJ 87 INTEGER IOPERA,IOPERB,IOPERC,IOPSOP,ISGNSOP,INUM 88 INTEGER ISYMA,ISYMB,ISYMC,ISYSOP 89 INTEGER ITAMPA, ITAMPB, ITAMPC 90 INTEGER IKAPA, IKAPB, IKAPC, ISYMAB 91 INTEGER KEND0, LEND0, KEND1, LEND1 92 INTEGER MXTRAN, MXVEC, ITRAN, IVEC, IORDER 93 INTEGER IOPT, IOPTRD, IFREQ, ISYM, ISYML, IOPER 94 INTEGER KBTERMS 95 INTEGER KGTRAN, KGDOTS, KGCONS, NGTRAN, MXGTRAN, MXGDOTS 96 INTEGER KFATRAN, KFADOTS, KFACONS, NFATRAN, MXFATRAN, MXFADOTS 97 INTEGER KEATRAN, KEADOTS, KEACONS, NEATRAN, MXEATRAN, MXEADOTS 98 INTEGER KEATRA1, KEADOT1, KEACON1, NEATRA1 99 INTEGER KEATRA2, KEADOT2, KEACON2, NEATRA2 100 INTEGER KFTRAN, KFDOTS, KFCONS, NFTRAN, MXFTRAN, MXFDOTS 101 INTEGER KFTRA1, KFDOT1, KFCON1, NFTRA1 102 INTEGER KXE2TRA, KX2DOTS, KE2DOTS, KX2CONS,KE2CONS 103 INTEGER NXE2TRA, MXXETRAN,MXXEDOTS 104 INTEGER KX2TRAN, KX2DOT1, KX2CON1, NX2TRAN 105 INTEGER IEXCIF, ISTATE, ISYMSF, ISTATF 106 INTEGER NTAMPA, NTAMPB, NTAMPC, NBCONTR 107 108 INTEGER IZETA0,IO2AB,IM1F,IER1A,IER1B,IEO1B,IEX1B,IEL1B,IPL1A 109 INTEGER IETAA,IETAB,IXKSIA,IXKSIB,IFTAMA,IFTAMB 110 INTEGER IETAC,IXKSIC,IFTAMC 111 112 INTEGER KRE_1,KRE_2,KM1F_1, KM1F_2, KTAMA_1,KTAMA_2 113 INTEGER KEX1B_1,KEX1B_2,KER1B_1,KER1B_2,KER1A_1, KER1A_2 114 INTEGER KPL1A_1,KPL1A_2,KEL1B_1,KEL1B_2,KFTA__1,KFTA__2 115 INTEGER KEO1B_1,KEO1B_2,KFTB__1,KFTB__2 116 INTEGER KLE_1,KLE_2,KO2AB_1,KO2AB_2 117 INTEGER KETAA_1,KETAA_2,KETAB_1,KETAB_2 118 INTEGER KXKSIA_1,KXKSIA_2 119 INTEGER KXKSIC_1,KXKSIC_2,KETAC_1,KETAC_2 120 121#if defined (SYS_CRAY) 122 REAL GCON, FACON1, FACON2, BCON1, BCON2 123 REAL EACON1, EACON2, EACON3, EACON4 124 REAL DOTCON1, DOTCON2, DOTCON3, DOTCON4, DOTCON5 125 REAL DWFB0, XE2CON1, XE2CON2, XE2CON3 126#else 127 DOUBLE PRECISION GCON, FACON1, FACON2, BCON1, BCON2 128 DOUBLE PRECISION EACON1, EACON2, EACON3, EACON4 129 DOUBLE PRECISION DOTCON1, DOTCON2, DOTCON3, DOTCON4, DOTCON5 130 DOUBLE PRECISION DWFB0, XE2CON1, XE2CON2, XE2CON3 131#endif 132* external functions 133 INTEGER IROPER 134 INTEGER IR1TAMP 135 INTEGER IR1KAPPA 136 INTEGER IRHSR2 137 INTEGER IER1AMP 138 INTEGER IEL1AMP 139 INTEGER ILRMAMP 140 INTEGER IETA1 141 INTEGER IRHSR1 142 INTEGER IPL1ZETA 143 144*---------------------------------------------------------------------* 145* print header for magnetic circular dichroism section 146*---------------------------------------------------------------------* 147 WRITE (LUPRI,'(7(/1X,2A),/)') 148 & '************************************', 149 & '*******************************', 150 & '* ', 151 & ' *', 152 & '*-------- OUTPUT FROM COUPLED CLUSTER Q', 153 & 'UADRATIC RESPONSE ---------*', 154 & '* ', 155 & ' *', 156 & '*-------- CALCULATION OF B TERMS IN MAG', 157 & 'NETIC CIRCULAR D. ---------*', 158 & '* ', 159 & ' *', 160 & '************************************', 161 & '*******************************' 162 163*---------------------------------------------------------------------* 164* print some debug/info output 165*---------------------------------------------------------------------* 166 IF (LOCDBG) THEN 167 WRITE(LUPRI,*) MSGDBG, ' Workspace = ',LWORK 168 WRITE(LUPRI,*) MSGDBG, ' LUSEPL1 = ', LUSEPL1 169 END IF 170*---------------------------------------------------------------------* 171* set IOPTRD for calls to CC_RDRSP 172*---------------------------------------------------------------------* 173 IF (CCS) THEN 174 MODEL = 'CCS ' 175 IOPTRD = 1 176 MODPRI = 'CCS ' 177 ELSE IF (CC2) THEN 178 MODEL = 'CC2 ' 179 IOPTRD = 3 180 MODPRI = 'CC2 ' 181 ELSE IF (CCSD) THEN 182 MODEL = 'CCSD ' 183 IOPTRD = 3 184 MODPRI = 'CCSD ' 185 ELSE 186 CALL QUIT('Unknown coupled cluster model in CC_MCD') 187 END IF 188*---------------------------------------------------------------------* 189* allocate & initialize work space for B terms 190*---------------------------------------------------------------------* 191 LPROJ = .FALSE. 192 193 NBCONTR = NMCDOPER*NMCDST 194 195 !max number of transformations: conservative estimate 196 MXTRAN = NLRTLBL * MAX(NLRM,NPL1LBL,NLRTLBL) 197 !max number of vector products.... 198 MXVEC = MAX(NLRTLBL,NER1LBL,NEL1LBL,NLRM,NO2LBL,NPL1LBL,NX1LBL) 199 200 MXGTRAN = MXDIM_GTRAN * MXTRAN 201 MXFATRAN = MXDIM_FATRAN * MXTRAN 202 MXEATRAN = MXDIM_XEVEC * MXTRAN 203 !for the use-left case 204 MXFTRAN = MXDIM_FTRAN * MXTRAN 205 !for the additional contrib from relax 206 MXXETRAN = MXDIM_XEVEC * MXTRAN 207 208 MXGDOTS = MXVEC * MXTRAN 209 MXFADOTS = MXVEC * MXTRAN 210 MXEADOTS = MXVEC * MXTRAN 211 !for the use-left case 212 MXFDOTS = MXVEC * MXTRAN 213 !for the additional contrib from relax 214 MXXEDOTS = MXVEC * MXTRAN 215 216 KBTERMS = 1 217 KGTRAN = KBTERMS + NBCONTR 218 KGDOTS = KGTRAN + MXGTRAN 219 KGCONS = KGDOTS + MXGDOTS 220 KFATRAN = KGCONS + MXGDOTS 221 KFADOTS = KFATRAN + MXFATRAN 222 KFACONS = KFADOTS + MXFADOTS 223 KEATRAN = KFACONS + MXFADOTS 224 KEADOTS = KEATRAN + MXEATRAN 225 KEACONS = KEADOTS + MXEADOTS 226 KEND0 = KEACONS + MXEADOTS 227 LEND0 = LWORK - KEND0 228 229 !additional allocations for integer arrays ETA{A} and BMAT 230 KEATRA1 = KEND0 231 KEADOT1 = KEATRA1 + MXEATRAN 232 KEACON1 = KEADOT1 + MXEADOTS 233 KEATRA2 = KEACON1 + MXEADOTS 234 KEADOT2 = KEATRA2 + MXEATRAN 235 KEACON2 = KEADOT2 + MXEADOTS 236 KFTRAN = KEACON2 + MXEADOTS 237 KFDOTS = KFTRAN + MXFTRAN 238 KFCONS = KFDOTS + MXFDOTS 239 KFTRA1 = KFCONS + MXFDOTS 240 KFDOT1 = KFTRA1 + MXFTRAN 241 KFCON1 = KFDOT1 + MXFDOTS 242 KEND0 = KFCON1 + MXFDOTS 243 244 !extra contributions from PDBS 245 KXE2TRA = KEND0 246 KX2DOTS = KXE2TRA + MXXETRAN !Xksi^{A^B} contrib. 247 KX2CONS = KX2DOTS + MXXEDOTS 248 KE2DOTS = KX2CONS + MXXEDOTS !Eta^{A^B} contrib. 249 KE2CONS = KE2DOTS + MXXEDOTS 250 KX2TRAN = KE2CONS + MXXEDOTS 251 KX2DOT1 = KX2TRAN + MXXETRAN 252 KX2CON1 = KX2DOT1 + MXXEDOTS 253 KEND0 = KX2CON1 + MXXEDOTS 254 LEND0 = LWORK - KEND0 255 256 IF (LEND0 .LT. 0 ) THEN 257 CALL QUIT('Insufficient work space in CC_MCDCAL (1)') 258 END IF 259 260 CALL DZERO(WORK(KBTERMS),NBCONTR) 261 262*---------------------------------------------------------------------* 263* set up lists for G and F{A} transformations and ETA{A} vectors: 264* as well as B (generalized F) 265*---------------------------------------------------------------------* 266 CALL CCMCD_SETUP(MXTRAN, MXVEC, 267 & WORK(KGTRAN), WORK(KGDOTS), NGTRAN, 268 & WORK(KFATRAN),WORK(KFADOTS),NFATRAN, 269 & WORK(KEATRAN),WORK(KEADOTS),NEATRAN, 270 & WORK(KEATRA1),WORK(KEADOT1),NEATRA1, 271 & WORK(KEATRA2),WORK(KEADOT2),NEATRA2, 272 & WORK(KFTRAN), WORK(KFDOTS), NFTRAN, 273 & WORK(KFTRA1), WORK(KFDOT1), NFTRA1, 274 & WORK(KXE2TRA),WORK(KX2DOTS), 275 & WORK(KE2DOTS),NXE2TRA, 276 & WORK(KX2TRAN),WORK(KX2DOT1),NX2TRAN, 277 & WORK(KEND0), LEND0 ) 278*---------------------------------------------------------------------* 279* calculate G matrix contributions: (G*T^A*T^B)*E^f GCON 280*---------------------------------------------------------------------* 281 IOPT = 5 !calculate ddot product 282 CALL CC_GMATRIX('L0 ','R1 ','R1 ','RE ',NGTRAN, MXVEC, 283 & WORK(KGTRAN),WORK(KGDOTS),WORK(KGCONS), 284 & WORK(KEND0), LEND0, IOPT ) 285*---------------------------------------------------------------------* 286* calculate F{A} matrix contributions: (F^A*T^B)*RE and (F^B*T^A)*RE 287* FACON1 FACON2 288*---------------------------------------------------------------------* 289 CALL CCQR_FADRV('L0 ','o1 ','R1 ','RE ',NFATRAN, MXVEC, 290 & WORK(KFATRAN),WORK(KFADOTS),WORK(KFACONS), 291 & WORK(KEND0), LEND0, 'DOTP' ) 292*---------------------------------------------------------------------* 293* calculate ETA{A} vector contributions: 294* ETAA(LE)*TB (always) EACON1 295* ETAA(M1)*TB,ETAB(M1)*TA EACON2, EACON3 296* ETAB(ZA)*RE EACON4 297*---------------------------------------------------------------------* 298 CALL DZERO(WORK(KEACONS),MXEADOTS) 299 IOPT = 5 300 IORDER = 1 301 CALL CC_XIETA( WORK(KEATRAN), NEATRAN, IOPT, IORDER, 'LE ', 302 & '---',IDUMMY, DUMMY, 303 & 'R1 ',WORK(KEADOTS),WORK(KEACONS), 304 & .FALSE.,MXVEC, WORK(KEND0), LEND0 ) 305 306 IF (LUSEPL1) THEN 307 CALL DZERO(WORK(KEACON1),MXEADOTS) 308 CALL DZERO(WORK(KEACON2),MXEADOTS) 309 IOPT = 5 310 IORDER = 1 311 CALL CC_XIETA( WORK(KEATRA1), NEATRA1, IOPT, IORDER, 'M1 ', 312 & '---',IDUMMY, DUMMY, 313 & 'R1 ',WORK(KEADOT1),WORK(KEACON1), 314 & .FALSE.,MXVEC, WORK(KEND0), LEND0 ) 315 CALL FLSHFO(LUPRI) 316 IOPT = 5 317 IORDER = 1 318 CALL CC_XIETA( WORK(KEATRA2), NEATRA2, IOPT, IORDER, 'PL1', 319 & '---',IDUMMY, DUMMY, 320 & 'RE ',WORK(KEADOT2),WORK(KEACON2), 321 & .FALSE.,MXVEC, WORK(KEND0), LEND0 ) 322 323 END IF 324 325*---------------------------------------------------------------------* 326* calculate B matrix contributions (use generalized F matrix): 327* (M1*B*TA)*TB BCON1 328* (ZA*B*TB)*RE BCON2 329*---------------------------------------------------------------------* 330 IF (LUSEPL1) THEN 331 IOPT = 5 332 CALL CC_FMATRIX(WORK(KFTRAN),NFTRAN,'M1 ','R1 ',IOPT, 333 & 'R1 ',WORK(KFDOTS),WORK(KFCONS),MXVEC, 334 & WORK(KEND0), LEND0) 335 336 IOPT = 5 337 CALL CC_FMATRIX(WORK(KFTRA1),NFTRA1,'PL1','R1 ',IOPT, 338 & 'RE ',WORK(KFDOT1),WORK(KFCON1),MXVEC, 339 & WORK(KEND0), LEND0) 340 END IF 341*---------------------------------------------------------------------* 342* calculate PDBS extra dot product contributions: 343* M1 * Xksi{O2} and ETA{O2} * RE (for left moment) XECON1,XECON2 344* LE * Xksi{O2} (for right moment) XECON3 345* It would be better to have Xksi{O2} of file and reuse it... 346*---------------------------------------------------------------------* 347 CALL DZERO(WORK(KX2CONS),MXXEDOTS) 348 CALL DZERO(WORK(KE2CONS),MXXEDOTS) 349 CALL DZERO(WORK(KX2CON1),MXXEDOTS) 350 351 IF (LUSEPL1) THEN 352 IOPT = 5 353 IORDER = 2 354 CALL CC_XIETA( WORK(KXE2TRA), NXE2TRA, IOPT, IORDER, 'L0 ', 355 & 'M1 ',WORK(KX2DOTS),WORK(KX2CONS), 356 & 'RE ',WORK(KE2DOTS),WORK(KE2CONS), 357 & .FALSE.,MXVEC, WORK(KEND0), LEND0 ) 358 CALL CC_XIETA( WORK(KX2TRAN), NX2TRAN, IOPT, IORDER, 'L0 ', 359 & 'LE ',WORK(KX2DOT1),WORK(KX2CON1), 360 & '---',IDUMMY,DUMMY, 361 & .FALSE.,MXVEC, WORK(KEND0), LEND0 ) 362 END IF 363 364*---------------------------------------------------------------------* 365C Gather the different contributions and calculate remaining 366C dot products: 367* M1 * R1 DOTCON1 M^f*T^A 368* EX1 * RE DWFB0 DeltaW^B(0) 369* PL1 * RE DOTCON2 370* LE * EO1 DWFB0 DeltaW^B(0) 371* M1 * O2 DOTCON3 M^f*Xksi^AB 372C---------------------------------------------------------------------* 373 374 DO 100 IOPER = 1, NMCDOPER 375 376 LPROJ = .FALSE. 377 378 IOPERA = IAMCDOP(IOPER) 379 IOPERB = IBMCDOP(IOPER) 380 IOPERC = ICMCDOP(IOPER) 381 382 LABELA = LBLOPR(IAMCDOP(IOPER)) 383 LABELB = LBLOPR(IBMCDOP(IOPER)) 384 LABELC = LBLOPR(ICMCDOP(IOPER)) 385 386 ISYMA = ISYOPR(IAMCDOP(IOPER)) 387 ISYMB = ISYOPR(IBMCDOP(IOPER)) 388 ISYMC = ISYOPR(ICMCDOP(IOPER)) 389 ISYMAB = MULD2H(ISYMA,ISYMB) 390 391 IF (LOCDBG) THEN 392 WRITE(LUPRI,*) MSGDBG, 'A oper: ',LABELA, ' symm: ', ISYMA 393 WRITE(LUPRI,*) MSGDBG, 'B oper: ',LABELB, ' symm: ', ISYMB 394 WRITE(LUPRI,*) MSGDBG, 'C oper: ',LABELC, ' symm: ', ISYMC 395 END IF 396 397 CALL FLSHFO(LUPRI) 398 399 IF (ISYMAB.EQ.ISYMC) THEN 400 401 LORXA = LAMCDRX(IOPER) 402 LORXB = LBMCDRX(IOPER) 403 LORXC = LCMCDRX(IOPER) 404 405 LPDBSA = LPDBSOP(IOPERA) 406 LPDBSB = LPDBSOP(IOPERB) 407 LPDBSC = LPDBSOP(IOPERC) 408 409 IF ( LORXA) CRLXA = '(relax.)' 410 IF ( .NOT. LORXA) CRLXA = '(unrel.)' 411 IF ( LORXB) CRLXB = '(relax.)' 412 IF ( .NOT. LORXB) CRLXB = '(unrel.)' 413 IF ( LORXC) CRLXC = '(relax.)' 414 IF ( .NOT. LORXC) CRLXC = '(unrel.)' 415 416 417 DO ISTATE = 1, NMCDST 418 ISYMSF = IMCDSTSY(ISTATE) 419 ISTATF = IMCDSTNR(ISTATE) !RE,LE index within sym 420 IEXCIF = ISYOFE(ISYMSF) + ISTATF !RE,LE absolute index 421 FREQEX = EIGVAL(IEXCIF) 422 423 IF (ISYMSF.EQ.ISYMC) THEN 424 425C------------------------------------------------------------------ 426C Calculate residues. 427C - request indices of vectors in the vector list 428C------------------------------------------------------------------ 429 430 IZETA0 = 0 431 ITAMPA = IR1TAMP(LABELA,LORXA,-FREQEX,ISYMA) 432 ITAMPB = IR1TAMP(LABELB,LORXB, ZERO, ISYMB) 433 IKAPA = 0 434 IKAPB = 0 435 IF(LORXA)IKAPA=IR1KAPPA(LABELA,-FREQEX,ISYMA) 436 IF(LORXB)IKAPB=IR1KAPPA(LABELB, ZERO, ISYMB) 437 438 IF (.NOT.LUSE2N1) THEN 439 ITAMPC = IR1TAMP(LABELC,LORXC, FREQEX,ISYMC) 440 IKAPC = 0 441 IF (LORXC) IKAPC = IR1KAPPA(LABELC, FREQEX,ISYMC) 442 IFTAMC = ITAMPC 443 END IF 444 445 IM1F = ILRMAMP(IEXCIF,FREQEX,ISYMC) 446 IER1A = IER1AMP(IEXCIF,FREQEX,ISYMC, 447 & LABELA,-FREQEX,ISYMA,.FALSE.) 448 IETAB = IETA1(LABELB,LORXB,ZERO,ISYMB) 449 IETAC = IETA1(LABELC,LORXC,FREQEX,ISYMC) 450 IFTAMB = ITAMPB 451 IF (ISYMB .EQ. 1) LPROJ = .TRUE. 452 IEL1B = IEL1AMP(IEXCIF,FREQEX,ISYMC, 453 & LABELB,ZERO,ISYMB,LORXB,LPROJ) 454 IEX1B = IEL1B 455 IXKSIA = IRHSR1(LABELA,LORXA,-FREQEX,ISYMA) 456 IXKSIC = IRHSR1(LABELC,LORXC,FREQEX,ISYMC) 457 458 IF (LUSEPL1) THEN 459 IF (ISYMB .EQ. 1) LPROJ = .TRUE. 460 IPL1A = IPL1ZETA(LABELA,LORXA,-FREQEX,ISYMA, 461 & LPROJ,IEXCIF, FREQEX,ISYMC) 462 ELSE 463 IF (ISYMB .EQ. 1) LPROJ = .TRUE. 464 IER1B = IER1AMP(IEXCIF,FREQEX,ISYMC, 465 & LABELB,ZERO,ISYMB,LPROJ) 466 IEO1B = IER1B 467 IO2AB = IRHSR2(LABELA,LORXA,-FREQEX,ISYMA, 468 & LABELB,LORXB, ZERO,ISYMB) 469 IETAA = IETA1(LABELA,LORXA,-FREQEX,ISYMA) 470 IFTAMA = ITAMPA 471 472 END IF 473*-----------------------------------------------------------------* 474* Calculate linear response residues. * 475*-----------------------------------------------------------------* 476 XLEFTC = ZERO 477 XRIGHTC = ZERO 478 FORZADIP = ZERO 479 SQRFORZA = ZERO 480 481 NTAMPC = NT1AM(ISYMC) + NT2AM(ISYMC) 482 IF (CCS) NTAMPC = NT1AM(ISYMC) 483 484 KXKSIC_1 = KEND0 485 KXKSIC_2 = KXKSIC_1 + NT1AM(ISYMC) 486 KLE_1 = KXKSIC_2 + NT2AM(ISYMC) 487 KLE_2 = KLE_1 + NT1AM(ISYMC) 488 KEND1 = KLE_2 + NT2AM(ISYMC) 489 LEND1 = LWORK - KEND1 490 491 IF (LEND1.LE.0) THEN 492 WRITE(LUPRI,*)'Insuff. work in CC_MCDCAL (LR1)' 493 WRITE(LUPRI,*)'Need: ', KEND1, 'Available: ',LWORK 494 CALL QUIT('Insufficient work space in CC_MCDCAL') 495 END IF 496 497 CALL CC_RDRSP('O1 ',IXKSIC,ISYMC,IOPTRD,MODEL, 498 & WORK(KXKSIC_1), WORK(KXKSIC_2)) 499 CALL CC_RDRSP('LE ',IEXCIF,ISYMC,IOPTRD,MODEL, 500 & WORK(KLE_1), WORK(KLE_2)) 501 502 XRIGHTC = DDOT(NTAMPC,WORK(KLE_1),1,WORK(KXKSIC_1),1) 503 504 IF (LUSE2N1) THEN 505 506 KM1F_1 = KLE_1 507 KM1F_2 = KM1F_1 + NT1AM(ISYMC) 508 KEND1 = KM1F_2 + NT2AM(ISYMC) 509 LEND1 = LWORK - KEND1 510 511 CALL CC_RDRSP('M1 ',IM1F,ISYMC,IOPTRD,MODEL, 512 & WORK(KM1F_1), WORK(KM1F_2)) 513 514 XLEFTC = DDOT(NTAMPC,WORK(KM1F_1),1,WORK(KXKSIC_1),1) 515 516 KETAC_1 = KM1F_1 517 KETAC_2 = KETAC_1 + NT1AM(ISYMC) 518 KRE_1 = KETAC_2 + NT2AM(ISYMC) 519 KRE_2 = KRE_1 + NT1AM(ISYMC) 520 KEND1 = KRE_2 + NT2AM(ISYMC) 521 LEND1 = LWORK - KEND1 522 523 CALL CC_RDRSP('RE ',IEXCIF,ISYMC,IOPTRD,MODEL, 524 & WORK(KRE_1), WORK(KRE_2)) 525 CALL CC_RDRSP('X1 ',IETAC,ISYMC,IOPTRD,MODEL, 526 & WORK(KETAC_1), WORK(KETAC_2)) 527 XLEFTC = XLEFTC + 528 & DDOT(NTAMPC,WORK(KETAC_1),1,WORK(KRE_1),1) 529 ELSE 530 !Unfinished 531 WRITE(LUPRI,*)'MCD: LUSE2N1=FALSE not yet implemented' 532 CALL QUIT('MCD: LUSE2N1=FALSE not yet implemented') 533 END IF 534c 535c Dipole strength: S^of_C,C = 0.5(M^C_of M^C_fo+[M^C_of M^C_fo]^*) 536c 537 FORZADIP = XLEFTC*XRIGHTC 538c 539 IF (FORZADIP.GE.ZERO) THEN 540 SQRFORZA = SQRT(FORZADIP) 541 ELSE 542 SQRFORZA = -SQRT(-FORZADIP) 543 ENDIF 544C------------------------------------------------------------------ 545C Calculate quadratic response residues. 546C------------------------------------------------------------------ 547C - LEFT PART 548C - RIGHT PART 549C the call to *SET* routines recover ITRAN,IVEC 550*----------------------------------------------------------------------* 551* First contribution: {G*T^A(-w_f)*T^B(0)}*E^f (see at the beginning) 552*----------------------------------------------------------------------* 553 CALL CC_SETG112(WORK(KGTRAN),WORK(KGDOTS),NGTRAN,MXVEC, 554 & IZETA0,ITAMPA,ITAMPB,IEXCIF,ITRAN,IVEC) 555 GCON = WORK(KGCONS-1 + (ITRAN-1)*MXVEC + IVEC) 556*----------------------------------------------------------------------* 557* Second contribution: {F^A * T^B} * E^f 558*----------------------------------------------------------------------* 559c CALL CC_SETFA12(WORK(KFATRAN),WORK(KFADOTS),NFATRAN,MXVEC, 560c & IZETA0,IOPERA,ITAMPB,IEXCIF,ITRAN,IVEC) 561 CALL CC_SETFB12(WORK(KFATRAN),WORK(KFADOTS),NFATRAN,MXVEC, 562 & IZETA0,IOPERA,IKAPA,ITAMPB,IEXCIF,ITRAN,IVEC) 563 FACON1 = WORK(KFACONS-1 + (ITRAN-1)*MXVEC + IVEC) 564*----------------------------------------------------------------------* 565* Third contribution: {F^B*T^A} * E^f 566*----------------------------------------------------------------------* 567c CALL CC_SETFA12(WORK(KFATRAN),WORK(KFADOTS),NFATRAN,MXVEC, 568c & IZETA0,IOPERB,ITAMPA,IEXCIF,ITRAN,IVEC) 569 CALL CC_SETFB12(WORK(KFATRAN),WORK(KFADOTS),NFATRAN,MXVEC, 570 & IZETA0,IOPERB,IKAPB,ITAMPA,IEXCIF,ITRAN,IVEC) 571 FACON2 = WORK(KFACONS-1 + (ITRAN-1)*MXVEC + IVEC) 572*----------------------------------------------------------------------* 573* Fourth contribution M^f*{various terms} (excluding relaxation and DW)* 574*----------------------------------------------------------------------* 575 IF (LUSEPL1) THEN 576 CALL CC_SETXE('Eta',WORK(KEATRA1),WORK(KEADOT1), 577 & NEATRA1,MXVEC,IM1F,IOPERA,IKAPA,0,0,0, 578 & ITAMPB,ITRAN,IVEC) 579 EACON2 = WORK(KEACON1-1 + (ITRAN-1)*MXVEC + IVEC) 580* 581 CALL CC_SETXE('Eta',WORK(KEATRA1),WORK(KEADOT1), 582 & NEATRA1,MXVEC,IM1F,IOPERB,IKAPB,0,0,0, 583 & ITAMPA,ITRAN,IVEC) 584 EACON3 = WORK(KEACON1-1 + (ITRAN-1)*MXVEC + IVEC) 585* 586 CALL CCQR_SETF(WORK(KFTRAN),WORK(KFDOTS),NFTRAN, 587 & MXVEC,IM1F,ITAMPA,ITAMPB,ITRAN,IVEC) 588 BCON1 = WORK(KFCONS-1 + (ITRAN-1)*MXVEC + IVEC) 589* 590 ELSE 591 NTAMPC = NT1AM(ISYMC) + NT2AM(ISYMC) 592 IF (CCS) NTAMPC = NT1AM(ISYMC) 593 KM1F_1 = KEND0 594 KM1F_2 = KM1F_1 + NT1AM(ISYMC) 595 KO2AB_1 = KM1F_2 + NT2AM(ISYMC) 596 KO2AB_2 = KO2AB_1 + NT1AM(ISYMC) 597 KEND1 = KO2AB_2 + NT2AM(ISYMC) 598 LEND1 = LWORK - KEND1 599 600 IF (LEND1.LE.0) THEN 601 WRITE(LUPRI,*)'Insuff. work sp. in CC_MCDCAL (Mf)' 602 WRITE(LUPRI,*)'Need: ', KEND1,'Available: ', LWORK 603 CALL QUIT('Insufficient work space in CC_MCDCAL') 604 END IF 605 606 CALL CC_RDRSP('M1 ',IM1F,ISYMC,IOPTRD,MODEL, 607 & WORK(KM1F_1), WORK(KM1F_2)) 608 609 CALL CC_RDRSP('O2 ',IO2AB,ISYMC,IOPTRD,MODEL, 610 & WORK(KO2AB_1), WORK(KO2AB_2)) 611 CALL CCLR_DIASCL(WORK(KO2AB_2),D05,ISYMC) 612 DOTCON3 = DDOT(NTAMPC,WORK(KM1F_1),1,WORK(KO2AB_1),1) 613 614 END IF 615*----------------------------------------------------------------------* 616* Fifth contribution: Xksibar^A(-wf) * E^fB (if not LUSEPL1: DOTCON4) 617*----------------------------------------------------------------------* 618 IF (LUSEPL1) THEN 619 620 CALL CC_SETF12(WORK(KFTRA1),WORK(KFDOT1),NFTRA1,MXVEC, 621 & IPL1A,ITAMPB,IEXCIF,ITRAN,IVEC) 622 BCON2 = WORK(KFCON1-1 + (ITRAN-1)*MXVEC + IVEC) 623 624 CALL CC_SETXE('Eta',WORK(KEATRA2),WORK(KEADOT2), 625 & NEATRA2,MXVEC,IPL1A,IOPERB,IKAPB,0,0,0, 626 & IEXCIF,ITRAN,IVEC) 627 EACON4 = WORK(KEACON2-1 + (ITRAN-1)*MXVEC + IVEC) 628 ELSE 629 NTAMPA = NT1AM(ISYMA) + NT2AM(ISYMA) 630 IF (CCS) NTAMPA = NT1AM(ISYMA) 631 KER1B_1 = KEND0 632 KER1B_2 = KER1B_1 + NT1AM(ISYMA) 633 KETAA_1 = KER1B_2 + NT2AM(ISYMA) 634 KETAA_2 = KETAA_1 + NT1AM(ISYMA) 635 KFTA__1 = KETAA_2 + NT2AM(ISYMA) 636 KFTA__2 = KFTA__1 + NT1AM(ISYMA) 637 KEND1 = KFTA__2 + NT2AM(ISYMA) 638 LEND1 = LWORK - KEND1 639 640 IF (LEND1.LE.0) THEN 641 WRITE(LUPRI,*)'Insuff. work in CC_MCDCAL (XibA)' 642 WRITE(LUPRI,*)'Need: ', KEND1, 'Available: ', LWORK 643 CALL QUIT('Insufficient work space in CC_MCDCAL') 644 END IF 645 646 CALL CC_RDRSP('ER1',IER1B,ISYMA,IOPTRD,MODEL, 647 & WORK(KER1B_1), WORK(KER1B_2)) 648 CALL CC_RDRSP('X1 ',IETAA,ISYMA,IOPTRD,MODEL, 649 & WORK(KETAA_1), WORK(KETAA_2)) 650 CALL CC_RDRSP('F1 ',IFTAMA,ISYMA,IOPTRD,MODEL, 651 & WORK(KFTA__1), WORK(KFTA__2)) 652 CALL DAXPY(NTAMPA,ONE,WORK(KFTA__1),1,WORK(KETAA_1),1) 653 654 DOTCON4 = DDOT(NTAMPA,WORK(KETAA_1),1,WORK(KER1B_1),1) 655 656 END IF 657*----------------------------------------------------------------------* 658* Sixth contribution Xksibar^B(-wf) * E^fA DOTCON5 659*----------------------------------------------------------------------* 660 NTAMPB = NT1AM(ISYMB) + NT2AM(ISYMB) 661 IF (CCS) NTAMPB = NT1AM(ISYMB) 662 KER1A_1 = KEND0 663 KER1A_2 = KER1A_1 + NT1AM(ISYMB) 664 KETAB_1 = KER1A_2 + NT2AM(ISYMB) 665 KETAB_2 = KETAB_1 + NT1AM(ISYMB) 666 KFTB__1 = KETAB_2 + NT2AM(ISYMB) 667 KFTB__2 = KFTB__1 + NT1AM(ISYMB) 668 KEND1 = KFTB__2 + NT2AM(ISYMB) 669 LEND1 = LWORK - KEND1 670 671 IF (LEND1.LT.NT2AM(ISYMA)) THEN 672 WRITE(LUPRI,*)'Insuff. work sp. in CC_MCDCAL (XibB*)' 673 WRITE(LUPRI,*)'Need: ', KEND1,'Available: ', LWORK 674 CALL QUIT('Insufficient work space in CC_MCDCAL') 675 END IF 676 677 CALL CC_RDRSP('ER1',IER1A,ISYMB,IOPTRD,MODEL, 678 & WORK(KER1A_1), WORK(KER1A_2)) 679 680 CALL CC_RDRSP('X1 ',IETAB,ISYMB,IOPTRD,MODEL, 681 & WORK(KETAB_1), WORK(KETAB_2)) 682 CALL CC_RDRSP('F1 ',IFTAMB,ISYMB,IOPTRD,MODEL, 683 & WORK(KFTB__1), WORK(KFTB__2)) 684 CALL DAXPY(NTAMPB,ONE,WORK(KFTB__1),1,WORK(KETAB_1),1) 685 686 DOTCON5 = DDOT(NTAMPB,WORK(KETAB_1),1,WORK(KER1A_1),1) 687 688*----------------------------------------------------------------------* 689* Extra contribution: {DeltaW}_of^B(0) * Mbar^f(w_f) * T^A(-w_f) 690* and : {DeltaW}_of^B(0) * PL1A * RE 691* 692* where DeltaW^B_of(0) = Ebar^f(-wf) * Csi^fB(w_f) = 693* = Ebar^f(-wf) * (B*T^B + A^B) 694* ---------------------------------------------------------------------* 695 IF (ISYMB.EQ.1) THEN 696 697 NTAMPA = NT1AM(ISYMA) + NT2AM(ISYMA) 698 IF (CCS) NTAMPA = NT1AM(ISYMA) 699 KM1F_1 = KEND0 700 KM1F_2 = KM1F_1 + NT1AM(ISYMA) 701 KTAMA_1 = KM1F_2 + NT2AM(ISYMA) 702 KTAMA_2 = KTAMA_1 + NT1AM(ISYMA) 703 KEND1 = KTAMA_2 + NT2AM(ISYMA) 704 LEND1 = LWORK - KEND1 705 706 IF (LEND1.LE.0) THEN 707 WRITE(LUPRI,*)'Insuff. work sp.in CC_MCDCAL (DeltaW)' 708 WRITE(LUPRI,*)'Need: ', KEND1,'Available: ', LWORK 709 CALL QUIT('Insufficient work space in CC_MCDCAL') 710 END IF 711 712 CALL CC_RDRSP('M1 ',IM1F,ISYMA,IOPTRD,MODEL, 713 & WORK(KM1F_1), WORK(KM1F_2)) 714 715 CALL CC_RDRSP('R1 ',ITAMPA,ISYMA,IOPTRD,MODEL, 716 & WORK(KTAMA_1), WORK(KTAMA_2)) 717 DOTCON1 = DDOT(NTAMPA,WORK(KM1F_1),1,WORK(KTAMA_1),1) 718 719 IF (LUSEPL1) THEN 720 NTAMPA = NT1AM(ISYMA) + NT2AM(ISYMA) 721 IF (CCS) NTAMPA = NT1AM(ISYMA) 722 KRE_1 = KEND0 723 KRE_2 = KRE_1 + NT1AM(ISYMA) 724 KEX1B_1 = KRE_2 + NT2AM(ISYMA) 725 KEX1B_2 = KEX1B_1 + NT1AM(ISYMA) 726 KPL1A_1 = KEX1B_2 + NT2AM(ISYMA) 727 KPL1A_2 = KPL1A_1 + NT1AM(ISYMA) 728 KEND1 = KPL1A_2 + NT2AM(ISYMA) 729 LEND1 = LWORK - KEND1 730 IF (LEND1.LE.0) THEN 731 WRITE(LUPRI,*)'Insuff. work sp.in CC_MCDCAL (LPL1)' 732 WRITE(LUPRI,*)'Need: ', KEND1,' Available: ', LWORK 733 CALL QUIT('Insufficient work space in CC_MCDCAL') 734 END IF 735 736 CALL CC_RDRSP('RE ',IEXCIF,ISYMA,IOPTRD,MODEL, 737 & WORK(KRE_1), WORK(KRE_2)) 738 739 CALL CC_RDRSP('EX1',IEX1B,ISYMA,IOPTRD,MODEL, 740 & WORK(KEX1B_1), WORK(KEX1B_2)) 741 742 CALL CC_RDRSP('PL1',IPL1A,ISYMA,IOPTRD,MODEL, 743 & WORK(KPL1A_1), WORK(KPL1A_2)) 744 745 DWFB0 = DDOT(NTAMPA,WORK(KEX1B_1),1,WORK(KRE_1),1) 746 DOTCON2 = DDOT(NTAMPA,WORK(KPL1A_1),1,WORK(KRE_1),1) 747 ELSE 748 NTAMPA = NT1AM(ISYMA) + NT2AM(ISYMA) 749 IF (CCS) NTAMPA = NT1AM(ISYMA) 750 KLE_1 = KEND0 751 KLE_2 = KRE_1 + NT1AM(ISYMA) 752 KEO1B_1 = KRE_2 + NT2AM(ISYMA) 753 KEO1B_2 = KEO1B_1 + NT1AM(ISYMA) 754 KEND1 = KEO1B_2 + NT2AM(ISYMA) 755 LEND1 = LWORK - KEND1 756 IF (LEND1.LE.0) THEN 757 WRITE(LUPRI,*)'Insuff. work sp.in CC_MCDCAL (NPL1)' 758 WRITE(LUPRI,*)'Need: ', KEND1,'Available: ', LWORK 759 CALL QUIT('Insufficient work space in CC_MCDCAL') 760 END IF 761 762 CALL CC_RDRSP('LE ',IEXCIF,ISYMA,IOPTRD,MODEL, 763 & WORK(KLE_1), WORK(KLE_2)) 764 765 CALL CC_RDRSP('EO1',IEO1B,ISYMA,IOPTRD,MODEL, 766 & WORK(KEO1B_1), WORK(KEO1B_2)) 767 768 DWFB0 = DDOT(NTAMPA,WORK(KLE_1),1,WORK(KEO1B_1),1) 769 DOTCON2 = ZERO 770 END IF 771 ELSE 772 DWFB0 = ZERO 773 DOTCON1 = ZERO 774 DOTCON2 = ZERO 775 END IF 776*----------------------------------------------------------------------* 777* PDBS/Relax contribution to Left 2nd order moment 778*----------------------------------------------------------------------* 779 IF ((LUSEPL1).AND.(LORXB.OR.LPDBSB)) THEN 780 CALL CC_FIND_SO_OP(LBLOPR(IOPERA),LBLOPR(IOPERB), 781 & LABSOP,ISYSOP,ISGNSOP,INUM, 782 & WORK(KEND0),LEND0) 783 IOPSOP = IROPER(LABSOP,ISYSOP) 784 CALL CC_SETXE('Xi ',WORK(KXE2TRA),WORK(KX2DOTS), 785 & MXTRAN,MXVEC, 786 & IZETA0,IOPSOP,IKAPA,IKAPB,0,0,IM1F, 787 & ITRAN,IVEC) 788 XE2CON1 = WORK(KX2CONS-1 + (ITRAN-1)*MXVEC + IVEC) 789 CALL CC_SETXE('Eta',WORK(KXE2TRA),WORK(KE2DOTS), 790 & MXTRAN,MXVEC, 791 & IZETA0,IOPSOP,IKAPA,IKAPB,0,0,IEXCIF, 792 & ITRAN,IVEC) 793 XE2CON2 = WORK(KE2CONS-1 + (ITRAN-1)*MXVEC + IVEC) 794 IF (LOCDBG) THEN 795 WRITE(LUPRI,*) MSGDBG, 796 & 'XE2CON1: ',XE2CON1,' XE2CON2: ',XE2CON2 797 END IF 798 ELSE 799 XE2CON1 = ZERO 800 XE2CON2 = ZERO 801 END IF 802 803*----------------------------------------------------------------------* 804* TOTAL LEFT TRANSITION M^AB{n<-f} * M^C{f<-n} = XMAB_MC 805*----------------------------------------------------------------------* 806 XLEFTAB = ZERO 807 XMAB_MC = ZERO 808c 809 XLEFTAB = GCON + FACON1 + FACON2 + DOTCON5 + 810 & DWFB0*(DOTCON1-DOTCON2) + 811 & XE2CON1 + XE2CON2 812 IF (LUSEPL1) THEN 813 XLEFTAB = XLEFTAB + EACON2 + EACON3 + BCON1 + 814 & EACON4 + BCON2 815 ELSE 816 XLEFTAB = XLEFTAB + DOTCON3 + DOTCON4 817 END IF 818 819 XMAB_MC = XLEFTAB * XRIGHTC 820 821 822************************************************************************ 823* RIGHT PART * 824************************************************************************ 825*----------------------------------------------------------------------* 826* First contribution: Ebar^fB(-w_f,0) * Xksi^A * 827*----------------------------------------------------------------------* 828 829 NTAMPA = NT1AM(ISYMA) + NT2AM(ISYMA) 830 KEL1B_1 = KEND0 831 KEL1B_2 = KEL1B_1 + NT1AM(ISYMA) 832 KXKSIA_1 = KEL1B_2 + NT2AM(ISYMA) 833 KXKSIA_2 = KXKSIA_1 + NT1AM(ISYMA) 834 KEND1 = KXKSIA_2 + NT2AM(ISYMA) 835 LEND1 = LWORK - KEND1 836 837 CALL CC_RDRSP('EL1',IEL1B,ISYMA,IOPTRD,MODEL, 838 & WORK(KEL1B_1),WORK(KEL1B_2)) 839 CALL CC_RDRSP('O1 ',IXKSIA,ISYMA,IOPTRD,MODEL, 840 & WORK(KXKSIA_1),WORK(KXKSIA_2)) 841 842C CALL CC_XKSI(WORK(KXKSIA_1),LABELA,ISYMA,DUMMY, 843C * WORK(KEND1),LEND1) 844 845 DOTCON5 = DDOT(NTAMPA,WORK(KXKSIA_1),1,WORK(KEL1B_1),1) 846* ---------------------------------------------------------------------* 847* Second contribution: Ebar^f * A^A * T^B(0) 848* computed as Eta^A(EL1) * T^B(0) 849* ---------------------------------------------------------------------* 850 CALL CC_SETXE('Eta',WORK(KEATRAN),WORK(KEADOTS),NEATRAN, 851 & MXVEC,IEXCIF,IOPERA,IKAPA,0,0,0,ITAMPB, 852 & ITRAN,IVEC) 853 EACON1 = WORK(KEACONS-1 + (ITRAN-1)*MXVEC + IVEC) 854*----------------------------------------------------------------------* 855* PDBS/Relax contribution to Right 2nd order moment 856*----------------------------------------------------------------------* 857 IF ((LUSEPL1).AND.(LORXB.OR.LPDBSB)) THEN 858c CALL CC_FIND_SO_OP(LBLOPR(IOPERA),LBLOPR(IOPERB), 859c & LABSOP,ISYSOP,ISGNSOP,INUM, 860c & WORK(KEND0),LEND0) 861 CALL CC_SETXE('Xi ',WORK(KX2TRAN),WORK(KX2DOT1), 862 & MXTRAN,MXVEC, 863 & IZETA0,IOPSOP,IKAPA,IKAPB,0,0,IEXCIF, 864 & ITRAN,IVEC) 865 XE2CON3 = WORK(KX2CON1-1 + (ITRAN-1)*MXVEC + IVEC) 866 IF (LOCDBG) THEN 867 WRITE(LUPRI,*) MSGDBG,'XE2CON3: ',XE2CON3 868 END IF 869 ELSE 870 XE2CON3 = ZERO 871 END IF 872* 873*----------------------------------------------------------------------* 874* TOTAL RIGHT TRANSITION M^C{o<-f} * M^AB{f<-o} = XMC_MAB 875*----------------------------------------------------------------------* 876* 877 XRIGHTAB = ZERO 878 XMC_MAB = ZERO 879 880 XRIGHTAB = DOTCON5 + EACON1 + XE2CON3 881 XMC_MAB = XLEFTC*XRIGHTAB 882*----------------------------------------------------------------------* 883* Debug information 884*----------------------------------------------------------------------* 885 IF (LOCDBG) THEN 886 WRITE(LUPRI,'(/1X,A,A,2(/1X,A,F16.8))') 887 * 'For C operator ', LABELC, 888 * ' Left 1st single moment = ', XLEFTC, 889 * ' Right 1st single moment = ', XRIGHTC 890 WRITE(LUPRI,'(1X,A,F16.8,A1,F16.8,A1)') 891 * ' Dipole Strength (a.u.) = ',FORZADIP,'(',SQRFORZA,')' 892 WRITE(LUPRI,'(/1X,A,A,A,A,3(/1X,A,F16.8))') 893 * 'For A, B operators ', LABELA,', ', LABELB, 894 * ' Left 2nd residue AB = ', XLEFTAB, 895 * ' Right 2nd residue AB = ', XRIGHTAB, 896 * ' S^of_AB,AB = ', XLEFTAB*XRIGHTAB 897 CALL FLSHFO(LUPRI) 898 END IF 899* 900* ----------------------------------------------------------------------* 901* combine various terms to get B term contribution 902* ----------------------------------------------------------------------* 903* 904 BTERM_S = D05 * (XMAB_MC + XMC_MAB) 905 BTERM_A = D05 * (XMAB_MC - XMC_MAB) 906* 907* ----------------------------------------------------------------------* 908* Write output 909* Note: final result divided by 2 because of (mu = -0.5 L) 910* ----------------------------------------------------------------------* 911* 912 WRITE(LUPRI,'(1X,58("-"))') 913 WRITE(LUPRI,'(/1x,a,f9.5,a,i1)') 914 & 'For transition |o> -> |f(',FREQEX,')>, of symm. ', ISYMC 915 WRITE(LUPRI,'(3(/1x,a,a,a,a1,f9.5,a,i1))') 916 & ' A oper.: ',LABELA,CRLXA,'(',-FREQEX,'), symm. ',ISYMA, 917 & ' B oper.: ',LABELB,CRLXB,'(',ZERO ,'), symm. ',ISYMB, 918 & ' C oper.: ',LABELC,CRLXC,'(',FREQEX ,'), symm. ',ISYMC 919C WRITE(LUPRI,'(3(/1X,A,F16.8))') 920C & ' S^of_C,C = (Dipole strength (au)) = ', FORZADIP, 921C & ' S^of_AB,C = M^AB_of(0.0) x M^C_fo = ', XMAB_MC, 922C & ' S^of_C,AB = M^C*_of x M^AB*_fo(0.0) = ', XMC_MAB 923 WRITE(LUPRI,'(3(/1X,A,F16.8))') 924 & ' Dipole strength (au) (C oper) = ',FORZADIP, 925 & ' M^AB_of(0.0) x M^C_fo = ',XMAB_MC, 926 & ' M^C*_of x M^AB*_fo(0.0) = ',XMC_MAB 927 IF (IPRINT.GT.5) THEN 928 WRITE(LUPRI,'(/1X,a5,A,F16.6,A,F16.8,A)') 929 & MODPRI, 930 & ' B term contribution (au): ', D05*BTERM_A, ' (antisym) ', 931 & D05*BTERM_S, ' (symm) ' 932 WRITE(LUPRI,'(1X,69("-"))') 933 ELSE 934 WRITE(LUPRI,'(/1X,A5,A,F16.8,A)') 935 & MODPRI, 936 & ' B term contribution (au): ', D05*BTERM_A, ' (antisym) ' 937 WRITE(LUPRI,'(1X,58("-"))') 938 END IF 939* ---------------------------------------------------------------------* 940 END IF 941 END DO 942 943 ENDIF 944 100 CONTINUE 945 946 RETURN 947 END 948*=====================================================================* 949* END OF SUBROUTINE CC_MCDCAL * 950*=====================================================================* 951