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 cc_hyppol */ 20*=====================================================================* 21 SUBROUTINE CC_HYPPOL(WORK,LWORK) 22*---------------------------------------------------------------------* 23* 24* Purpose: direct calculation of frequency-dependent 25* first hyperpolarizabilities for the Couple Cluster models 26* 27* CCS, CC2, CCSD 28* 29* Written by Christof Haettig summer/fall 1996. 30* Dispersion coefficients, october 1997, Christof Haettig. 31* Orbital relax. for one operator, june 1999, Christof Haettig. 32* 33*=====================================================================* 34 USE PELIB_INTERFACE, ONLY: USE_PELIB 35#if defined (IMPLICIT_NONE) 36 IMPLICIT NONE 37#else 38# include "implicit.h" 39#endif 40#include "priunit.h" 41#include "ccsdinp.h" 42#include "ccorb.h" 43#include "ccsdsym.h" 44#include "ccqrinf.h" 45#include "ccr1rsp.h" 46#include "ccrc1rsp.h" 47#include "ccroper.h" 48#include "cclists.h" 49#include "second.h" 50#include "dummy.h" 51#include "ccsections.h" 52#include "ccslvinf.h" 53 54* local parameters: 55 CHARACTER*(19) MSGDBG 56 PARAMETER (MSGDBG = '[debug] CC_HYPPOL> ') 57 LOGICAL LOCDBG 58 PARAMETER (LOCDBG = .FALSE.) 59 60#if defined (SYS_CRAY) 61 REAL ZERO 62#else 63 DOUBLE PRECISION ZERO 64#endif 65 PARAMETER (ZERO = 0.0d0) 66 67 CHARACTER*10 MODEL 68 CHARACTER*8 LABSOP 69 LOGICAL LADD 70 LOGICAL LORXA, LORXB, LORXC, LPDBSA, LPDBSB, LPDBSC 71 INTEGER LWORK 72 INTEGER ISYMA, ITAMPA, IZETVA, IOPERA, IKAPA, IR2AB 73 INTEGER ISYMB, ITAMPB, IZETVB, IOPERB, IKAPB, IR2AC 74 INTEGER ISYMC, ITAMPC, IZETVC, IOPERC, IKAPC, IR2BC, IOPSOP 75 INTEGER KHYPPOL, NBHYPPOL, KEND0, LEND0, KEND1, LEND1 76 INTEGER KEXPCOF, NBEXPCOF, NBINOM, KBINOM, KDSPCF 77 INTEGER KDSPCFA, KSHGCFA, KORCOFA, KEOPECA, KABCDE 78 INTEGER KSHGCF, KORCOF, KEOPEC, ISGNSOP 79 INTEGER MXTRAN, MXVEC, ITRAN, IVEC, ISIGN, IORDER, ISYSOP, INUM 80 INTEGER IOPT, IDX, IFREQ, ISYM, ISYML, ISYML1, ISYML2, IOPER 81 INTEGER KGTRAN, KGDOTS, KGCONS, NGTRAN, MXGTRAN, MXGDOTS 82 INTEGER KFTRAN, KFDOTS, KFCONS, NFTRAN, MXFTRAN, MXFDOTS 83 INTEGER KBTRAN, KBDOTS, KBCONS, NBTRAN, MXBTRAN, MXBDOTS 84 INTEGER KKTRAN, KKDOTS, KKCONS, NKTRAN, MXKTRAN, MXKDOTS 85 INTEGER KFATRAN, KFADOTS, KFACONS, NFATRAN, MXFATRAN, MXFADOTS 86 INTEGER KEATRAN, KEADOTS, KEACONS, NEATRAN, MXEATRAN, MXEADOTS 87 INTEGER KR2TRAN, KR2DOTS, KR2CONS, NR2TRAN, MXR2TRAN, MXR2DOTS 88 INTEGER KXETRAN, KEDOTS, KECONS, NXETRAN, MXXETRAN, MXXEDOTS 89 INTEGER KAATRAN, KAADOTS, KAACONS, NAATRAN, MXAATRAN, MXAADOTS 90 INTEGER KXDOTS, KXCONS 91 92#if defined (SYS_CRAY) 93 REAL WORK(LWORK) 94 REAL GCON, FACON1,FACON2,FACON3, FCON1,FCON2,FCON3 95 REAL EACON1, EACON2, EACON3, EACON4, EACON5, EACON6 96 REAL BCON1, BCON2, BCON3, R2CON1, R2CON2, R2CON3 97 REAL SLVKCON1, SLVKCON2, SLVKCON3 98 REAL X2CON1, X2CON2, X2CON3, E2CON1, E2CON2, E2CON3 99 REAL SIGN 100 REAL TIMTOT, TIM0, TIM1, TIMG, TIMF, TIMB, TIMFA, TIMEA, TIM2 101 REAL TIMR2, TIMXE, TIMAA, TIMK 102#else 103 DOUBLE PRECISION WORK(LWORK) 104 DOUBLE PRECISION GCON, FACON1,FACON2,FACON3, FCON1,FCON2,FCON3 105 DOUBLE PRECISION EACON1, EACON2, EACON3, EACON4, EACON5, EACON6 106 DOUBLE PRECISION BCON1, BCON2, BCON3, R2CON1, R2CON2, R2CON3 107 DOUBLE PRECISION SLVKCON1, SLVKCON2, SLVKCON3 108 DOUBLE PRECISION X2CON1, X2CON2, X2CON3, E2CON1, E2CON2, E2CON3 109 DOUBLE PRECISION SIGN 110 DOUBLE PRECISION TIMTOT,TIM0,TIM1,TIMG,TIMF,TIMB,TIMFA,TIMEA,TIM2 111 DOUBLE PRECISION TIMR2, TIMXE, TIMAA, TIMK 112#endif 113 114* external functions 115 INTEGER IR1TAMP 116 INTEGER IR1KAPPA 117 INTEGER IL1ZETA 118 INTEGER IR2TAMP 119 INTEGER IROPER 120 121*---------------------------------------------------------------------* 122* print header for hyperpolarizability section 123*---------------------------------------------------------------------* 124 WRITE (LUPRI,'(7(/1X,2A),/)') 125 & '************************************', 126 & '*******************************', 127 & '* ', 128 & ' *', 129 & '*-------- OUTPUT FROM COUPLED CLUSTER Q', 130 & 'UADRATIC RESPONSE ---------*', 131 & '* ', 132 & ' *', 133 & '*-------- CALCULATION OF CC HYPER', 134 & 'POLARIZABILITIES ---------*', 135 & '* ', 136 & ' *', 137 & '************************************', 138 & '*******************************' 139 140*---------------------------------------------------------------------* 141 IF (.NOT. (CCS .OR. CC2 .OR. CCSD .OR. CC3) ) THEN 142 CALL QUIT('CC_HYPPOL called for unknown Coupled Cluster.') 143 END IF 144 145* print some debug/info output 146 IF (IPRINT .GT. 10) WRITE(LUPRI,*) 'CC_HYPPOL Workspace:',LWORK 147 148 TIM0 = SECOND() 149 150 call FLSHFO(LUPRI) 151 152*---------------------------------------------------------------------* 153* allocate & initialize work space for hyperpolarizabilities 154*---------------------------------------------------------------------* 155 NBHYPPOL = NQROPER * NQRFREQ 156 157 MXTRAN = 2 * NLRTLBL * NLRTLBL 158 MXVEC = NLRTLBL 159 160 MXGTRAN = MXDIM_GTRAN * MXTRAN 161 MXFTRAN = MXDIM_FTRAN * MXTRAN 162 MXBTRAN = MXDIM_BTRAN * MXTRAN 163 MXFATRAN = MXDIM_FATRAN * MXTRAN 164 MXAATRAN = MXDIM_AATRAN * MXTRAN 165 MXEATRAN = MXDIM_XEVEC * MXTRAN 166 MXR2TRAN = 1 * MXTRAN 167 MXXETRAN = MXDIM_XEVEC * MXTRAN 168 MXKTRAN = MXDIM_KTRAN * MXTRAN 169 170 MXGDOTS = MXVEC * MXTRAN 171 MXFDOTS = MXVEC * MXTRAN 172 MXBDOTS = MXVEC * MXTRAN 173 MXFADOTS = MXVEC * MXTRAN 174 MXEADOTS = MXVEC * MXTRAN 175 MXR2DOTS = MXVEC * MXTRAN 176 MXXEDOTS = MXVEC * MXTRAN 177 MXAADOTS = MXVEC * MXTRAN 178 MXKDOTS = MXVEC * MXTRAN 179 180 KHYPPOL = 1 181 KGTRAN = KHYPPOL + 2 * NBHYPPOL 182 KGDOTS = KGTRAN + MXGTRAN 183 KGCONS = KGDOTS + MXGDOTS 184 KFTRAN = KGCONS + MXGDOTS 185 KFDOTS = KFTRAN + MXFTRAN 186 KFCONS = KFDOTS + MXFDOTS 187 KBTRAN = KFCONS + MXFDOTS 188 KBDOTS = KBTRAN + MXBTRAN 189 KBCONS = KBDOTS + MXBDOTS 190 KKTRAN = KBCONS + MXBDOTS 191 KKDOTS = KKTRAN + MXKTRAN 192 KKCONS = KKDOTS + MXKDOTS 193 KFATRAN = KKCONS + MXKDOTS 194 KFADOTS = KFATRAN + MXFATRAN 195 KFACONS = KFADOTS + MXFADOTS 196 KEATRAN = KFACONS + MXFADOTS 197 KEADOTS = KEATRAN + MXEATRAN 198 KEACONS = KEADOTS + MXEADOTS 199 KR2TRAN = KEACONS + MXEADOTS 200 KR2DOTS = KR2TRAN + MXR2TRAN 201 KR2CONS = KR2DOTS + MXR2DOTS 202 KXETRAN = KR2CONS + MXR2DOTS 203 KXDOTS = KXETRAN + MXXETRAN 204 KXCONS = KXDOTS + MXXEDOTS 205 KEDOTS = KXCONS + MXXEDOTS 206 KECONS = KEDOTS + MXXEDOTS 207 KAATRAN = KECONS + MXXEDOTS 208 KAADOTS = KAATRAN + MXAATRAN 209 KAACONS = KAADOTS + MXAADOTS 210 KEND0 = KAACONS + MXAADOTS 211 LEND0 = LWORK - KEND0 212 213 IF (LEND0 .LT. 0) THEN 214 CALL QUIT('Insufficient memory in CC_HYPPOL. (1)') 215 END IF 216 217 CALL DZERO(WORK(KHYPPOL),2*NBHYPPOL) 218 219*---------------------------------------------------------------------* 220* set up lists for G, F and F{A} transformations and ETA{O} vectors: 221*---------------------------------------------------------------------* 222 CALL CCQR_SETUP(MXTRAN, MXVEC, 223 & WORK(KGTRAN), WORK(KGDOTS), NGTRAN, 224 & WORK(KFTRAN), WORK(KFDOTS), NFTRAN, 225 & WORK(KBTRAN), WORK(KBDOTS), NBTRAN, 226 & WORK(KKTRAN), WORK(KKDOTS), NKTRAN, 227 & WORK(KFATRAN),WORK(KFADOTS),NFATRAN, 228 & WORK(KAATRAN),WORK(KAADOTS),NAATRAN, 229 & WORK(KEATRAN),WORK(KEADOTS),NEATRAN, 230 & WORK(KR2TRAN),WORK(KR2DOTS),NR2TRAN, 231 & WORK(KXETRAN),WORK(KXDOTS),WORK(KEDOTS),NXETRAN, 232 & WORK(KEND0), LEND0) 233 234*---------------------------------------------------------------------* 235* calculate G matrix contributions: 236*---------------------------------------------------------------------* 237 238 TIM1 = SECOND() 239 240 IOPT = 5 241 CALL CC_GMATRIX('L0 ','R1 ','R1 ','R1 ',NGTRAN, MXVEC, 242 & WORK(KGTRAN),WORK(KGDOTS),WORK(KGCONS), 243 & WORK(KEND0), LEND0, IOPT ) 244 245 246 TIMG = SECOND() - TIM1 247 248 WRITE (LUPRI,'(/A,I5,A,F12.2," seconds.")') 249 & ' Time used for',NGTRAN,' G matrix transformations:',TIMG 250 CALL FLSHFO(LUPRI) 251 252 253*---------------------------------------------------------------------* 254* calculate B matrix contributions: 255*---------------------------------------------------------------------* 256 IF ( USEBTRAN .AND. (.NOT.USE_R2) ) THEN 257 TIM1 = SECOND() 258 259 IOPT = 5 260 CALL CC_BMATRIX(WORK(KBTRAN),NBTRAN,'R1 ','R1 ',IOPT,'L1 ', 261 & WORK(KBDOTS),WORK(KBCONS),MXVEC, 262 & .FALSE.,WORK(KEND0),LEND0) 263 264 TIMB = SECOND() - TIM1 265 266 WRITE (LUPRI,'(/A,I5,A,F12.2," seconds.")') 267 & ' Time used for',NBTRAN,' B matrix transformations:',TIMB 268 CALL FLSHFO(LUPRI) 269 270 ELSE 271 TIMB = 0.0d0 272 END IF 273*---------------------------------------------------------------------* 274* calculate K matrix contributions: 275*---------------------------------------------------------------------* 276 IF (CCSLV.OR.USE_PELIB()) THEN 277 TIM1 = SECOND() 278 279 IOPT = 5 280 CALL CC_KMATRIX(WORK(KKTRAN),NKTRAN,'L1 ','L1 ',IOPT,'R1 ', 281 & WORK(KKDOTS),WORK(KKCONS),MXVEC,WORK(KEND0),LEND0) 282 283 TIMK = SECOND() - TIM1 284 285 WRITE (LUPRI,'(/A,I5,A,F12.2," seconds.")') 286 & ' Time used for',NKTRAN,' K matrix transformations:',TIMK 287 CALL FLSHFO(LUPRI) 288 289 ELSE 290 TIMK = 0.0d0 291 END IF 292 293*---------------------------------------------------------------------* 294* calculate F matrix contributions: 295*---------------------------------------------------------------------* 296 IF ( (.NOT.USEBTRAN) .AND. (.NOT.USE_R2) ) THEN 297 TIM1 = SECOND() 298 299 IOPT = 5 300 CALL CC_FMATRIX(WORK(KFTRAN),NFTRAN,'L1 ','R1 ',IOPT,'R1 ', 301 & WORK(KFDOTS),WORK(KFCONS),MXVEC, 302 & WORK(KEND0), LEND0) 303 304 TIMF = SECOND() - TIM1 305 306 WRITE (LUPRI,'(/A,I5,A,F12.2," seconds.")') 307 & ' Time used for',NFTRAN,' F matrix transformations:',TIMF 308 CALL FLSHFO(LUPRI) 309 310 ELSE 311 TIMF = 0.0d0 312 END IF 313 314*---------------------------------------------------------------------* 315* calculate F{O} matrix contributions: 316*---------------------------------------------------------------------* 317 TIM1 = SECOND() 318 319 CALL CCQR_FADRV('L0 ','o1 ','R1 ','R1 ',NFATRAN, MXVEC, 320 & WORK(KFATRAN),WORK(KFADOTS),WORK(KFACONS), 321 & WORK(KEND0), LEND0, 'DOTP' ) 322 323 TIMFA = SECOND() - TIM1 324 325 WRITE (LUPRI,'(/A,I5,A,F12.2," seconds.")') 326 & ' Time used for',NFATRAN,' F{O} matrix transform.: ',TIMFA 327 CALL FLSHFO(LUPRI) 328 329*---------------------------------------------------------------------* 330* calculate ETA{O} vector contributions: 331*---------------------------------------------------------------------* 332 CALL DZERO(WORK(KEACONS),MXEADOTS) 333 334 IF ((.NOT. USE_R2) .AND. (.NOT. USE_AAMAT)) THEN 335 TIM1 = SECOND() 336 337 IOPT = 5 338 IORDER = 1 339 CALL CC_XIETA( WORK(KEATRAN), NEATRAN, IOPT, IORDER, 'L1 ', 340 & '---',IDUMMY, DUMMY, 341 & 'R1 ',WORK(KEADOTS),WORK(KEACONS), 342 & .FALSE.,MXVEC, WORK(KEND0), LEND0 ) 343 344 TIMEA = SECOND() - TIM1 345 WRITE (LUPRI,'(/A,I5,A,F12.2," seconds.")') 346 & ' Time used for', 347 & NEATRAN,' ETA{O} vector calculat.: ', SECOND() - TIM1 348 CALL FLSHFO(LUPRI) 349 350 ELSE 351 TIMEA = 0.0d0 352 END IF 353 354*---------------------------------------------------------------------* 355* calculate A{O} matrix contributions: 356*---------------------------------------------------------------------* 357 CALL DZERO(WORK(KAACONS),MXAADOTS) 358 359 IF ((.NOT. USE_R2) .AND. USE_AAMAT) THEN 360 TIM1 = SECOND() 361 362 IOPT = 5 363 CALL CC_AAMAT(WORK(KAATRAN),NAATRAN,'o1 ','R1 ',IOPT,'L1 ', 364 & WORK(KAADOTS),WORK(KAACONS), 365 & MXVEC, WORK(KEND0), LEND0 ) 366 367 TIMAA = SECOND() - TIM1 368 WRITE (LUPRI,'(/A,I5,A,F12.2," seconds.")') 369 & ' Time used for', 370 & NAATRAN,' A{O} matrix transform. : ', SECOND() - TIM1 371 CALL FLSHFO(LUPRI) 372 373 ELSE 374 TIMAA = 0.0d0 375 END IF 376 377*---------------------------------------------------------------------* 378* calculate X1 x R2 dots products: 379*---------------------------------------------------------------------* 380 IF ( USE_R2 ) THEN 381 TIM1 = SECOND() 382 383C CALL CC_DOTDRV('O2 ','L1 ',NR2TRAN,MXVEC, 384C & WORK(KR2TRAN), WORK(KR2DOTS), WORK(KR2CONS), 385C & WORK(KEND0), LEND0 ) 386 CALL CC_DOTDRV('R2 ','X1B',NR2TRAN,MXVEC, 387 & WORK(KR2TRAN), WORK(KR2DOTS), WORK(KR2CONS), 388 & WORK(KEND0), LEND0 ) 389 390 TIMR2 = SECOND() - TIM1 391 WRITE (LUPRI,'(/A,I5,A,F12.2," seconds.")') 392 & ' Time used for', 393 & NR2TRAN,' X1 x R2 dot products: ', SECOND() - TIM1 394 CALL FLSHFO(LUPRI) 395 396 ELSE 397 TIMR2 = 0.0d0 398 END IF 399 400*---------------------------------------------------------------------* 401* calculate ETA{O2} x R1 and Xksi{O2} x L1 dots products: 402*---------------------------------------------------------------------* 403 TIM1 = SECOND() 404 405 CALL DZERO(WORK(KXCONS),MXXEDOTS) 406 CALL DZERO(WORK(KECONS),MXXEDOTS) 407 408 IOPT = 5 409 IORDER = 2 410 CALL CC_XIETA( WORK(KXETRAN), NXETRAN, IOPT, IORDER, 'L0 ', 411 & 'L1 ',WORK(KXDOTS),WORK(KXCONS), 412 & 'R1 ',WORK(KEDOTS),WORK(KECONS), 413 & .FALSE.,MXVEC, WORK(KEND0), LEND0 ) 414 415 TIMXE = SECOND() - TIM1 416 WRITE (LUPRI,'(/A,I5,A,F12.2," seconds.")') ' Time used for', 417 & NXETRAN,' Xksi{O} and Eta{O2} vector calculat.: ', SECOND()-TIM1 418 CALL FLSHFO(LUPRI) 419 420*=====================================================================* 421* loop over triples of operator labels and over frequencies and 422* collect all contributions to the final hyperpolarizabilities: 423*=====================================================================* 424 DO IOPER = 1, NQROPER 425 426 IOPERA = IAQROP(IOPER) 427 IOPERB = IBQROP(IOPER) 428 IOPERC = ICQROP(IOPER) 429 430 ISYMB = ISYOPR(IOPERB) 431 ISYMC = ISYOPR(IOPERC) 432 ISYMA = ISYOPR(IOPERA) 433 434 435 IF (MULD2H(ISYMB,ISYMC) .EQ. ISYMA) THEN 436 437 LORXA = LAQLRX(IOPER) 438 LORXB = LBQLRX(IOPER) 439 LORXC = LCQLRX(IOPER) 440 441 LPDBSA = LPDBSOP(IOPERA) 442 LPDBSB = LPDBSOP(IOPERB) 443 LPDBSC = LPDBSOP(IOPERC) 444 445 DO IFREQ = 1, NQRFREQ 446 447 DO ISIGN = +1, -1, -2 448 SIGN = DBLE(ISIGN) 449 450 ITAMPA = IR1TAMP(LBLOPR(IOPERA),LORXA,SIGN*AQRFR(IFREQ),ISYML) 451 ITAMPB = IR1TAMP(LBLOPR(IOPERB),LORXB,SIGN*BQRFR(IFREQ),ISYML) 452 ITAMPC = IR1TAMP(LBLOPR(IOPERC),LORXC,SIGN*CQRFR(IFREQ),ISYML) 453 454 IZETVA = IL1ZETA(LBLOPR(IOPERA),LORXA,SIGN*AQRFR(IFREQ),ISYML) 455 IZETVB = IL1ZETA(LBLOPR(IOPERB),LORXB,SIGN*BQRFR(IFREQ),ISYML) 456 IZETVC = IL1ZETA(LBLOPR(IOPERC),LORXC,SIGN*CQRFR(IFREQ),ISYML) 457 458 IKAPA = 0 459 IKAPB = 0 460 IKAPC = 0 461 IF(LORXA)IKAPA=IR1KAPPA(LBLOPR(IOPERA),SIGN*AQRFR(IFREQ),ISYML) 462 IF(LORXB)IKAPB=IR1KAPPA(LBLOPR(IOPERB),SIGN*BQRFR(IFREQ),ISYML) 463 IF(LORXC)IKAPC=IR1KAPPA(LBLOPR(IOPERC),SIGN*CQRFR(IFREQ),ISYML) 464 465 IF ( USE_R2 ) THEN 466 IR2AB=IR2TAMP(LBLOPR(IOPERA),.FALSE.,SIGN*AQRFR(IFREQ),ISYML1, 467 & LBLOPR(IOPERB),.FALSE.,SIGN*BQRFR(IFREQ),ISYML2) 468 IR2AC=IR2TAMP(LBLOPR(IOPERA),.FALSE.,SIGN*AQRFR(IFREQ),ISYML1, 469 & LBLOPR(IOPERC),.FALSE.,SIGN*CQRFR(IFREQ),ISYML2) 470 IR2BC=IR2TAMP(LBLOPR(IOPERB),.FALSE.,SIGN*BQRFR(IFREQ),ISYML1, 471 & LBLOPR(IOPERC),.FALSE.,SIGN*CQRFR(IFREQ),ISYML2) 472 END IF 473 474*---------------------------------------------------------------------* 475 476 CALL CCQR_SETG(WORK(KGTRAN),WORK(KGDOTS),NGTRAN,MXVEC, 477 & 0,ITAMPA,ITAMPB,ITAMPC,ITRAN,IVEC) 478 GCON = WORK(KGCONS-1 + (ITRAN-1)*MXVEC + IVEC) 479 480*---------------------------------------------------------------------* 481 482 CALL CC_SETFA(WORK(KFATRAN),WORK(KFADOTS),NFATRAN,MXVEC, 483 & 0,IOPERA,IKAPA,ITAMPB,ITAMPC,ITRAN,IVEC) 484 FACON1 = WORK(KFACONS-1 + (ITRAN-1)*MXVEC + IVEC) 485 486 CALL CC_SETFA(WORK(KFATRAN),WORK(KFADOTS),NFATRAN,MXVEC, 487 & 0,IOPERB,IKAPB,ITAMPC,ITAMPA,ITRAN,IVEC) 488 FACON2 = WORK(KFACONS-1 + (ITRAN-1)*MXVEC + IVEC) 489 490 CALL CC_SETFA(WORK(KFATRAN),WORK(KFADOTS),NFATRAN,MXVEC, 491 & 0,IOPERC,IKAPC,ITAMPA,ITAMPB,ITRAN,IVEC) 492 FACON3 = WORK(KFACONS-1 + (ITRAN-1)*MXVEC + IVEC) 493 494C--------------------------------------------------------------------* 495C K-Matrix contribution: 496C--------------------------------------------------------------------* 497C 498 IF (CCSLV.OR.USE_PELIB()) THEN 499 CALL CC_SETB11(WORK(KKTRAN),WORK(KKDOTS),NKTRAN,MXVEC, 500 & ITAMPA,IZETVB,IZETVC,ITRAN,IVEC) 501 502 SLVKCON1 = WORK(KKCONS-1 + (ITRAN-1)*MXVEC + IVEC) 503 504 CALL CC_SETB11(WORK(KKTRAN),WORK(KKDOTS),NKTRAN,MXVEC, 505 & ITAMPB,IZETVC,IZETVA,ITRAN,IVEC) 506 SLVKCON2 = WORK(KKCONS-1 + (ITRAN-1)*MXVEC + IVEC) 507 508 CALL CC_SETB11(WORK(KKTRAN),WORK(KKDOTS),NKTRAN,MXVEC, 509 & ITAMPC,IZETVA,IZETVB,ITRAN,IVEC) 510 SLVKCON3 = WORK(KKCONS-1 + (ITRAN-1)*MXVEC + IVEC) 511 ELSE 512 SLVKCON1 = 0.0d0 513 SLVKCON2 = 0.0d0 514 SLVKCON3 = 0.0d0 515 END IF 516C 517C---------------------------------------------------------------------* 518C 519 IF ( USEBTRAN .AND. (.NOT.USE_R2) ) THEN 520 CALL CC_SETB11(WORK(KBTRAN),WORK(KBDOTS),NBTRAN,MXVEC, 521 & IZETVA,ITAMPB,ITAMPC,ITRAN,IVEC) 522 BCON1 = WORK(KBCONS-1 + (ITRAN-1)*MXVEC + IVEC) 523 524 CALL CC_SETB11(WORK(KBTRAN),WORK(KBDOTS),NBTRAN,MXVEC, 525 & IZETVB,ITAMPC,ITAMPA,ITRAN,IVEC) 526 BCON2 = WORK(KBCONS-1 + (ITRAN-1)*MXVEC + IVEC) 527 528 CALL CC_SETB11(WORK(KBTRAN),WORK(KBDOTS),NBTRAN,MXVEC, 529 & IZETVC,ITAMPA,ITAMPB,ITRAN,IVEC) 530 BCON3 = WORK(KBCONS-1 + (ITRAN-1)*MXVEC + IVEC) 531 532 FCON1 = 0.0d0 533 FCON2 = 0.0d0 534 FCON3 = 0.0d0 535 R2CON1 = 0.0d0 536 R2CON2 = 0.0d0 537 R2CON3 = 0.0d0 538 539 ELSE IF ( (.NOT.USEBTRAN) .AND. (.NOT.USE_R2) ) THEN 540 541 CALL CCQR_SETF(WORK(KFTRAN),WORK(KFDOTS),NFTRAN,MXVEC, 542 & IZETVA,ITAMPB,ITAMPC,ITRAN,IVEC) 543 FCON1 = WORK(KFCONS-1 + (ITRAN-1)*MXVEC + IVEC) 544 545 CALL CCQR_SETF(WORK(KFTRAN),WORK(KFDOTS),NFTRAN,MXVEC, 546 & IZETVB,ITAMPC,ITAMPA,ITRAN,IVEC) 547 FCON2 = WORK(KFCONS-1 + (ITRAN-1)*MXVEC + IVEC) 548 549 CALL CCQR_SETF(WORK(KFTRAN),WORK(KFDOTS),NFTRAN,MXVEC, 550 & IZETVC,ITAMPA,ITAMPB,ITRAN,IVEC) 551 FCON3 = WORK(KFCONS-1 + (ITRAN-1)*MXVEC + IVEC) 552 553 BCON1 = 0.0d0 554 BCON2 = 0.0d0 555 BCON3 = 0.0d0 556 R2CON1 = 0.0d0 557 R2CON2 = 0.0d0 558 R2CON3 = 0.0d0 559 560 ELSE IF ( USE_R2 ) THEN 561 562 CALL CC_SETDOT(WORK(KR2TRAN),WORK(KR2DOTS),NR2TRAN,MXVEC, 563 & IR2AB,IZETVC,ITRAN,IVEC) 564 R2CON1 = WORK(KR2CONS-1 + (ITRAN-1)*MXVEC + IVEC) 565 566 CALL CC_SETDOT(WORK(KR2TRAN),WORK(KR2DOTS),NR2TRAN,MXVEC, 567 & IR2AC,IZETVB,ITRAN,IVEC) 568 R2CON2 = WORK(KR2CONS-1 + (ITRAN-1)*MXVEC + IVEC) 569 570 CALL CC_SETDOT(WORK(KR2TRAN),WORK(KR2DOTS),NR2TRAN,MXVEC, 571 & IR2BC,IZETVA,ITRAN,IVEC) 572 R2CON3 = WORK(KR2CONS-1 + (ITRAN-1)*MXVEC + IVEC) 573 574 BCON1 = 0.0d0 575 BCON2 = 0.0d0 576 BCON3 = 0.0d0 577 FCON1 = 0.0d0 578 FCON2 = 0.0d0 579 FCON3 = 0.0d0 580 ELSE 581 CALL QUIT('Error in CC_HYPPOL.') 582 END IF 583 584*---------------------------------------------------------------------* 585 586 IF ( .NOT. USE_R2 ) THEN 587 IF (USE_AAMAT) THEN 588 IF (LOCDBG) WRITE(LUPRI,*) 'Use A{O} matrix transformation.' 589 IF (IKAPA.NE.0 .OR. IKAPB.NE.0 .OR. IKAPC.NE.0) THEN 590 CALL QUIT('A{O} matrix transformation not yet relaxed.') 591 END IF 592 593 CALL CC_SETAA(WORK(KAATRAN),WORK(KAADOTS),MXTRAN,MXVEC, 594 & IZETVA,IOPERB,ITAMPC,ITRAN,IVEC) 595 EACON1 = WORK(KAACONS-1 + (ITRAN-1)*MXVEC + IVEC) 596 597 CALL CC_SETAA(WORK(KAATRAN),WORK(KAADOTS),MXTRAN,MXVEC, 598 & IZETVA,IOPERC,ITAMPB,ITRAN,IVEC) 599 EACON2 = WORK(KAACONS-1 + (ITRAN-1)*MXVEC + IVEC) 600 601 CALL CC_SETAA(WORK(KAATRAN),WORK(KAADOTS),MXTRAN,MXVEC, 602 & IZETVB,IOPERC,ITAMPA,ITRAN,IVEC) 603 EACON3 = WORK(KAACONS-1 + (ITRAN-1)*MXVEC + IVEC) 604 605 CALL CC_SETAA(WORK(KAATRAN),WORK(KAADOTS),MXTRAN,MXVEC, 606 & IZETVB,IOPERA,ITAMPC,ITRAN,IVEC) 607 EACON4 = WORK(KAACONS-1 + (ITRAN-1)*MXVEC + IVEC) 608 609 CALL CC_SETAA(WORK(KAATRAN),WORK(KAADOTS),MXTRAN,MXVEC, 610 & IZETVC,IOPERA,ITAMPB,ITRAN,IVEC) 611 EACON5 = WORK(KAACONS-1 + (ITRAN-1)*MXVEC + IVEC) 612 613 CALL CC_SETAA(WORK(KAATRAN),WORK(KAADOTS),MXTRAN,MXVEC, 614 & IZETVC,IOPERB,ITAMPA,ITRAN,IVEC) 615 EACON6 = WORK(KAACONS-1 + (ITRAN-1)*MXVEC + IVEC) 616 ELSE 617 CALL CC_SETXE('Eta',WORK(KEATRAN),WORK(KEADOTS),NEATRAN,MXVEC, 618 & IZETVA,IOPERB,IKAPB,0,0,0,ITAMPC,ITRAN,IVEC) 619 EACON1 = WORK(KEACONS-1 + (ITRAN-1)*MXVEC + IVEC) 620 621 CALL CC_SETXE('Eta',WORK(KEATRAN),WORK(KEADOTS),NEATRAN,MXVEC, 622 & IZETVA,IOPERC,IKAPC,0,0,0,ITAMPB,ITRAN,IVEC) 623 EACON2 = WORK(KEACONS-1 + (ITRAN-1)*MXVEC + IVEC) 624 625 CALL CC_SETXE('Eta',WORK(KEATRAN),WORK(KEADOTS),NEATRAN,MXVEC, 626 & IZETVB,IOPERC,IKAPC,0,0,0,ITAMPA,ITRAN,IVEC) 627 EACON3 = WORK(KEACONS-1 + (ITRAN-1)*MXVEC + IVEC) 628 629 CALL CC_SETXE('Eta',WORK(KEATRAN),WORK(KEADOTS),NEATRAN,MXVEC, 630 & IZETVB,IOPERA,IKAPA,0,0,0,ITAMPC,ITRAN,IVEC) 631 EACON4 = WORK(KEACONS-1 + (ITRAN-1)*MXVEC + IVEC) 632 633 CALL CC_SETXE('Eta',WORK(KEATRAN),WORK(KEADOTS),NEATRAN,MXVEC, 634 & IZETVC,IOPERA,IKAPA,0,0,0,ITAMPB,ITRAN,IVEC) 635 EACON5 = WORK(KEACONS-1 + (ITRAN-1)*MXVEC + IVEC) 636 637 CALL CC_SETXE('Eta',WORK(KEATRAN),WORK(KEADOTS),NEATRAN,MXVEC, 638 & IZETVC,IOPERB,IKAPB,0,0,0,ITAMPA,ITRAN,IVEC) 639 EACON6 = WORK(KEACONS-1 + (ITRAN-1)*MXVEC + IVEC) 640 END IF 641 ELSE 642 643 EACON1 = 0.0d0 644 EACON2 = 0.0d0 645 EACON3 = 0.0d0 646 EACON4 = 0.0d0 647 EACON5 = 0.0d0 648 EACON6 = 0.0d0 649 650 END IF 651 652 653*---------------------------------------------------------------------* 654 IF (LPDBSA.OR.LORXA .OR. LPDBSB.OR.LORXB) THEN 655 CALL CC_FIND_SO_OP(LBLOPR(IOPERA),LBLOPR(IOPERB), 656 & LABSOP,ISYSOP,ISGNSOP,INUM,WORK,LWORK) 657 IOPSOP = IROPER(LABSOP,ISYSOP) 658 659 CALL CC_SETXE('Xi ',WORK(KXETRAN),WORK(KXDOTS),MXTRAN,MXVEC, 660 & 0,IOPSOP,IKAPA,IKAPB,0,0,IZETVC,ITRAN,IVEC) 661 X2CON1 = WORK(KXCONS-1 + (ITRAN-1)*MXVEC + IVEC) 662 663 CALL CC_SETXE('Eta',WORK(KXETRAN),WORK(KEDOTS),MXTRAN,MXVEC, 664 & 0,IOPSOP,IKAPA,IKAPB,0,0,ITAMPC,ITRAN,IVEC) 665 E2CON1 = WORK(KECONS-1 + (ITRAN-1)*MXVEC + IVEC) 666 ELSE 667 X2CON1 = 0.0d0 668 E2CON1 = 0.0d0 669 END IF 670 671 IF (LPDBSA.OR.LORXA .OR. LPDBSC.OR.LORXC) THEN 672 CALL CC_FIND_SO_OP(LBLOPR(IOPERA),LBLOPR(IOPERC), 673 & LABSOP,ISYSOP,ISGNSOP,INUM,WORK,LWORK) 674 IOPSOP = IROPER(LABSOP,ISYSOP) 675 676 CALL CC_SETXE('Xi ',WORK(KXETRAN),WORK(KXDOTS),MXTRAN,MXVEC, 677 & 0,IOPSOP,IKAPA,IKAPC,0,0,IZETVB,ITRAN,IVEC) 678 X2CON2 = WORK(KXCONS-1 + (ITRAN-1)*MXVEC + IVEC) 679 680 CALL CC_SETXE('Eta',WORK(KXETRAN),WORK(KEDOTS),MXTRAN,MXVEC, 681 & 0,IOPSOP,IKAPA,IKAPC,0,0,ITAMPB,ITRAN,IVEC) 682 E2CON2 = WORK(KECONS-1 + (ITRAN-1)*MXVEC + IVEC) 683 ELSE 684 X2CON2 = 0.0d0 685 E2CON2 = 0.0d0 686 END IF 687 688 IF (LPDBSC.OR.LORXC .OR. LPDBSB.OR.LORXB) THEN 689 CALL CC_FIND_SO_OP(LBLOPR(IOPERC),LBLOPR(IOPERB), 690 & LABSOP,ISYSOP,ISGNSOP,INUM,WORK,LWORK) 691 IOPSOP = IROPER(LABSOP,ISYSOP) 692 693 CALL CC_SETXE('Xi ',WORK(KXETRAN),WORK(KXDOTS),MXTRAN,MXVEC, 694 & 0,IOPSOP,IKAPC,IKAPB,0,0,IZETVA,ITRAN,IVEC) 695 X2CON3 = WORK(KXCONS-1 + (ITRAN-1)*MXVEC + IVEC) 696 697 CALL CC_SETXE('Eta',WORK(KXETRAN),WORK(KEDOTS),MXTRAN,MXVEC, 698 & 0,IOPSOP,IKAPC,IKAPB,0,0,ITAMPA,ITRAN,IVEC) 699 E2CON3 = WORK(KECONS-1 + (ITRAN-1)*MXVEC + IVEC) 700 ELSE 701 X2CON3 = 0.0d0 702 E2CON3 = 0.0d0 703 END IF 704*---------------------------------------------------------------------* 705 706 IDX = NQRFREQ*(IOPER-1)+IFREQ 707 708 IF (ISIGN.EQ.-1) IDX = IDX + NBHYPPOL 709 710 WORK(KHYPPOL-1+IDX) = WORK(KHYPPOL-1+IDX) + 711 & GCON + 712 & FACON1 + FACON2 + FACON3 + 713 & FCON1 + FCON2 + FCON3 + 714 & BCON1 + BCON2 + BCON3 + 715 & SLVKCON1 + SLVKCON2 + SLVKCON3 + 716 & R2CON1 + R2CON2 + R2CON3 + 717 & EACON1 + EACON2 + EACON3 + EACON4 + EACON5 + EACON6 + 718 & E2CON1 + E2CON2 + E2CON3 + X2CON1 + X2CON2 + X2CON3 719 720*---------------------------------------------------------------------* 721 722 IF (LOCDBG .OR. IPRQHYP.GE.15) THEN 723 WRITE(LUPRI,*) 'ITAMPA,ITAMPB,ITAMPC:',ITAMPA,ITAMPB,ITAMPC 724 WRITE(LUPRI,*) 'IZETVA,IZETVB,IZETVC:',IZETVA,IZETVB,IZETVC 725 WRITE(LUPRI,*) 'IOPERA,IOPERB,IOPERC:',IOPERA,IOPERB,IOPERC 726 WRITE(LUPRI,*) 'IKAPA, IKAPB, IKAPC :',IKAPA, IKAPB, IKAPC 727 WRITE(LUPRI,*) 'GCON :',GCON 728 WRITE(LUPRI,*) 'FACON:',FACON1,FACON2,FACON3 729 WRITE(LUPRI,*) 'FCON :',FCON1,FCON2,FCON3 730 WRITE(LUPRI,*) 'BCON :',BCON1,BCON2,BCON3 731 WRITE(LUPRI,*) 'SLVKCON :',SLVKCON1,SLVKCON2,SLVKCON3 732 WRITE(LUPRI,*) 'EACON:',EACON1,EACON2,EACON3, 733 & EACON4,EACON5,EACON6 734 WRITE(LUPRI,*) 'E2CON:',E2CON1,E2CON2,E2CON3 735 WRITE(LUPRI,*) 'X2CON:',X2CON1,X2CON2,X2CON3 736 WRITE(LUPRI,*) 'R2CON:',R2CON1,R2CON2,R2CON3 737 WRITE(LUPRI,*) 'HYPPOL:',WORK(KHYPPOL-1+IDX) 738 END IF 739 740 END DO 741 END DO 742 END IF 743 END DO 744 745*---------------------------------------------------------------------* 746* print timing: 747*---------------------------------------------------------------------* 748 WRITE (LUPRI,'(/A,I4,A,F12.2," seconds.")') ' Total time for', 749 & NBHYPPOL,' quadratic response func.:', SECOND() - TIM0 750 751*---------------------------------------------------------------------* 752* print frequency-dependent hyperpolarizabilities: 753*---------------------------------------------------------------------* 754 755 CALL CCQHYPPRT(WORK(KHYPPOL)) 756 757 CALL FLSHFO(LUPRI) 758 759 IF (NQRDISP.EQ.0) THEN 760 RETURN 761 END IF 762 763*---------------------------------------------------------------------* 764* allocate & initialize work space for expansion coefficients 765*---------------------------------------------------------------------* 766 NBEXPCOF = NQROPER * NQRDISP 767 768 MXTRAN = 2 * NLRCLBL * NLRCLBL 769 MXVEC = NLRCLBL 770 771 MXGTRAN = 4 * MXTRAN 772 MXFTRAN = 3 * MXTRAN 773 MXBTRAN = 3 * MXTRAN 774 MXFATRAN = 5 * MXTRAN 775 MXEATRAN = 3 * MXTRAN 776 MXR2TRAN = 1 * MXTRAN 777 778 MXGDOTS = MXVEC * MXTRAN 779 MXFDOTS = MXVEC * MXTRAN 780 MXBDOTS = MXVEC * MXTRAN 781 MXFADOTS = MXVEC * MXTRAN 782 MXEADOTS = MXVEC * MXTRAN 783 MXR2DOTS = MXVEC * MXTRAN 784 785 KEXPCOF = 1 786 KGTRAN = KEXPCOF + NBEXPCOF 787 KGDOTS = KGTRAN + MXGTRAN 788 KGCONS = KGDOTS + MXGDOTS 789 KFTRAN = KGCONS + MXGDOTS 790 KFDOTS = KFTRAN + MXFTRAN 791 KFCONS = KFDOTS + MXFDOTS 792 KBTRAN = KFCONS + MXFDOTS 793 KBDOTS = KBTRAN + MXBTRAN 794 KBCONS = KBDOTS + MXBDOTS 795 KFATRAN = KBCONS + MXBDOTS 796 KFADOTS = KFATRAN + MXFATRAN 797 KFACONS = KFADOTS + MXFADOTS 798 KEATRAN = KFACONS + MXFADOTS 799 KEADOTS = KEATRAN + MXEATRAN 800 KEACONS = KEADOTS + MXEADOTS 801 KR2TRAN = KEACONS + MXEADOTS 802 KR2DOTS = KR2TRAN + MXR2TRAN 803 KR2CONS = KR2DOTS + MXR2DOTS 804 KEND0 = KR2CONS + MXR2DOTS 805 LEND0 = LWORK - KEND0 806 807 IF (LEND0 .LT. 0) THEN 808 CALL QUIT('Insufficient memory in CC_HYPPOL. (2)') 809 END IF 810 811* initialze hyperpolarizabilities and the lists of contributions: 812 CALL DZERO(WORK(KEXPCOF),NBEXPCOF) 813 CALL DZERO(WORK(KGTRAN), MXGTRAN +2*MXGDOTS) 814 CALL DZERO(WORK(KFTRAN), MXFTRAN +2*MXFDOTS) 815 CALL DZERO(WORK(KBTRAN), MXBTRAN +2*MXBDOTS) 816 CALL DZERO(WORK(KFATRAN),MXFATRAN+2*MXFADOTS) 817 CALL DZERO(WORK(KEATRAN),MXEATRAN+2*MXEADOTS) 818 CALL DZERO(WORK(KR2TRAN),MXR2TRAN+2*MXR2DOTS) 819 820* initialize timing: 821 TIM2 = SECOND() 822 823*---------------------------------------------------------------------* 824* set up lists for G, F/B and F{A} transformations and ETA{O} vectors: 825*---------------------------------------------------------------------* 826 LADD = .FALSE. 827 828 CALL CCQR_DISP_SETUP(MXTRAN, MXVEC, 829 & WORK(KGTRAN), WORK(KGDOTS), WORK(KGCONS), NGTRAN, 830 & WORK(KFTRAN), WORK(KFDOTS), WORK(KFCONS), NFTRAN, 831 & WORK(KBTRAN), WORK(KBDOTS), WORK(KBCONS), NBTRAN, 832 & WORK(KFATRAN),WORK(KFADOTS),WORK(KFACONS),NFATRAN, 833 & WORK(KEATRAN),WORK(KEADOTS),WORK(KEACONS),NEATRAN, 834 & WORK(KR2TRAN),WORK(KR2DOTS),WORK(KR2CONS),NR2TRAN, 835 & WORK(KEXPCOF),NBEXPCOF, LADD ) 836 837*---------------------------------------------------------------------* 838* calculate G matrix contributions: 839*---------------------------------------------------------------------* 840 TIM1 = SECOND() 841 842 IOPT = 5 843 CALL CC_GMATRIX('L0 ','RC ','RC ','RC ',NGTRAN, MXVEC, 844 & WORK(KGTRAN),WORK(KGDOTS),WORK(KGCONS), 845 & WORK(KEND0), LEND0, IOPT ) 846 847 TIMG = TIMG + SECOND() - TIM1 848 849 WRITE (LUPRI,'(/A,I5,A,F12.2," seconds.")') ' Time used for', 850 & NGTRAN,' G matrix transformations:', SECOND() - TIM1 851 CALL FLSHFO(LUPRI) 852 853*---------------------------------------------------------------------* 854* calculate B matrix contributions: 855*---------------------------------------------------------------------* 856 IF ( USEBTRAN .AND. (.NOT.USE_R2) ) THEN 857 858 TIM1 = SECOND() 859 860 IOPT = 5 861 CALL CC_BMATRIX(WORK(KBTRAN),NBTRAN,'RC ','RC ',IOPT,'LC ', 862 & WORK(KBDOTS),WORK(KBCONS),MXVEC, 863 & .FALSE.,WORK(KEND0),LEND0) 864 865 TIMB = TIMB + SECOND() - TIM1 866 867 WRITE (LUPRI,'(/A,I5,A,F12.2," seconds.")') ' Time used for', 868 & NBTRAN,' B matrix transformations:', SECOND() - TIM1 869 CALL FLSHFO(LUPRI) 870 871 END IF 872 873*---------------------------------------------------------------------* 874* calculate F matrix contributions: 875*---------------------------------------------------------------------* 876 IF ( (.NOT.USEBTRAN) .AND. (.NOT.USE_R2) ) THEN 877 878 TIM1 = SECOND() 879 880 IOPT = 5 881 CALL CC_FMATRIX(WORK(KFTRAN),NFTRAN,'LC ','RC ',IOPT,'RC ', 882 & WORK(KFDOTS),WORK(KFCONS),MXVEC, 883 & WORK(KEND0), LEND0) 884 885 TIMF = TIMF + SECOND() - TIM1 886 887 WRITE (LUPRI,'(/A,I5,A,F12.2," seconds.")') ' Time used for', 888 & NFTRAN,' F matrix transformations:', SECOND() - TIM1 889 CALL FLSHFO(LUPRI) 890 891 END IF 892 893*---------------------------------------------------------------------* 894* calculate F{O} matrix contributions: 895*---------------------------------------------------------------------* 896 TIM1 = SECOND() 897 898 CALL CCQR_FADRV('L0 ','o1 ','RC ','RC ',NFATRAN, MXVEC, 899 & WORK(KFATRAN),WORK(KFADOTS),WORK(KFACONS), 900 & WORK(KEND0), LEND0, 'DOTP' ) 901 902 TIMFA = TIMFA + SECOND() - TIM1 903 904 WRITE (LUPRI,'(/A,I5,A,F12.2," seconds.")') ' Time used for', 905 & NFATRAN,' F{O} matrix transform.: ', SECOND() - TIM1 906 CALL FLSHFO(LUPRI) 907 908*---------------------------------------------------------------------* 909* calculate ETA{O} vector contributions: 910*---------------------------------------------------------------------* 911 IF ( .NOT. USE_R2 ) THEN 912 TIM1 = SECOND() 913 914 CALL CCQR_EADRV('LC ','o1 ','RC ',NEATRAN, MXVEC, 915 & WORK(KEATRAN),WORK(KEADOTS),WORK(KEACONS), 916 & WORK(KEND0), LEND0, 'DOTP' ) 917 918 TIMEA = TIMEA + SECOND() - TIM1 919 920 WRITE (LUPRI,'(/A,I5,A,F12.2," seconds.")') ' Time used for', 921 & NEATRAN,' ETA{O} vector calculat.: ', SECOND() - TIM1 922 CALL FLSHFO(LUPRI) 923 924 ELSE 925 TIMEA = 0.0d0 926 END IF 927 928*---------------------------------------------------------------------* 929* calculate X1 x R2 dots products: 930*---------------------------------------------------------------------* 931 IF ( USE_R2 ) THEN 932 TIM1 = SECOND() 933 934 CALL CC_DOTDRV('CR2','XC ',NR2TRAN,MXVEC, 935 & WORK(KR2TRAN), WORK(KR2DOTS), WORK(KR2CONS), 936 & WORK(KEND0), LEND0 ) 937 938 TIMR2 = SECOND() - TIM1 939 WRITE (LUPRI,'(/A,I5,A,F12.2," seconds.")') ' Time used for', 940 & NR2TRAN,' XC x CR2 dot products: ', SECOND() - TIM1 941 CALL FLSHFO(LUPRI) 942 943 ELSE 944 TIMR2 = 0.0d0 945 END IF 946 947*---------------------------------------------------------------------* 948* collect contributions and add them to hyperpolarizabilities 949*---------------------------------------------------------------------* 950 LADD = .TRUE. 951 952 CALL CCQR_DISP_SETUP(MXTRAN, MXVEC, 953 & WORK(KGTRAN), WORK(KGDOTS), WORK(KGCONS), NGTRAN, 954 & WORK(KFTRAN), WORK(KFDOTS), WORK(KFCONS), NFTRAN, 955 & WORK(KBTRAN), WORK(KBDOTS), WORK(KBCONS), NBTRAN, 956 & WORK(KFATRAN),WORK(KFADOTS),WORK(KFACONS),NFATRAN, 957 & WORK(KEATRAN),WORK(KEADOTS),WORK(KEACONS),NEATRAN, 958 & WORK(KR2TRAN),WORK(KR2DOTS),WORK(KR2CONS),NR2TRAN, 959 & WORK(KEXPCOF),NBEXPCOF, LADD ) 960 961*---------------------------------------------------------------------* 962* print timing: 963*---------------------------------------------------------------------* 964 WRITE (LUPRI,'(/A,I4,A,F12.2," seconds.")') ' Total time for', 965 & NBEXPCOF,' expansion coeffients: ', SECOND() - TIM2 966 967*---------------------------------------------------------------------* 968* print expansion coefficients for first hyperpolarizabilities: 969*---------------------------------------------------------------------* 970 971 IF (IPRQHYP .GE. 15) THEN 972 CALL CCQEXPPRT(WORK(KEXPCOF)) 973 END IF 974 975 CALL FLSHFO(LUPRI) 976 977*---------------------------------------------------------------------* 978* calculate from the d_ABC expansion coefficients the D_ABC disp. cof. 979*---------------------------------------------------------------------* 980 NBINOM = (NQRDSPE+2)*(NQRDSPE+1)/2 981 982 KBINOM = KEND0 983 KDSPCF = KBINOM + NBINOM 984 KSHGCF = KDSPCF + NBINOM*NQROPER 985 KORCOF = KSHGCF + (NQRDSPE+1)*NQROPER 986 KEOPEC = KORCOF + (NQRDSPE+1)*NQROPER 987 KEND1 = KEOPEC + (NQRDSPE+1)*NQROPER 988 LEND1 = LWORK - KEND1 989 990 IF (BETA_AVERAGE) THEN 991 KDSPCFA = KEND1 992 KSHGCFA = KDSPCFA + 4*NBINOM 993 KORCOFA = KSHGCFA + 4*(NQRDSPE+1) 994 KEOPECA = KORCOFA + 4*(NQRDSPE+1) 995 KABCDE = KEOPECA + 4*(NQRDSPE+1) 996 KEND1 = KABCDE + 12*2 997 LEND1 = LWORK - KEND1 998 END IF 999 1000 IF (LEND1 .LT. 0) THEN 1001 CALL QUIT('Insufficient memory in CC_HYPPOL. (3)') 1002 END IF 1003 1004 CALL DZERO(WORK(KBINOM),KEND1-KBINOM) 1005 1006 CALL CCQRDISP(WORK(KDSPCF), WORK(KDSPCFA),WORK(KABCDE), 1007 & WORK(KSHGCF), WORK(KORCOF), WORK(KEOPEC), 1008 & WORK(KSHGCFA),WORK(KORCOFA),WORK(KEOPECA), 1009 & WORK(KEXPCOF),WORK(KBINOM)) 1010 1011*---------------------------------------------------------------------* 1012* print dispersion coefficients for first hyperpolarizabilities: 1013*---------------------------------------------------------------------* 1014 1015 CALL CCQDISPRT(WORK(KDSPCF), WORK(KDSPCFA),WORK(KABCDE), 1016 & WORK(KSHGCF), WORK(KORCOF), WORK(KEOPEC), 1017 & WORK(KSHGCFA),WORK(KORCOFA),WORK(KEOPECA)) 1018 1019*---------------------------------------------------------------------* 1020* print timing & return: 1021*---------------------------------------------------------------------* 1022 WRITE (LUPRI,'(/A,A,F12.2," seconds.")') ' Total time used ', 1023 & 'in quadratic response section:', SECOND() - TIM0 1024 1025 CALL FLSHFO(LUPRI) 1026 1027 RETURN 1028 END 1029 1030*=====================================================================* 1031* END OF SUBROUTINE CC_HYPPOL * 1032*=====================================================================* 1033c /* deck ccqrdisp */ 1034*=====================================================================* 1035 SUBROUTINE CCQRDISP(DISPCF,AVEDISP,ABCDE, 1036 & SHGCF, ORCOF,EOPECF, 1037 & AVESHG,AVEOR,AVEEOPE, 1038 & EXPCOF,BINOM) 1039*---------------------------------------------------------------------* 1040* 1041* Purpose: calculate from the expansion coefficients defined by 1042* beta_ABC = sum_{klm} w_A^k w_B^l w_C^m d_ABC(klm) the 1043* physical more relevant dispersion coefficients defined 1044* by beta_ABC = sum_{mn} w_B^m w_C^n D_ABC(mn) 1045* 1046* 1047* Written by Christof Haettig in October 1997. 1048* 1049*=====================================================================* 1050#if defined (IMPLICIT_NONE) 1051 IMPLICIT NONE 1052#else 1053# include "implicit.h" 1054#endif 1055#include "priunit.h" 1056#include "ccqrinf.h" 1057#include "ccroper.h" 1058 1059 LOGICAL LOCDBG 1060 PARAMETER (LOCDBG = .FALSE.) 1061 1062 INTEGER KA, KB, KC, KCP, KD, KDP, KE, KEP 1063 PARAMETER (KA=1,KB=2,KC=3,KCP=4,KD=5,KDP=6,KE=7,KEP=8) 1064 INTEGER KBP, KDP2, KEP2, KEP3 1065 PARAMETER (KBP=9,KDP2=10,KEP2=11,KEP3=12) 1066 1067 INTEGER IDISP, ICOMP, ISAMA, ISAMB, ISAMC, ISAPROP, MN0, MNS 1068 1069 1070#if defined (SYS_CRAY) 1071 REAL DISPCF((NQRDSPE+2)*(NQRDSPE+1)/2,NQROPER) 1072 REAL AVEDISP((NQRDSPE+2)*(NQRDSPE+1)/2,4) 1073 REAL ABCDE(12,2) 1074 REAL SHGCF(NQRDSPE+1,NQROPER) 1075 REAL ORCOF(NQRDSPE+1,NQROPER) 1076 REAL EOPECF(NQRDSPE+1,NQROPER) 1077 REAL AVESHG(NQRDSPE+1,4) 1078 REAL AVEOR(NQRDSPE+1,4) 1079 REAL AVEEOPE(NQRDSPE+1,4) 1080 REAL EXPCOF(NQRDISP,NQROPER) 1081 REAL BINOM((NQRDSPE+2)*(NQRDSPE+1)/2) 1082 REAL SIGNP, SIGNPQ, SIGNN, ERROR 1083 REAL BETA0,A0,A1,A2,B0,B1, B2, B3, B4, C0, C1, C2, C3 1084 REAL CP0,CP1,CP2,D0,D1,D2,D3, DP0, DP1, DP2, DP3, DP4 1085 REAL E0,E1,E2,E3, EP0, EP1, EP2, EP3, EP4, EP5, EP6 1086 REAL BP0,BP1,CP3,DPP0, DPP1, DPP2, DPP3 1087 REAL EPP0, EPP1, EPP2, EPPP0, EPPP1, EPPP2, EPPP3 1088 REAL BORTFAC(7), BM23FAC(7) 1089#else 1090 DOUBLE PRECISION DISPCF((NQRDSPE+2)*(NQRDSPE+1)/2,NQROPER) 1091 DOUBLE PRECISION AVEDISP((NQRDSPE+2)*(NQRDSPE+1)/2,4) 1092 DOUBLE PRECISION ABCDE(12,2) 1093 DOUBLE PRECISION SHGCF(NQRDSPE+1,NQROPER) 1094 DOUBLE PRECISION ORCOF(NQRDSPE+1,NQROPER) 1095 DOUBLE PRECISION EOPECF(NQRDSPE+1,NQROPER) 1096 DOUBLE PRECISION AVESHG(NQRDSPE+1,4) 1097 DOUBLE PRECISION AVEOR(NQRDSPE+1,4) 1098 DOUBLE PRECISION AVEEOPE(NQRDSPE+1,4) 1099 DOUBLE PRECISION EXPCOF(NQRDISP,NQROPER) 1100 DOUBLE PRECISION BINOM((NQRDSPE+2)*(NQRDSPE+1)/2) 1101 DOUBLE PRECISION SIGNP, SIGNPQ, SIGNN, ERROR 1102 DOUBLE PRECISION BETA0,A0,A1,A2,B0,B1, B2, B3, B4, C0, C1, C2, C3 1103 DOUBLE PRECISION CP0,CP1,CP2,D0,D1,D2,D3, DP0, DP1, DP2, DP3, DP4 1104 DOUBLE PRECISION E0,E1,E2,E3, EP0, EP1, EP2, EP3, EP4, EP5, EP6 1105 DOUBLE PRECISION BP0,BP1,CP3,DPP0, DPP1, DPP2, DPP3 1106 DOUBLE PRECISION EPP0, EPP1, EPP2, EPPP0, EPPP1, EPPP2, EPPP3 1107 DOUBLE PRECISION BORTFAC(7), BM23FAC(7), BKERFAC(7) 1108#endif 1109C 1110 DATA BORTFAC / 0.2d0, 0.4d0,-0.6d0,0.4d0, 0.4d0,-0.6d0,0.4d0 / 1111 DATA BM23FAC / 0.0d0, 0.5d0,-1.0d0,0.5d0, 0.5d0,-1.0d0,0.5d0 / 1112 DATA BKERFAC / 0.6d0,-0.3d0,1.2d0,-0.3d0,-0.3d0,1.2d0,-0.3d0 / 1113C 1114 INTEGER N, M, MN, P, Q, I, J, ITRI, IOPER, OFFEX, IEXPCF 1115 ITRI(I,J) = (MAX(I,J)+1)*(MAX(I,J)-2)/2 + I + J + 2 1116C 1117 1118 DO M = 0, NQRDSPE 1119 BINOM(ITRI(M,0)) = 1.0d0 1120 DO N = 1, M-1 1121 BINOM(ITRI(M,N)) = 1122 & BINOM(ITRI(M-1,N-1)) + BINOM(ITRI(M-1,N)) 1123 END DO 1124 BINOM(ITRI(M,M)) = 1.0d0 1125 END DO 1126 1127 IF (LOCDBG) THEN 1128 WRITE (LUPRI,*) 'binomial coefficients:' 1129 DO M = 0, NQRDSPE 1130 WRITE (LUPRI,'(10F8.2)') (BINOM(ITRI(M,N)),N=0,M) 1131 END DO 1132 END IF 1133 1134 DO IOPER = 1, NQROPER 1135 ISAMA = ISYMAT(IAQROP(IOPER)) 1136 ISAMB = ISYMAT(IBQROP(IOPER)) 1137 ISAMC = ISYMAT(ICQROP(IOPER)) 1138 ISAPROP = ISAMA * ISAMB * ISAMC 1139 IF (ALLDSPCF) ISAPROP = 0 1140 1141 MN0 = 0 1142 MNS = 2 1143 1144 IF (ISAPROP.EQ.-1) MN0 = 1 1145 IF (ISAPROP.EQ. 0) MNS = 1 1146 1147 IF (LOCDBG) THEN 1148 WRITE (LUPRI,*) 'LABELS:',LBLOPR(IAQROP(IOPER)), 1149 & LBLOPR(IBQROP(IOPER)), LBLOPR(ICQROP(IOPER)) 1150 WRITE (LUPRI,*) 'ISAPROP,MN0,MNS:',ISAPROP,MN0,MNS 1151 END IF 1152 1153 OFFEX = 0 1154 IF (MN0.EQ.1) OFFEX = 1 1155 DO MN = MN0, NQRDSPE, MNS 1156 IF (LOCDBG) THEN 1157 WRITE (LUPRI,'(//A,62X,A,10X,A,10X,A)') 1158 & ' M N ', ' SHGCF ','ORCOF','EOPECF,' 1159 END IF 1160 SHGCF(MN+1,IOPER) = 0.0d0 1161 ORCOF(MN+1,IOPER) = 0.0d0 1162 EOPECF(MN+1,IOPER) = 0.0d0 1163 SIGNN = 1.0d0 1164 DO N = 0, MN 1165 IF (LOCDBG) THEN 1166 WRITE (LUPRI,'(//A,A)') 1167 & ' M N P Q SIGN BINOM', 1168 & ' I J K EXPCOF DISPCF' 1169 END IF 1170 M = MN - N 1171 DISPCF(ITRI(M+N,N),IOPER) = 0.0d0 1172 SIGNP = +1.0d0 1173 DO P = 0, M 1174 SIGNPQ = SIGNP 1175 DO Q = 0, N 1176 IEXPCF = OFFEX + ITRI(M+N-P-Q,N-Q) 1177 DISPCF(ITRI(M+N,N),IOPER) = 1178 & DISPCF(ITRI(M+N,N),IOPER) + 1179 & SIGNPQ * BINOM(ITRI(P+Q,P)) * EXPCOF(IEXPCF,IOPER) 1180 IF (LOCDBG) THEN 1181 WRITE (LUPRI,'(4I3,2F8.2,3X,3I2,4E16.8)') 1182 & M,N,P,Q,SIGNPQ,BINOM(ITRI(P+Q,P)),P+Q,M-P,N-Q, 1183 & EXPCOF(IEXPCF,IOPER),DISPCF(ITRI(M+N,N),IOPER) 1184 END IF 1185 SIGNPQ = -SIGNPQ 1186 END DO 1187 SIGNP = -SIGNP 1188 END DO 1189 SHGCF(MN+1,IOPER) = SHGCF(MN+1,IOPER) + 1190 & DISPCF(ITRI(M+N,N),IOPER) 1191 ORCOF(MN+1,IOPER) = ORCOF(MN+1,IOPER) + 1192 & SIGNN * DISPCF(ITRI(M+N,N),IOPER) 1193 IF (N.EQ.0) EOPECF(MN+1,IOPER)=DISPCF(ITRI(M+N,N),IOPER) 1194 SIGNN = -SIGNN 1195 IF (LOCDBG) THEN 1196 WRITE (LUPRI,'(2I3,66X,3E16.8)') M,N,SHGCF(MN+1,IOPER), 1197 & ORCOF(MN+1,IOPER),EOPECF(MN+1,IOPER) 1198 END IF 1199 END DO 1200 OFFEX = OFFEX + (MN+2)*(MN+1)/2 1201 IF (MNS.EQ.2) OFFEX = OFFEX + (MN+3)*(MN+2)/2 1202 END DO 1203 END DO 1204 1205 IF (BETA_AVERAGE) THEN 1206 1207* calculate averaged vector components parallel/orthogonal to z axis: 1208 DO IDISP = 1, (NQRDSPE+2)*(NQRDSPE+1)/2 1209 AVEDISP(IDISP,1) = 0.6d0*DISPCF(IDISP,1) 1210 AVEDISP(IDISP,2) = BORTFAC(1)*DISPCF(IDISP,1) 1211 AVEDISP(IDISP,3) = BM23FAC(1)*DISPCF(IDISP,1) 1212 AVEDISP(IDISP,4) = BKERFAC(1)*DISPCF(IDISP,1) 1213 END DO 1214 DO IDISP = 1, NQRDSPE+1 1215 AVESHG(IDISP,1) = 0.6d0*SHGCF(IDISP,1) 1216 AVEOR(IDISP,1) = 0.6d0*ORCOF(IDISP,1) 1217 AVEEOPE(IDISP,1) = 0.6d0*EOPECF(IDISP,1) 1218 AVESHG(IDISP,2) = BORTFAC(1)*SHGCF(IDISP,1) 1219 AVEOR(IDISP,2) = BORTFAC(1)*ORCOF(IDISP,1) 1220 AVEEOPE(IDISP,2) = BORTFAC(1)*EOPECF(IDISP,1) 1221 AVESHG(IDISP,3) = BM23FAC(1)*SHGCF(IDISP,1) 1222 AVEOR(IDISP,3) = BM23FAC(1)*ORCOF(IDISP,1) 1223 AVEEOPE(IDISP,3) = BM23FAC(1)*EOPECF(IDISP,1) 1224 AVESHG(IDISP,4) = BKERFAC(1)*SHGCF(IDISP,1) 1225 AVEOR(IDISP,4) = BKERFAC(1)*ORCOF(IDISP,1) 1226 AVEEOPE(IDISP,4) = BKERFAC(1)*EOPECF(IDISP,1) 1227 END DO 1228 IF (XY_DEGENERAT) THEN 1229 DO ICOMP = 2, 4 1230 DO IDISP = 1, (NQRDSPE+2)*(NQRDSPE+1)/2 1231 AVEDISP(IDISP,1)= AVEDISP(IDISP,1)+0.4d0*DISPCF(IDISP,ICOMP) 1232 AVEDISP(IDISP,2)= AVEDISP(IDISP,2) + 1233 & 2.0d0*BORTFAC(ICOMP)*DISPCF(IDISP,ICOMP) 1234 AVEDISP(IDISP,3)= AVEDISP(IDISP,3) + 1235 & 2.0d0*BM23FAC(ICOMP)*DISPCF(IDISP,ICOMP) 1236 AVEDISP(IDISP,4)= AVEDISP(IDISP,4) + 1237 & 2.0d0*BKERFAC(ICOMP)*DISPCF(IDISP,ICOMP) 1238 END DO 1239 DO IDISP = 1, NQRDSPE+1 1240 AVESHG(IDISP,1) = AVESHG(IDISP,1) +0.4d0*SHGCF(IDISP,ICOMP) 1241 AVEOR(IDISP,1) = AVEOR(IDISP,1) +0.4d0*ORCOF(IDISP,ICOMP) 1242 AVEEOPE(IDISP,1)= AVEEOPE(IDISP,1)+0.4d0*EOPECF(IDISP,ICOMP) 1243 AVESHG(IDISP,2) = AVESHG(IDISP,2) + 1244 & 2.0d0*BORTFAC(ICOMP)*SHGCF(IDISP,ICOMP) 1245 AVEOR(IDISP,2) = AVEOR(IDISP,2) + 1246 & 2.0d0*BORTFAC(ICOMP)*ORCOF(IDISP,ICOMP) 1247 AVEEOPE(IDISP,2)= AVEEOPE(IDISP,2)+ 1248 & 2.0d0*BORTFAC(ICOMP)*EOPECF(IDISP,ICOMP) 1249 AVESHG(IDISP,3) = AVESHG(IDISP,3) + 1250 & 2.0d0*BM23FAC(ICOMP)*SHGCF(IDISP,ICOMP) 1251 AVEOR(IDISP,3) = AVEOR(IDISP,3) + 1252 & 2.0d0*BM23FAC(ICOMP)*ORCOF(IDISP,ICOMP) 1253 AVEEOPE(IDISP,3)= AVEEOPE(IDISP,3)+ 1254 & 2.0d0*BM23FAC(ICOMP)*EOPECF(IDISP,ICOMP) 1255 AVESHG(IDISP,4) = AVESHG(IDISP,4) + 1256 & 2.0d0*BKERFAC(ICOMP)*SHGCF(IDISP,ICOMP) 1257 AVEOR(IDISP,4) = AVEOR(IDISP,4) + 1258 & 2.0d0*BKERFAC(ICOMP)*ORCOF(IDISP,ICOMP) 1259 AVEEOPE(IDISP,4)= AVEEOPE(IDISP,4)+ 1260 & 2.0d0*BKERFAC(ICOMP)*EOPECF(IDISP,ICOMP) 1261 END DO 1262 END DO 1263 ELSE 1264 DO ICOMP = 2, 7 1265 DO IDISP = 1, (NQRDSPE+2)*(NQRDSPE+1)/2 1266 AVEDISP(IDISP,1)= AVEDISP(IDISP,1)+0.2d0*DISPCF(IDISP,ICOMP) 1267 AVEDISP(IDISP,2)= AVEDISP(IDISP,2) + 1268 & BORTFAC(ICOMP)*DISPCF(IDISP,ICOMP) 1269 AVEDISP(IDISP,3)= AVEDISP(IDISP,3) + 1270 & BM23FAC(ICOMP)*DISPCF(IDISP,ICOMP) 1271 AVEDISP(IDISP,4)= AVEDISP(IDISP,4) + 1272 & BKERFAC(ICOMP)*DISPCF(IDISP,ICOMP) 1273 END DO 1274 DO IDISP = 1, NQRDSPE+1 1275 AVESHG(IDISP,1) = AVESHG(IDISP,1) +0.2d0*SHGCF(IDISP,ICOMP) 1276 AVEOR(IDISP,1) = AVEOR(IDISP,1) +0.2d0*ORCOF(IDISP,ICOMP) 1277 AVEEOPE(IDISP,1)= AVEEOPE(IDISP,1)+0.2d0*EOPECF(IDISP,ICOMP) 1278 AVESHG(IDISP,2) = AVESHG(IDISP,2) + 1279 & BORTFAC(ICOMP)*SHGCF(IDISP,ICOMP) 1280 AVEOR(IDISP,2) = AVEOR(IDISP,2) + 1281 & BORTFAC(ICOMP)*ORCOF(IDISP,ICOMP) 1282 AVEEOPE(IDISP,2)= AVEEOPE(IDISP,2)+ 1283 & BORTFAC(ICOMP)*EOPECF(IDISP,ICOMP) 1284 AVESHG(IDISP,3) = AVESHG(IDISP,3) + 1285 & BM23FAC(ICOMP)*SHGCF(IDISP,ICOMP) 1286 AVEOR(IDISP,3) = AVEOR(IDISP,3) + 1287 & BM23FAC(ICOMP)*ORCOF(IDISP,ICOMP) 1288 AVEEOPE(IDISP,3)= AVEEOPE(IDISP,3)+ 1289 & BM23FAC(ICOMP)*EOPECF(IDISP,ICOMP) 1290 AVESHG(IDISP,4) = AVESHG(IDISP,4) + 1291 & BKERFAC(ICOMP)*SHGCF(IDISP,ICOMP) 1292 AVEOR(IDISP,4) = AVEOR(IDISP,4) + 1293 & BKERFAC(ICOMP)*ORCOF(IDISP,ICOMP) 1294 AVEEOPE(IDISP,4)= AVEEOPE(IDISP,4)+ 1295 & BKERFAC(ICOMP)*EOPECF(IDISP,ICOMP) 1296 END DO 1297 END DO 1298 END IF 1299 1300 1301 IF (LOCDBG) THEN 1302 WRITE (LUPRI,'(//,A)') 1303 & 'DISPERSION COEFF. FOR PARALLEL/ORTHOGONAL/23/Kerr AVERAGE:' 1304 DO IDISP = 1, (NQRDSPE+2)*(NQRDSPE+1)/2 1305 END DO 1306 DO MN = 0, NQRDSPE 1307 DO M = 0, MN 1308 N = MN - M 1309 WRITE (LUPRI,'(2I5,4G16.8)') N, M, 1310 & (AVEDISP(ITRI(MN,M),I),I=1,4) 1311 END DO 1312 END DO 1313 1314 WRITE (LUPRI,'(//,A)') 1315 & 'SHG COEFFICIENTS FOR PARALLEL/ORTOGONAL/23/Kerr AVERAGE:' 1316 DO IDISP = 1, NQRDSPE+1 1317 WRITE (LUPRI,'(I5,4G16.8)') IDISP-1, (AVESHG(IDISP,I),I=1,4) 1318 END DO 1319 1320 WRITE (LUPRI,'(//,A)') 1321 & 'OR COEFFICIENTS FOR PARALLEL/ORTHOGONAL/23/Kerr AVERAGE:' 1322 DO IDISP = 1, NQRDSPE+1 1323 WRITE (LUPRI,'(I5,4G16.8)') IDISP-1, (AVEOR(IDISP,I),I=1,4) 1324 END DO 1325 1326 WRITE (LUPRI,'(//,A)') 1327 & 'EOPE COEFFICIENTS FOR PARALLEL/ORTHOGONAL/23/Kerr AVERAGE:' 1328 DO IDISP = 1, NQRDSPE+1 1329 WRITE (LUPRI,'(I5,4G16.8)') IDISP-1,(AVEEOPE(IDISP,I),I=1,4) 1330 END DO 1331 END IF 1332 1333* calculate A, B, C, etc. coefficients: 1334 BETA0 = AVEDISP(ITRI(0,0),1) 1335 IF (NQRDSPE.GE.2) THEN 1336 A0 = AVEDISP(ITRI(2,0),1) / (2.0d0 * BETA0) 1337 A1 = AVEDISP(ITRI(2,1),1) / (2.0d0 * BETA0) 1338 A2 = AVEDISP(ITRI(2,2),1) / (2.0d0 * BETA0) 1339 ERROR = ( DABS(A1-A0) + DABS(A2-A0) ) / DABS(A0) 1340 IF ( ERROR .GT. 1.0d-10 ) THEN 1341 WRITE (LUPRI,*) 'ERROR DURING CALCULATION OF A COEFFICIENT:' 1342 WRITE (LUPRI,*) 'A0:',A0 1343 WRITE (LUPRI,*) 'A1:',A1 1344 WRITE (LUPRI,*) 'A2:',A2 1345 CALL QUIT('ERROR DURING CALCULATION OF A COEFFICIENT.') 1346 END IF 1347 ABCDE(KA,1) = A0 1348 IF (LOCDBG) THEN 1349 WRITE (LUPRI,'(//,A)') 1350 & 'A,B,C,D,E COEFF. FOR PARALLEL/23 AVERAGES:' 1351 WRITE (LUPRI,*) ' A coefficient: ',A0 1352 END IF 1353 1354 A0 = AVEDISP(ITRI(2,0),3) / ( 1.0d0 * BETA0) 1355 A1 = AVEDISP(ITRI(2,1),3) / (-2.0d0 * BETA0) 1356 A2 = AVEDISP(ITRI(2,2),3) / (-2.0d0 * BETA0) 1357 ERROR = ( DABS(A1-A0) + DABS(A2-A0) ) / DABS(A0) 1358 IF ( ERROR .GT. 1.0d-10 ) THEN 1359 WRITE (LUPRI,*) 1360 & 'ERROR DURING CALCULATION OF A_23 COEFFICIENT:' 1361 WRITE (LUPRI,*) 'A0:',A0 1362 WRITE (LUPRI,*) 'A1:',A1 1363 WRITE (LUPRI,*) 'A2:',A2 1364 CALL QUIT('ERROR DURING CALCULATION OF A_23 COEFFICIENT.') 1365 END IF 1366 ABCDE(KA,2) = A0 1367 IF (LOCDBG) THEN 1368 WRITE (LUPRI,*) ' A_23 coefficient: ',A0 1369 END IF 1370 END IF 1371 IF (NQRDSPE.GE.4) THEN 1372 B0 = AVEDISP(ITRI(4,0),1) / ( 4.0d0 * BETA0) 1373 B1 = AVEDISP(ITRI(4,1),1) / ( 8.0d0 * BETA0) 1374 B2 = AVEDISP(ITRI(4,2),1) / (12.0d0 * BETA0) 1375 B3 = AVEDISP(ITRI(4,3),1) / ( 8.0d0 * BETA0) 1376 B4 = AVEDISP(ITRI(4,4),1) / ( 4.0d0 * BETA0) 1377 ERROR = ( DABS(B1-B0)+DABS(B2-B0)+DABS(B3-B0)+DABS(B4-B0) )/ 1378 & DABS(B0) 1379 IF ( ERROR .GT. 1.0d-10 ) THEN 1380 WRITE (LUPRI,*) 'ERROR DURING CALCULATION OF B COEFFICIENT:' 1381 WRITE (LUPRI,*) 'B0:',B0 1382 WRITE (LUPRI,*) 'B1:',B1 1383 WRITE (LUPRI,*) 'B2:',B2 1384 WRITE (LUPRI,*) 'B3:',B3 1385 WRITE (LUPRI,*) 'B4:',B4 1386 CALL QUIT('ERROR DURING CALCULATION OF B COEFFICIENT.') 1387 END IF 1388 ABCDE(KB,1) = B0 1389 IF (LOCDBG) THEN 1390 WRITE (LUPRI,*) ' B coefficient: ',B0 1391 END IF 1392 1393 B0 = AVEDISP(ITRI(4,0),3) / ( 2.0d0 * BETA0) 1394 B1 = AVEDISP(ITRI(4,3),3) / (-8.0d0 * BETA0) 1395 B2 = AVEDISP(ITRI(4,4),3) / (-4.0d0 * BETA0) 1396 BP0 = AVEDISP(ITRI(4,1),3) + 1.0d0 * AVEDISP(ITRI(4,0),3) 1397 BP1 = AVEDISP(ITRI(4,2),3) + 3.0d0 * AVEDISP(ITRI(4,0),3) 1398 BP0 = BP0 / (-9.0d0 * BETA0) 1399 BP1 = BP1 / (-9.0d0 * BETA0) 1400 ERROR = ( DABS(B1-B0)+DABS(B2-B0)+DABS(BP1-BP0))/ DABS(B0) 1401 IF ( ERROR .GT. 1.0d-10 ) THEN 1402 WRITE (LUPRI,*) 1403 & 'ERROR DURING CALCULATION OF B_23 COEFFICIENT:' 1404 WRITE (LUPRI,*) 'B0:',B0 1405 WRITE (LUPRI,*) 'B1:',B1 1406 WRITE (LUPRI,*) 'B2:',B2 1407 WRITE (LUPRI,*) 'BP0:',BP0 1408 WRITE (LUPRI,*) 'BP1:',BP1 1409 CALL QUIT('ERROR DURING CALCULATION OF B_23 COEFFICIENT.') 1410 END IF 1411 ABCDE(KB,2) = B0 1412 ABCDE(KBP,2) = BP0 1413 IF (LOCDBG) THEN 1414 WRITE (LUPRI,*) " B_23, B_23' coefficients: ",B0, BP0 1415 END IF 1416 END IF 1417 IF (NQRDSPE.GE.6) THEN 1418 C0 = AVEDISP(ITRI(6,0),1) / ( 8.0d0 * BETA0) 1419 C1 = AVEDISP(ITRI(6,1),1) / (24.0d0 * BETA0) 1420 C2 = AVEDISP(ITRI(6,5),1) / (24.0d0 * BETA0) 1421 C3 = AVEDISP(ITRI(6,6),1) / ( 8.0d0 * BETA0) 1422 CP0 = AVEDISP(ITRI(6,2),1) - 6.0d0 * AVEDISP(ITRI(6,0),1) 1423 CP1 = AVEDISP(ITRI(6,3),1) - 7.0d0 * AVEDISP(ITRI(6,0),1) 1424 CP2 = AVEDISP(ITRI(6,4),1) - 6.0d0 * AVEDISP(ITRI(6,0),1) 1425 CP0 = CP0 / ( 9.0d0 * BETA0) 1426 CP1 = CP1 / (18.0d0 * BETA0) 1427 CP2 = CP2 / ( 9.0d0 * BETA0) 1428 ERROR = ( DABS(C1-C0) + DABS(C2-C0) + DABS(C3-C0) + 1429 & DABS(CP1-CP0) + DABS(CP2-CP0) ) / DABS(C0) 1430 IF ( ERROR .GT. 1.0d-10 ) THEN 1431 WRITE (LUPRI,*) 1432 & 'ERROR DURING CALCULATION OF C COEFFICIENTS:' 1433 WRITE (LUPRI,*) 'C0:',C0 1434 WRITE (LUPRI,*) 'C1:',C1 1435 WRITE (LUPRI,*) 'C2:',C2 1436 WRITE (LUPRI,*) 'C3:',C3 1437 WRITE (LUPRI,*) 'CP0:',CP0 1438 WRITE (LUPRI,*) 'CP1:',CP1 1439 WRITE (LUPRI,*) 'CP2:',CP2 1440 CALL QUIT('ERROR DURING CALCULATION OF C COEFFICIENTS.') 1441 END IF 1442 ABCDE(KC,1) = C0 1443 ABCDE(KCP,1) = CP0 1444 IF (LOCDBG) THEN 1445 WRITE (LUPRI,*) " C, C' coefficients:",C0,CP0 1446 END IF 1447 1448 C0 = AVEDISP(ITRI(6,0),3) / ( 4.0d0 * BETA0) 1449 C1 = AVEDISP(ITRI(6,5),3) / (-24.0d0 * BETA0) 1450 C2 = AVEDISP(ITRI(6,6),3) / ( -8.0d0 * BETA0) 1451 CP0 = AVEDISP(ITRI(6,1),3) 1452 CP1 = AVEDISP(ITRI(6,2),3) + 3.0d0 * AVEDISP(ITRI(6,0),3) 1453 CP2 = AVEDISP(ITRI(6,3),3) + 8.0d0 * AVEDISP(ITRI(6,0),3) 1454 CP3 = AVEDISP(ITRI(6,4),3) + 9.0d0 * AVEDISP(ITRI(6,0),3) 1455 CP0 = CP0 / (-18.0d0 * BETA0) 1456 CP1 = CP1 / (-36.0d0 * BETA0) 1457 CP2 = CP2 / (-36.0d0 * BETA0) 1458 CP3 = CP3 / (-18.0d0 * BETA0) 1459 ERROR = ( DABS(C1-C0) + DABS(C2-C0) + DABS(CP3-CP0) + 1460 & DABS(CP1-CP0) + DABS(CP2-CP0) ) / DABS(C0) 1461 IF ( ERROR .GT. 1.0d-10 ) THEN 1462 WRITE (LUPRI,*) 1463 & 'ERROR DURING CALCULATION OF C_23 COEFFICIENTS:' 1464 WRITE (LUPRI,*) 'C0:',C0 1465 WRITE (LUPRI,*) 'C1:',C1 1466 WRITE (LUPRI,*) 'C2:',C2 1467 WRITE (LUPRI,*) 'CP0:',CP0 1468 WRITE (LUPRI,*) 'CP1:',CP1 1469 WRITE (LUPRI,*) 'CP2:',CP2 1470 WRITE (LUPRI,*) 'CP3:',CP3 1471 CALL QUIT('ERROR DURING CALCULATION OF C_23 COEFFICIENTS.') 1472 END IF 1473 ABCDE(KC,2) = C0 1474 ABCDE(KCP,2) = CP0 1475 IF (LOCDBG) THEN 1476 WRITE (LUPRI,*) " C_23, C_23' coefficients:",C0,CP0 1477 END IF 1478 END IF 1479 IF (NQRDSPE.GE.8) THEN 1480 D0 = AVEDISP(ITRI(8,0),1) / (16.0d0 * BETA0) 1481 D1 = AVEDISP(ITRI(8,1),1) / (64.0d0 * BETA0) 1482 D2 = AVEDISP(ITRI(8,7),1) / (64.0d0 * BETA0) 1483 D3 = AVEDISP(ITRI(8,8),1) / (16.0d0 * BETA0) 1484 DP0 = AVEDISP(ITRI(8,2),1) - 10.0d0 * AVEDISP(ITRI(8,0),1) 1485 DP1 = AVEDISP(ITRI(8,3),1) - 16.0d0 * AVEDISP(ITRI(8,0),1) 1486 DP2 = AVEDISP(ITRI(8,4),1) - 19.0d0 * AVEDISP(ITRI(8,0),1) 1487 DP3 = AVEDISP(ITRI(8,5),1) - 16.0d0 * AVEDISP(ITRI(8,0),1) 1488 DP4 = AVEDISP(ITRI(8,6),1) - 10.0d0 * AVEDISP(ITRI(8,0),1) 1489 DP0 = DP0 / (18.0d0 * BETA0) 1490 DP1 = DP1 / (54.0d0 * BETA0) 1491 DP2 = DP2 / (72.0d0 * BETA0) 1492 DP3 = DP3 / (54.0d0 * BETA0) 1493 DP4 = DP4 / (18.0d0 * BETA0) 1494 ERROR = ( DABS(D1-D0) + DABS(D2-D0) + DABS(D3-D0) + 1495 & DABS(DP1-DP0) + DABS(DP2-DP0) + 1496 & DABS(DP3-DP0) + DABS(DP4-DP0) ) / DABS(D0) 1497 IF ( ERROR .GT. 1.0d-10 ) THEN 1498 WRITE (LUPRI,*) 1499 & 'ERROR DURING CALCULATION OF D COEFFICIENTS:' 1500 WRITE (LUPRI,*) 'D0:',D0 1501 WRITE (LUPRI,*) 'D1:',D1 1502 WRITE (LUPRI,*) 'D2:',D2 1503 WRITE (LUPRI,*) 'D3:',D3 1504 WRITE (LUPRI,*) 'DP0:',DP0 1505 WRITE (LUPRI,*) 'DP1:',DP1 1506 WRITE (LUPRI,*) 'DP2:',DP2 1507 WRITE (LUPRI,*) 'DP3:',DP3 1508 WRITE (LUPRI,*) 'DP4:',DP4 1509 CALL QUIT('ERROR DURING CALCULATION OF D COEFFICIENTS.') 1510 END IF 1511 ABCDE(KD,1) = D0 1512 ABCDE(KDP,1) = DP0 1513 IF (LOCDBG) THEN 1514 WRITE (LUPRI,*) " D, D' coefficients:",D0,DP0 1515 END IF 1516 1517 D0 = AVEDISP(ITRI(8,0),3) / ( 8.0d0 * BETA0) 1518 D1 = AVEDISP(ITRI(8,7),3) / (-64.0d0 * BETA0) 1519 D2 = AVEDISP(ITRI(8,8),3) / (-16.0d0 * BETA0) 1520 DP0 = AVEDISP(ITRI(8,1),3) - 1.0d0 * AVEDISP(ITRI(8,0),3) 1521 DP1 = AVEDISP(ITRI(8,3),3) + 11.0d0 * AVEDISP(ITRI(8,0),3) 1522 DPP0 = AVEDISP(ITRI(8,2),3) - 3.0d0 * AVEDISP(ITRI(8,1),3) + 1523 & 5.0d0 * AVEDISP(ITRI(8,0),3) 1524 DPP1 = AVEDISP(ITRI(8,4),3) - 5.0d0 * AVEDISP(ITRI(8,1),3) + 1525 & 25.0d0 * AVEDISP(ITRI(8,0),3) 1526 DPP2 = AVEDISP(ITRI(8,5),3) - 3.0d0 * AVEDISP(ITRI(8,1),3) + 1527 & 26.0d0 * AVEDISP(ITRI(8,0),3) 1528 DPP3 = AVEDISP(ITRI(8,6),3) - 1.0d0 * AVEDISP(ITRI(8,1),3) + 1529 & 18.0d0 * AVEDISP(ITRI(8,0),3) 1530 DP0 = DP0 / ( -36.0d0 * BETA0) 1531 DP1 = DP1 / ( -180.0d0 * BETA0) 1532 DPP0 = DPP0 / ( 9.0d0 * BETA0) 1533 DPP1 = DPP1 / ( -45.0d0 * BETA0) 1534 DPP2 = DPP2 / ( -54.0d0 * BETA0) 1535 DPP3 = DPP3 / ( -18.0d0 * BETA0) 1536 ERROR = ( DABS(D1-D0) + DABS(D2-D0) + 1537 & DABS(DP1-DP0) + DABS(DPP1-DPP0) + 1538 & DABS(DPP2-DPP0) + DABS(DPP3-DPP0) ) / DABS(D0) 1539 IF ( ERROR .GT. 1.0d-10 ) THEN 1540 WRITE (LUPRI,*) 1541 & 'ERROR DURING CALCULATION OF D COEFFICIENTS:' 1542 WRITE (LUPRI,*) 'D0:',D0 1543 WRITE (LUPRI,*) 'D1:',D1 1544 WRITE (LUPRI,*) 'D2:',D2 1545 WRITE (LUPRI,*) 'DP0:',DP0 1546 WRITE (LUPRI,*) 'DP1:',DP1 1547 WRITE (LUPRI,*) 'DPP0:',DPP0 1548 WRITE (LUPRI,*) 'DPP1:',DPP1 1549 WRITE (LUPRI,*) 'DPP2:',DPP2 1550 WRITE (LUPRI,*) 'DPP3:',DPP3 1551 CALL QUIT('ERROR DURING CALCULATION OF D COEFFICIENTS.') 1552 END IF 1553 ABCDE(KD,2) = D0 1554 ABCDE(KDP,2) = DP0 1555 ABCDE(KDP2,2) = DPP0 1556 IF (LOCDBG) THEN 1557 WRITE (LUPRI,'(A,3G16.8)') 1558 & " D_23, D_23', D_23'' coefficients:",D0,DP0,DPP0 1559 END IF 1560 END IF 1561 IF (NQRDSPE.GE.10) THEN 1562 E0 = AVEDISP(ITRI(10,0),1) / ( 32.0d0 * BETA0) 1563 E1 = AVEDISP(ITRI(10,1),1) / (160.0d0 * BETA0) 1564 E2 = AVEDISP(ITRI(10,9),1) / (160.0d0 * BETA0) 1565 E3 = AVEDISP(ITRI(10,10),1) / ( 32.0d0 * BETA0) 1566 EP0 = AVEDISP(ITRI(10,2),1) - 15.0d0 * AVEDISP(ITRI(10,0),1) 1567 EP1 = AVEDISP(ITRI(10,3),1) - 30.0d0 * AVEDISP(ITRI(10,0),1) 1568 EP2 = AVEDISP(ITRI(10,4),1) - 45.0d0 * AVEDISP(ITRI(10,0),1) 1569 EP3 = AVEDISP(ITRI(10,5),1) - 51.0d0 * AVEDISP(ITRI(10,0),1) 1570 EP4 = AVEDISP(ITRI(10,6),1) - 45.0d0 * AVEDISP(ITRI(10,0),1) 1571 EP5 = AVEDISP(ITRI(10,7),1) - 30.0d0 * AVEDISP(ITRI(10,0),1) 1572 EP6 = AVEDISP(ITRI(10,8),1) - 15.0d0 * AVEDISP(ITRI(10,0),1) 1573 EP0 = EP0 / ( 36.0d0 * BETA0) 1574 EP1 = EP1 / (144.0d0 * BETA0) 1575 EP2 = EP2 / (288.0d0 * BETA0) 1576 EP3 = EP3 / (360.0d0 * BETA0) 1577 EP4 = EP4 / (288.0d0 * BETA0) 1578 EP5 = EP5 / (144.0d0 * BETA0) 1579 EP6 = EP6 / ( 36.0d0 * BETA0) 1580 ERROR = ( DABS(E1-E0) + DABS(E2-E0) + DABS(E3-E0) + 1581 & DABS(EP1-EP0) + DABS(EP2-EP0) + 1582 & DABS(EP3-EP0) + DABS(EP4-EP0) + 1583 & DABS(EP5-EP0) + DABS(EP6-EP0) ) / DABS(E0) 1584 IF ( ERROR .GT. 1.0d-10 ) THEN 1585 WRITE (LUPRI,*) 1586 & 'ERROR DURING CALCULATION OF E COEFFICIENTS:' 1587 WRITE (LUPRI,*) 'E0: ',E0 1588 WRITE (LUPRI,*) 'E1: ',E1 1589 WRITE (LUPRI,*) 'E2: ',E2 1590 WRITE (LUPRI,*) 'E3: ',E3 1591 WRITE (LUPRI,*) 'EP0:',EP0 1592 WRITE (LUPRI,*) 'EP1:',EP1 1593 WRITE (LUPRI,*) 'EP2:',EP2 1594 WRITE (LUPRI,*) 'EP3:',EP3 1595 WRITE (LUPRI,*) 'EP4:',EP4 1596 WRITE (LUPRI,*) 'EP5:',EP5 1597 WRITE (LUPRI,*) 'EP6:',EP6 1598 CALL QUIT('ERROR DURING CALCULATION OF E COEFFICIENTS.') 1599 END IF 1600 ABCDE(KE,1) = E0 1601 ABCDE(KEP,1) = EP0 1602 IF (LOCDBG) THEN 1603 WRITE (LUPRI,*) " E, E' coefficients:",E0,EP0 1604 END IF 1605 1606 E0 = AVEDISP(ITRI(10,0),3) / ( 16.0d0 * BETA0) 1607 E1 = AVEDISP(ITRI(10,9),3) / (-160.0d0 * BETA0) 1608 E2 = AVEDISP(ITRI(10,10),3) / ( -32.0d0 * BETA0) 1609 EP0 = AVEDISP(ITRI(10,1),3) - 2.0d0 * AVEDISP(ITRI(10,0),3) 1610 EP1 = AVEDISP(ITRI(10,8),3) +27.0d0 * AVEDISP(ITRI(10,0),3)+ 1611 & 2.0d0 * AVEDISP(ITRI(10,2),3) 1612 EP0 = EP0 / (-72.0d0 * BETA0) 1613 EP1 = EP1 / (-648.0d0 * BETA0) 1614 1615 EPP0 = AVEDISP(ITRI(10,2),3) - 4.0d0*AVEDISP(ITRI(10,1),3) + 1616 & 8.0d0*AVEDISP(ITRI(10,0),3) 1617 EPP1 = AVEDISP(ITRI(10,7),3) - 4.0d0*AVEDISP(ITRI(10,1),3) + 1618 & 56.0d0*AVEDISP(ITRI(10,0),3) 1619 EPP2 = AVEDISP(ITRI(10,8),3) - 1.0d0*AVEDISP(ITRI(10,1),3) + 1620 & 29.0d0*AVEDISP(ITRI(10,0),3) 1621 EPP0 = EPP0 / ( 18.0d0 * BETA0) 1622 EPP1 = EPP1 / (-144.0d0 * BETA0) 1623 EPP2 = EPP2 / ( -36.0d0 * BETA0) 1624 1625 EPPP0= AVEDISP(ITRI(10,3),3) - 1.0d0 * AVEDISP(ITRI(10,2),3) 1626 & - 5.0d0*AVEDISP(ITRI(10,1),3) + 22.0d0*AVEDISP(ITRI(10,0),3) 1627 EPPP1= AVEDISP(ITRI(10,4),3) + 4.0d0 * AVEDISP(ITRI(10,2),3) 1628 & -29.0d0*AVEDISP(ITRI(10,1),3) + 91.0d0*AVEDISP(ITRI(10,0),3) 1629 EPPP2= AVEDISP(ITRI(10,5),3) + 11.0d0 * AVEDISP(ITRI(10,2),3) 1630 & -57.0d0*AVEDISP(ITRI(10,1),3) +168.0d0*AVEDISP(ITRI(10,0),3) 1631 EPPP3= AVEDISP(ITRI(10,6),3) + 13.0d0 * AVEDISP(ITRI(10,2),3) 1632 & -61.0d0*AVEDISP(ITRI(10,1),3) +182.0d0*AVEDISP(ITRI(10,0),3) 1633 EPPP0 = EPPP0 / ( -81.0d0 * BETA0) 1634 EPPP1 = EPPP1 / (-243.0d0 * BETA0) 1635 EPPP2 = EPPP2 / (-243.0d0 * BETA0) 1636 EPPP3 = EPPP3 / ( -81.0d0 * BETA0) 1637 ERROR = ( DABS(E1-E0) + DABS(E2-E0) + DABS(EP1-EP0) + 1638 & DABS(EPP1-EPP0) + DABS(EPP2-EPP0) + 1639 & DABS(EPPP1-EPPP0) + DABS(EPPP2-EPPP0) + 1640 & DABS(EPPP3-EPPP0) ) / DABS(E0) 1641 IF ( ERROR .GT. 1.0d-10 ) THEN 1642 WRITE (LUPRI,*) 1643 & 'ERROR DURING CALCULATION OF E_23 COEFFICIENTS:' 1644 WRITE (LUPRI,*) 'E0: ',E0 1645 WRITE (LUPRI,*) 'E1: ',E1 1646 WRITE (LUPRI,*) 'E2: ',E2 1647 WRITE (LUPRI,*) 'EP0:',EP0 1648 WRITE (LUPRI,*) 'EP1:',EP1 1649 WRITE (LUPRI,*) 'EPP0:',EPP0 1650 WRITE (LUPRI,*) 'EPP1:',EPP1 1651 WRITE (LUPRI,*) 'EPP2:',EPP2 1652 WRITE (LUPRI,*) 'EPPP0:',EPPP0 1653 WRITE (LUPRI,*) 'EPPP1:',EPPP1 1654 WRITE (LUPRI,*) 'EPPP2:',EPPP2 1655 WRITE (LUPRI,*) 'EPPP3:',EPPP3 1656 CALL QUIT('ERROR DURING CALCULATION OF E_23 COEFFICIENTS.') 1657 END IF 1658 ABCDE(KE,2) = E0 1659 ABCDE(KEP,2) = EP0 1660 ABCDE(KEP2,2) = EPP0 1661 ABCDE(KEP3,2) = EPPP0 1662 IF (LOCDBG) THEN 1663 WRITE (LUPRI,'(A,4G16.8)') 1664 & " E_23, E_23', E_23'', E_23''' coefficients:", 1665 & E0, EP0, EPP0, EPPP0 1666 END IF 1667 END IF 1668 1669 END IF 1670 1671 RETURN 1672 END 1673*---------------------------------------------------------------------* 1674* END OF SUBROUTINE CCQRDISP * 1675*---------------------------------------------------------------------* 1676 1677c /* deck ccqhypprt */ 1678*=====================================================================* 1679 SUBROUTINE CCQHYPPRT(HYPERPOL) 1680*---------------------------------------------------------------------* 1681* 1682* Purpose: print frequency-dependent first hyperpolarizabilities 1683* 1684* 1685* Written by Christof Haettig in October 1996. 1686* 1687*=====================================================================* 1688#if defined (IMPLICIT_NONE) 1689 IMPLICIT NONE 1690#else 1691# include "implicit.h" 1692#endif 1693#include "priunit.h" 1694#include "ccorb.h" 1695#include "ccsdinp.h" 1696#include "ccqrinf.h" 1697#include "ccroper.h" 1698 1699 1700 CHARACTER*15 BLANKS 1701 CHARACTER*10 RLXA, RLXB, RLXC 1702 CHARACTER*80 STRING 1703 CHARACTER*8 BETALAB 1704 INTEGER ISYMA, ISYMB, ISYMC, ISAMA, ISAMB, ISAMC, ISAPROP 1705 INTEGER IFREQ, IOPER, ICMP, I, J, K, L 1706 1707 1708#if defined (SYS_CRAY) 1709 REAL HYPERPOL(NQRFREQ,NQROPER,2) 1710 REAL HALF, ONE, SIGN, BETA, PROP, ERROR, THIRD, ZERO 1711 REAL FACTOR 1712 REAL BETANEW(NQRFREQ,3,3,3) 1713 REAL MEANBETA, MUBETANORM, MU_DOT_BETA 1714 REAL BETAVEC(3), MUBETA(3) 1715#else 1716 DOUBLE PRECISION HYPERPOL(NQRFREQ,NQROPER,2) 1717 DOUBLE PRECISION HALF, ONE, SIGN, BETA, PROP, ERROR, THIRD, ZERO 1718 DOUBLE PRECISION FACTOR 1719 DOUBLE PRECISION BETANEW(NQRFREQ,3,3,3) 1720 DOUBLE PRECISION MEANBETA, MUBETANORM, MU_DOT_BETA 1721 DOUBLE PRECISION BETAVEC(3), MUBETA(3) 1722#endif 1723 PARAMETER (HALF = 0.5d0, ONE = 1.0d0, THIRD = 1.0d0/3.0d0) 1724 PARAMETER (ZERO = 0.0d0, FACTOR = 1.0d0/6.0d0) 1725 1726 CHARACTER*10 MODEL,MOPRPC 1727 INTEGER ISYMAB,ISYMABC 1728*---------------------------------------------------------------------* 1729* print header for first hyperpolarizabilities: 1730*---------------------------------------------------------------------* 1731 BLANKS = ' ' 1732 STRING = ' RESULTS FOR THE FIRST HYPERPOLARIZABILITIES ' 1733 1734 IF (CCS) THEN 1735 CALL AROUND120( BLANKS//'FINAL CCS'//STRING(1:45)//BLANKS ) 1736 ELSE IF (CC2) THEN 1737 CALL AROUND120( BLANKS//'FINAL CC2'//STRING(1:45)//BLANKS ) 1738 ELSE IF (CCSD) THEN 1739 CALL AROUND120( BLANKS//'FINAL CCSD'//STRING(1:45)//BLANKS ) 1740 ELSE IF (CC3) THEN 1741 CALL AROUND120( BLANKS//'FINAL CC3'//STRING(1:45)//BLANKS ) 1742 ELSE 1743 CALL QUIT('CCQHYPPRT called for an unknown '// 1744 & 'Coupled Cluster model.') 1745 END IF 1746 1747 IF (CCSD) THEN 1748 MODEL = 'CCSD' 1749 MOPRPC = 'CCSD ' 1750 ELSE IF (CIS) THEN 1751 MODEL = 'CIS' 1752 MOPRPC = 'CIS ' 1753 ELSE IF (CCS) THEN 1754 MODEL = 'CCS' 1755 MOPRPC = 'CCS ' 1756 ELSE IF (CC2) THEN 1757 MODEL = 'CC2' 1758 MOPRPC = 'CC2 ' 1759 ELSE IF (CC3) THEN 1760 MODEL = 'CC3' 1761 MOPRPC = 'CC3 ' 1762 END IF 1763 1764 IF (IPRQHYP.GT.5) THEN 1765 WRITE(LUPRI,'(/1X,3(1X,A," operator",17X),4X,A,6X,A,/,130("-"))') 1766 & 'A','B','C','property ','(asy./sym. Resp.)' 1767 ELSE 1768 WRITE(LUPRI,'(/1X,3(1X,A," operator",17X),4X,A,/,102("-"))') 1769 & 'A','B','C','property ' 1770 END IF 1771 1772 1773 DO IOPER = 1, NQROPER 1774 ISYMA = ISYOPR(IAQROP(IOPER)) 1775 ISYMB = ISYOPR(IBQROP(IOPER)) 1776 ISYMC = ISYOPR(ICQROP(IOPER)) 1777 1778 ISAMA = ISYMAT(IAQROP(IOPER)) 1779 ISAMB = ISYMAT(IBQROP(IOPER)) 1780 ISAMC = ISYMAT(ICQROP(IOPER)) 1781 1782 IF ( LAQLRX(IOPER) ) RLXA = ' (relax.) ' 1783 IF ( .NOT. LAQLRX(IOPER) ) RLXA = ' (unrel.) ' 1784 IF ( LBQLRX(IOPER) ) RLXB = ' (relax.) ' 1785 IF ( .NOT. LBQLRX(IOPER) ) RLXB = ' (unrel.) ' 1786 IF ( LCQLRX(IOPER) ) RLXC = ' (relax.) ' 1787 IF ( .NOT. LCQLRX(IOPER) ) RLXC = ' (unrel.) ' 1788 1789 ISAPROP = ISAMA * ISAMB * ISAMC 1790 SIGN = DBLE(ISAPROP) 1791 IF (ISAPROP.EQ.0) SIGN = +ONE 1792 1793 IFREQ = 1 1794cmbh 1795 PROP=0.0d0 1796cmbh end 1797 IF (MULD2H(ISYMA,ISYMB).EQ.ISYMC) THEN 1798 1799 PROP = -HALF*(HYPERPOL(IFREQ,IOPER,1) + 1800 & SIGN * HYPERPOL(IFREQ,IOPER,2)) 1801 ERROR = -HALF*(HYPERPOL(IFREQ,IOPER,1) - 1802 & SIGN * HYPERPOL(IFREQ,IOPER,2)) 1803 1804 IF (IPRQHYP.GT.5 .OR. ISAPROP.EQ.0) THEN 1805 IF (ISAPROP.NE.0) THEN 1806 WRITE(LUPRI, 1807 & '(/1X,3(A8,A10,F7.4,3X),G18.10," (",G18.10,")")') 1808 & LBLOPR(IAQROP(IOPER)),RLXA,AQRFR(IFREQ), 1809 & LBLOPR(IBQROP(IOPER)),RLXB,BQRFR(IFREQ), 1810 & LBLOPR(ICQROP(IOPER)),RLXC,CQRFR(IFREQ), PROP, ERROR 1811 ELSE 1812 WRITE(LUPRI,'(/1X,2A)') 1813 & 'Cannot determine if real or imaginary property... ', 1814 & 'print: symmetric (antisymmetric) parts in +/- w' 1815 WRITE(LUPRI, 1816 & '(1X,3(A8,A10,F7.4,3X),G18.10," (",G18.10,")")') 1817 & LBLOPR(IAQROP(IOPER)),RLXA,AQRFR(IFREQ), 1818 & LBLOPR(IBQROP(IOPER)),RLXB,BQRFR(IFREQ), 1819 & LBLOPR(ICQROP(IOPER)),RLXC,CQRFR(IFREQ), PROP, ERROR 1820 END IF 1821 ELSE 1822 WRITE(LUPRI,'(/1X,3(A8,A10,F7.4,3X),G16.8)') 1823CCN WRITE(LUPRI,'(/1X,3(A8,A10,F7.4,3X),F24.20)') 1824 & LBLOPR(IAQROP(IOPER)),RLXA,AQRFR(IFREQ), 1825 & LBLOPR(IBQROP(IOPER)),RLXB,BQRFR(IFREQ), 1826 & LBLOPR(ICQROP(IOPER)),RLXC,CQRFR(IFREQ), PROP 1827 ENDIF 1828 ELSE 1829 IF (IPRQHYP.GT.5 .OR. ISAPROP.EQ.0) THEN 1830 WRITE(LUPRI, 1831 & '(/1X,3(A8,A10,A7,3X),7X,A,8X," (",9X,A,10X,")")') 1832 & LBLOPR(IAQROP(IOPER)),RLXA,' -.- ', 1833 & LBLOPR(IBQROP(IOPER)),RLXB,' -.- ', 1834 & LBLOPR(ICQROP(IOPER)),RLXC,' -.- ', 1835 & '---', 1836 & '---' 1837 ELSE IF (IPRQHYP.GT.0) THEN 1838 WRITE(LUPRI,'(/1X,3(A8,A10,A7,3X),6X,A,7X)') 1839 & LBLOPR(IAQROP(IOPER)),RLXA,' -.- ', 1840 & LBLOPR(IBQROP(IOPER)),RLXB,' -.- ', 1841 & LBLOPR(ICQROP(IOPER)),RLXC,' -.- ', 1842 & '---' 1843 END IF 1844 END IF 1845 ISYMAB = MULD2H(ISYMA,ISYMB) 1846 ISYMABC = MULD2H(ISYMAB,ISYMC) 1847 CALL WRIPRO(PROP,MOPRPC,3, 1848 & LBLOPR(IAQROP(IOPER)),LBLOPR(IBQROP(IOPER)), 1849 & LBLOPR(ICQROP(IOPER)),LBLOPR(ICQROP(IOPER)), 1850 & BQRFR(IFREQ),CQRFR(IFREQ),CQRFR(IFREQ),ISYMABC, 1851 & 0,0,0) 1852 DO IFREQ = 2, NQRFREQ 1853 IF (MULD2H(ISYMA,ISYMB).EQ.ISYMC) THEN 1854 1855 PROP = -HALF*(HYPERPOL(IFREQ,IOPER,1) + 1856 & SIGN * HYPERPOL(IFREQ,IOPER,2)) 1857 ERROR = -HALF*(HYPERPOL(IFREQ,IOPER,1) - 1858 & SIGN * HYPERPOL(IFREQ,IOPER,2)) 1859 1860 IF (IPRQHYP.GT.5 .OR. ISAPROP.EQ.0) THEN 1861 WRITE(LUPRI,'(1X,3(18X,F7.4,3X),G18.10," (",G18.10,")")') 1862 & AQRFR(IFREQ), BQRFR(IFREQ), CQRFR(IFREQ), PROP, ERROR 1863 ELSE 1864 WRITE(LUPRI,'(1X,3(18X,F7.4,3X),G16.8)') 1865 & AQRFR(IFREQ), BQRFR(IFREQ), CQRFR(IFREQ), PROP 1866 END IF 1867 END IF 1868 CALL WRIPRO(PROP,MOPRPC,3, 1869 & LBLOPR(IAQROP(IOPER)),LBLOPR(IBQROP(IOPER)), 1870 & LBLOPR(ICQROP(IOPER)),LBLOPR(ICQROP(IOPER)), 1871 & BQRFR(IFREQ),CQRFR(IFREQ),CQRFR(IFREQ),ISYMABC, 1872 & 0,0,0) 1873 END DO 1874 1875 END DO 1876 1877 IF (BETA_AVERAGE) THEN 1878 WRITE(LUPRI,'(/////A,/,55("-"))')' average frequencies' 1879 1880* calculate & print beta_{||}: 1881 IFREQ = 1 1882 BETA = -0.3d0*(HYPERPOL(IFREQ,1,1)+HYPERPOL(IFREQ,1,2)) 1883 IF (XY_DEGENERAT) THEN 1884 DO ICMP = 2, 4 1885 BETA = BETA- 1886 & 0.2d0*(HYPERPOL(IFREQ,ICMP,1)+HYPERPOL(IFREQ,ICMP,2)) 1887 END DO 1888 ELSE 1889 DO ICMP = 2, 7 1890 BETA = BETA- 1891 & 0.1d0*(HYPERPOL(IFREQ,ICMP,1)+HYPERPOL(IFREQ,ICMP,2)) 1892 END DO 1893 END IF 1894 1895 WRITE(LUPRI,'(/1X,A8,3X,3(F7.4,2X),G16.8)') 1896 & ' beta_||', AQRFR(IFREQ), BQRFR(IFREQ), CQRFR(IFREQ), BETA 1897C 1898 BETALAB='beta_|| ' 1899 CALL WRIPRO(BETA,MOPRPC,3, 1900 * BETALAB,BETALAB,BETALAB,BETALAB, 1901 * AQRFR(IFREQ),BQRFR(IFREQ),CQRFR(IFREQ),ISYMABC, 1902 * 0,0,0) 1903C 1904 DO IFREQ = 2, NQRFREQ 1905 BETA = -0.3d0*(HYPERPOL(IFREQ,1,1)+HYPERPOL(IFREQ,1,2)) 1906 IF (XY_DEGENERAT) THEN 1907 DO ICMP = 2, 4 1908 BETA = BETA- 1909 & 0.2d0*(HYPERPOL(IFREQ,ICMP,1)+HYPERPOL(IFREQ,ICMP,2)) 1910 END DO 1911 ELSE 1912 DO ICMP = 2, 7 1913 BETA = BETA- 1914 & 0.1d0*(HYPERPOL(IFREQ,ICMP,1)+HYPERPOL(IFREQ,ICMP,2)) 1915 END DO 1916 END IF 1917 1918 WRITE(LUPRI,'(1X,A8,3X,3(F7.4,2X),G16.8)') 1919 & ' ', AQRFR(IFREQ), BQRFR(IFREQ), CQRFR(IFREQ), BETA 1920C 1921 BETALAB='beta_|| ' 1922 CALL WRIPRO(BETA,MOPRPC,3, 1923 * BETALAB,BETALAB,BETALAB,BETALAB, 1924 * AQRFR(IFREQ),BQRFR(IFREQ),CQRFR(IFREQ),ISYMABC, 1925 * 0,0,0) 1926C 1927 END DO 1928 1929* calculate & print beta^K: 1930 IFREQ = 1 1931 BETA = -0.3d0*(HYPERPOL(IFREQ,1,1)+HYPERPOL(IFREQ,1,2)) 1932 BETA = BETA + 0.15d0*(HYPERPOL(IFREQ,2,1)+HYPERPOL(IFREQ,2,2)) 1933 BETA = BETA - 0.60d0*(HYPERPOL(IFREQ,3,1)+HYPERPOL(IFREQ,3,2)) 1934 BETA = BETA + 0.15d0*(HYPERPOL(IFREQ,4,1)+HYPERPOL(IFREQ,4,2)) 1935 IF (XY_DEGENERAT) THEN 1936 BETA = BETA + 0.15d0*(HYPERPOL(IFREQ,2,1)+HYPERPOL(IFREQ,2,2)) 1937 BETA = BETA - 0.60d0*(HYPERPOL(IFREQ,3,1)+HYPERPOL(IFREQ,3,2)) 1938 BETA = BETA + 0.15d0*(HYPERPOL(IFREQ,4,1)+HYPERPOL(IFREQ,4,2)) 1939 ELSE 1940 BETA = BETA + 0.15d0*(HYPERPOL(IFREQ,5,1)+HYPERPOL(IFREQ,5,2)) 1941 BETA = BETA - 0.60d0*(HYPERPOL(IFREQ,6,1)+HYPERPOL(IFREQ,6,2)) 1942 BETA = BETA + 0.15d0*(HYPERPOL(IFREQ,7,1)+HYPERPOL(IFREQ,7,2)) 1943 END IF 1944 1945 WRITE(LUPRI,'(/1X,A8,3X,3(F7.4,2X),G16.8)') 1946 & ' beta^K ', AQRFR(IFREQ), BQRFR(IFREQ), CQRFR(IFREQ), BETA 1947C 1948 BETALAB='beta^K ' 1949 CALL WRIPRO(BETA,MOPRPC,3, 1950 * BETALAB,BETALAB,BETALAB,BETALAB, 1951 * AQRFR(IFREQ),BQRFR(IFREQ),CQRFR(IFREQ),ISYMABC, 1952 * 0,0,0) 1953C 1954 DO IFREQ = 2, NQRFREQ 1955 BETA = -0.3d0*(HYPERPOL(IFREQ,1,1)+HYPERPOL(IFREQ,1,2)) 1956 BETA = BETA + 0.15d0*(HYPERPOL(IFREQ,2,1)+HYPERPOL(IFREQ,2,2)) 1957 BETA = BETA - 0.60d0*(HYPERPOL(IFREQ,3,1)+HYPERPOL(IFREQ,3,2)) 1958 BETA = BETA + 0.15d0*(HYPERPOL(IFREQ,4,1)+HYPERPOL(IFREQ,4,2)) 1959 IF (XY_DEGENERAT) THEN 1960 BETA = BETA + 0.15d0*(HYPERPOL(IFREQ,2,1)+HYPERPOL(IFREQ,2,2)) 1961 BETA = BETA - 0.60d0*(HYPERPOL(IFREQ,3,1)+HYPERPOL(IFREQ,3,2)) 1962 BETA = BETA + 0.15d0*(HYPERPOL(IFREQ,4,1)+HYPERPOL(IFREQ,4,2)) 1963 ELSE 1964 BETA = BETA + 0.15d0*(HYPERPOL(IFREQ,5,1)+HYPERPOL(IFREQ,5,2)) 1965 BETA = BETA - 0.60d0*(HYPERPOL(IFREQ,6,1)+HYPERPOL(IFREQ,6,2)) 1966 BETA = BETA + 0.15d0*(HYPERPOL(IFREQ,7,1)+HYPERPOL(IFREQ,7,2)) 1967 END IF 1968 1969 WRITE(LUPRI,'(1X,A8,3X,3(F7.4,2X),G16.8)') 1970 & ' ', AQRFR(IFREQ), BQRFR(IFREQ), CQRFR(IFREQ), BETA 1971C 1972 BETALAB='beta^K ' 1973 CALL WRIPRO(BETA,MOPRPC,3, 1974 * BETALAB,BETALAB,BETALAB,BETALAB, 1975 * AQRFR(IFREQ),BQRFR(IFREQ),CQRFR(IFREQ),ISYMABC, 1976 * 0,0,0) 1977C 1978 END DO 1979 1980 1981* calculate & print beta_{_|_}: 1982 IFREQ = 1 1983 BETA = -0.1d0*(HYPERPOL(IFREQ,1,1)+HYPERPOL(IFREQ,1,2)) 1984 BETA = BETA - 0.2d0*(HYPERPOL(IFREQ,2,1)+HYPERPOL(IFREQ,2,2)) 1985 BETA = BETA + 0.3d0*(HYPERPOL(IFREQ,3,1)+HYPERPOL(IFREQ,3,2)) 1986 BETA = BETA - 0.2d0*(HYPERPOL(IFREQ,4,1)+HYPERPOL(IFREQ,4,2)) 1987 IF (XY_DEGENERAT) THEN 1988 BETA = BETA - 0.2d0*(HYPERPOL(IFREQ,2,1)+HYPERPOL(IFREQ,2,2)) 1989 BETA = BETA + 0.3d0*(HYPERPOL(IFREQ,3,1)+HYPERPOL(IFREQ,3,2)) 1990 BETA = BETA - 0.2d0*(HYPERPOL(IFREQ,4,1)+HYPERPOL(IFREQ,4,2)) 1991 ELSE 1992 BETA = BETA - 0.2d0*(HYPERPOL(IFREQ,5,1)+HYPERPOL(IFREQ,5,2)) 1993 BETA = BETA + 0.3d0*(HYPERPOL(IFREQ,6,1)+HYPERPOL(IFREQ,6,2)) 1994 BETA = BETA - 0.2d0*(HYPERPOL(IFREQ,7,1)+HYPERPOL(IFREQ,7,2)) 1995 END IF 1996 1997 WRITE(LUPRI,'(/1X,A8,3X,3(F7.4,2X),G16.8)') 1998 & ' beta_|_', AQRFR(IFREQ), BQRFR(IFREQ), CQRFR(IFREQ), BETA 1999C 2000 BETALAB='beta_|_ ' 2001 CALL WRIPRO(BETA,MOPRPC,3, 2002 * BETALAB,BETALAB,BETALAB,BETALAB, 2003 * AQRFR(IFREQ),BQRFR(IFREQ),CQRFR(IFREQ),ISYMABC, 2004 * 0,0,0) 2005C 2006 DO IFREQ = 2, NQRFREQ 2007 BETA = -0.1d0*(HYPERPOL(IFREQ,1,1)+HYPERPOL(IFREQ,1,2)) 2008 BETA = BETA - 0.2d0*(HYPERPOL(IFREQ,2,1)+HYPERPOL(IFREQ,2,2)) 2009 BETA = BETA + 0.3d0*(HYPERPOL(IFREQ,3,1)+HYPERPOL(IFREQ,3,2)) 2010 BETA = BETA - 0.2d0*(HYPERPOL(IFREQ,4,1)+HYPERPOL(IFREQ,4,2)) 2011 IF (XY_DEGENERAT) THEN 2012 BETA = 2.0d0 * BETA 2013 ELSE 2014 BETA = BETA - 0.2d0*(HYPERPOL(IFREQ,5,1)+HYPERPOL(IFREQ,5,2)) 2015 BETA = BETA + 0.3d0*(HYPERPOL(IFREQ,6,1)+HYPERPOL(IFREQ,6,2)) 2016 BETA = BETA - 0.2d0*(HYPERPOL(IFREQ,7,1)+HYPERPOL(IFREQ,7,2)) 2017 END IF 2018 2019 WRITE(LUPRI,'(1X,A8,3X,3(F7.4,2X),G16.8)') 2020 & ' ', AQRFR(IFREQ), BQRFR(IFREQ), CQRFR(IFREQ), BETA 2021C 2022 BETALAB='beta_|_ ' 2023 CALL WRIPRO(BETA,MOPRPC,3, 2024 * BETALAB,BETALAB,BETALAB,BETALAB, 2025 * AQRFR(IFREQ),BQRFR(IFREQ),CQRFR(IFREQ),ISYMABC, 2026 * 0,0,0) 2027C 2028 END DO 2029 2030* calculate & print beta_ms: 2031 IFREQ = 1 2032 BETA = 0.0d0 2033 BETA = BETA - 0.25d0*(HYPERPOL(IFREQ,2,1)+HYPERPOL(IFREQ,2,2)) 2034 BETA = BETA + 0.50d0*(HYPERPOL(IFREQ,3,1)+HYPERPOL(IFREQ,3,2)) 2035 BETA = BETA - 0.25d0*(HYPERPOL(IFREQ,4,1)+HYPERPOL(IFREQ,4,2)) 2036 IF (XY_DEGENERAT) THEN 2037 BETA = BETA - 0.25d0*(HYPERPOL(IFREQ,2,1)+HYPERPOL(IFREQ,2,2)) 2038 BETA = BETA + 0.50d0*(HYPERPOL(IFREQ,3,1)+HYPERPOL(IFREQ,3,2)) 2039 BETA = BETA - 0.25d0*(HYPERPOL(IFREQ,4,1)+HYPERPOL(IFREQ,4,2)) 2040 ELSE 2041 BETA = BETA - 0.25d0*(HYPERPOL(IFREQ,5,1)+HYPERPOL(IFREQ,5,2)) 2042 BETA = BETA + 0.50d0*(HYPERPOL(IFREQ,6,1)+HYPERPOL(IFREQ,6,2)) 2043 BETA = BETA - 0.25d0*(HYPERPOL(IFREQ,7,1)+HYPERPOL(IFREQ,7,2)) 2044 END IF 2045 2046 WRITE(LUPRI,'(/1X,A11,3(F7.4,2X),G16.8)') 2047 & ' beta_ms ', AQRFR(IFREQ), BQRFR(IFREQ), CQRFR(IFREQ), BETA 2048C 2049 BETALAB='beta_ms ' 2050 CALL WRIPRO(BETA,MOPRPC,3, 2051 * BETALAB,BETALAB,BETALAB,BETALAB, 2052 * AQRFR(IFREQ),BQRFR(IFREQ),CQRFR(IFREQ),ISYMABC, 2053 * 0,0,0) 2054C 2055 2056 2057 DO IFREQ = 2, NQRFREQ 2058 BETA = 0.0d0 2059 BETA = BETA - 0.25d0*(HYPERPOL(IFREQ,2,1)+HYPERPOL(IFREQ,2,2)) 2060 BETA = BETA + 0.50d0*(HYPERPOL(IFREQ,3,1)+HYPERPOL(IFREQ,3,2)) 2061 BETA = BETA - 0.25d0*(HYPERPOL(IFREQ,4,1)+HYPERPOL(IFREQ,4,2)) 2062 IF (XY_DEGENERAT) THEN 2063 BETA = 2.0d0 * BETA 2064 ELSE 2065 BETA = BETA - 0.25d0*(HYPERPOL(IFREQ,5,1)+HYPERPOL(IFREQ,5,2)) 2066 BETA = BETA + 0.50d0*(HYPERPOL(IFREQ,6,1)+HYPERPOL(IFREQ,6,2)) 2067 BETA = BETA - 0.25d0*(HYPERPOL(IFREQ,7,1)+HYPERPOL(IFREQ,7,2)) 2068 END IF 2069 2070 WRITE(LUPRI,'(1X,A8,3X,3(F7.4,2X),G16.8)') 2071 & ' ', AQRFR(IFREQ), BQRFR(IFREQ), CQRFR(IFREQ), BETA 2072C 2073 BETALAB='beta_ms ' 2074 CALL WRIPRO(BETA,MOPRPC,3, 2075 * BETALAB,BETALAB,BETALAB,BETALAB, 2076 * AQRFR(IFREQ),BQRFR(IFREQ),CQRFR(IFREQ),ISYMABC, 2077 * 0,0,0) 2078C 2079 END DO 2080 END IF 2081C----------------------------------------------------------------------- 2082C June 04: AO+JK 2083C .AVANEW in input block *CCQR -> LAVANEW =.TRUE. 2084C Calculates: 2085C i) Beta_i; i=x,y,z 2086C ii) Mu_i*Beta_i -> |Mu*Beta| 2087C iii) <Beta> 2088C----------------------------------------------------------------------- 2089C 2090 IF (LAVANEW) THEN 2091C 2092 BLANKS = ' ' 2093 STRING = ' AVERAGES FOR SECOND HARMONIC GENERATION ' 2094C 2095 IF (CCS) THEN 2096 CALL AROUND120( BLANKS//'FINAL CCS'//STRING(1:45)//BLANKS ) 2097 ELSE IF (CC2) THEN 2098 CALL AROUND120( BLANKS//'FINAL CC2'//STRING(1:45)//BLANKS ) 2099 ELSE IF (CCSD) THEN 2100 CALL AROUND120( BLANKS//'FINAL CCSD'//STRING(1:45)//BLANKS ) 2101 ELSE 2102 CALL QUIT('CCQHYPPRT called for an unknown '// 2103 & 'Coupled Cluster model.') 2104 END IF 2105C 2106C---------------------------------- 2107C Reset BETANEW(IFREQ,I,J,K): 2108C---------------------------------- 2109C 2110 DO IFREQ=1,NQRFREQ 2111 DO I=1,3 2112 DO J=1,3 2113 DO K=1,3 2114 BETANEW(IFREQ,I,J,K) = ZERO 2115 END DO 2116 END DO 2117 END DO 2118 END DO 2119C 2120C--------------------------------------- 2121C Rearrange the components of beta 2122C--------------------------------------- 2123C 2124 WRITE(LUPRI,'(/////20X,A)') 2125 * '+----------+----------+----------+----------+----------+' 2126 WRITE(LUPRI,'(20X,A)') 2127 * '| Property:| Freq_A: | Freq_B: | Freq_C: | Value: |' 2128 WRITE(LUPRI,'(20X,A)') 2129 * '|----------+----------+----------+----------+----------|' 2130 WRITE(LUPRI,'(20X,A,A,A,F8.4,A,F8.4,A,F8.4,A,F8.4,A)') 2131 * '| ','Dipole_x',' | ',ZERO,' | ',ZERO,' | ',ZERO,' | ', 2132 * DIPSAVE(1),' |' 2133 WRITE(LUPRI,'(20X,A,A,A,F8.4,A,F8.4,A,F8.4,A,F8.4,A)') 2134 * '| ','Dipole_y',' | ',ZERO,' | ',ZERO,' | ',ZERO,' | ', 2135 * DIPSAVE(2),' |' 2136 WRITE(LUPRI,'(20X,A,A,A,F8.4,A,F8.4,A,F8.4,A,F8.4,A)') 2137 * '| ','Dipole_z',' | ',ZERO,' | ',ZERO,' | ',ZERO,' | ', 2138 * DIPSAVE(3),' |' 2139 WRITE(LUPRI,'(20X,A)') 2140 * '|----------+----------+----------+----------+----------|' 2141C 2142 L = 0 2143 DO I = 1, 3 2144 DO J = 1, 3 2145 DO K = 1, 3 2146 L = L + 1 2147 ISAMA = ISYMAT(IAQROP(L)) 2148 ISAMB = ISYMAT(IBQROP(L)) 2149 ISAMC = ISYMAT(ICQROP(L)) 2150 ISYMA = ISYOPR(IAQROP(L)) 2151 ISYMB = ISYOPR(IBQROP(L)) 2152 ISYMC = ISYOPR(ICQROP(L)) 2153 ISYMAB = MULD2H(ISYMA,ISYMB) 2154 ISYMABC = MULD2H(ISYMAB,ISYMC) 2155 ISAPROP = ISAMA * ISAMB * ISAMC 2156 SIGN = DBLE(ISAPROP) 2157 IF (ISAPROP .EQ. 0) SIGN = +ONE 2158 DO IFREQ = 1, NQRFREQ 2159 BETANEW(IFREQ,I,J,K) = -HALF*(HYPERPOL(IFREQ,L,1) 2160 * + SIGN * HYPERPOL(IFREQ,L,2)) 2161 END DO 2162 END DO 2163 END DO 2164 END DO 2165C 2166C 2167C---------------------------------------------------------------- 2168C Calculate the components of the beta_i; i=x,y,z vector: 2169C---------------------------------------------------------------- 2170C 2171 DO IFREQ = 1, NQRFREQ 2172 MUBETANORM = ZERO 2173 MU_DOT_BETA= ZERO 2174 MEANBETA = ZERO 2175 DO I = 1, 3 2176 BETAVEC(I) = ZERO 2177 MUBETA(I) = ZERO 2178 DO J = 1, 3 2179 BETAVEC(I) = BETAVEC(I) + THIRD * (BETANEW(IFREQ,I,J,J) + 2180 * BETANEW(IFREQ,J,J,I) + BETANEW(IFREQ,J,I,J)) 2181 END DO 2182 MUBETA(I) = DIPSAVE(I) * BETAVEC(I) 2183C 2184 IF (I .EQ. 1) BETALAB='beta_x ' 2185 IF (I .EQ. 2) BETALAB='beta_y ' 2186 IF (I .EQ. 3) BETALAB='beta_z ' 2187C 2188 WRITE(LUPRI,'(20X,A,A,A,F8.4,A,F8.4,A,F8.4,A,F8.4,A)') 2189 * '| ',BETALAB,' | ',AQRFR(IFREQ),' | ', 2190 * BQRFR(IFREQ),' | ',CQRFR(IFREQ),' | ', 2191 * BETAVEC(I),' |' 2192C 2193 CALL WRIPRO(BETAVEC(I),MOPRPC,3, 2194 * BETALAB,BETALAB,BETALAB,BETALAB, 2195 * AQRFR(IFREQ),BQRFR(IFREQ),CQRFR(IFREQ),ISYMABC, 2196 * 0,0,0) 2197C 2198 END DO 2199C 2200 WRITE(LUPRI,'(20X,A)') 2201 * '+----------+----------+----------+----------+----------+' 2202C 2203 DO I = 1, 3 2204 MUBETANORM = MUBETANORM + MUBETA(I)*MUBETA(I) 2205 MU_DOT_BETA= MU_DOT_BETA+ MUBETA(I) 2206 IF (I .EQ. 1) BETALAB='mu*bet_x' 2207 IF (I .EQ. 2) BETALAB='mu*bet_y' 2208 IF (I .EQ. 3) BETALAB='mu*bet_z' 2209 WRITE(LUPRI,'(20X,A,A,A,F8.4,A,F8.4,A,F8.4,A,F8.4,A)') 2210 * '| ',BETALAB,' | ',AQRFR(IFREQ),' | ', 2211 * BQRFR(IFREQ),' | ',CQRFR(IFREQ),' | ', 2212 * MUBETA(I),' |' 2213C 2214 CALL WRIPRO(MUBETA(I),MOPRPC,3, 2215 * BETALAB,BETALAB,BETALAB,BETALAB, 2216 * AQRFR(IFREQ),BQRFR(IFREQ),CQRFR(IFREQ),ISYMABC, 2217 * 0,0,0) 2218 END DO 2219C 2220 WRITE(LUPRI,'(20X,A)') 2221 * '+----------+----------+----------+----------+----------+' 2222C 2223 BETALAB = ' mu*bet ' 2224 WRITE(LUPRI,'(20X,A,A,A,F8.4,A,F8.4,A,F8.4,A,F8.4,A)') 2225 * '| ',BETALAB,' | ',AQRFR(IFREQ),' | ', 2226 * BQRFR(IFREQ),' | ',CQRFR(IFREQ),' | ', 2227 * MU_DOT_BETA,' |' 2228C 2229 CALL WRIPRO(MU_DOT_BETA,MOPRPC,3, 2230 * BETALAB,BETALAB,BETALAB,BETALAB, 2231 * AQRFR(IFREQ),BQRFR(IFREQ),CQRFR(IFREQ),ISYMABC, 2232 * 0,0,0) 2233C 2234 WRITE(LUPRI,'(20X,A)') 2235 * '+----------+----------+----------+----------+----------+' 2236 MUBETANORM = SQRT(MUBETANORM) 2237 BETALAB = '|mu*bet|' 2238 WRITE(LUPRI,'(20X,A,A,A,F8.4,A,F8.4,A,F8.4,A,F8.4,A)') 2239 * '| ',BETALAB,' | ',AQRFR(IFREQ),' | ', 2240 * BQRFR(IFREQ),' | ',CQRFR(IFREQ),' | ', 2241 * MUBETANORM,' |' 2242C 2243 CALL WRIPRO(MUBETANORM,MOPRPC,3, 2244 * BETALAB,BETALAB,BETALAB,BETALAB, 2245 * AQRFR(IFREQ),BQRFR(IFREQ),CQRFR(IFREQ),ISYMABC, 2246 * 0,0,0) 2247C 2248 WRITE(LUPRI,'(20X,A)') 2249 * '+----------+----------+----------+----------+----------+' 2250 MEANBETA = FACTOR * ( BETANEW(IFREQ,1,2,3) - 2251 * BETANEW(IFREQ,1,3,2) + BETANEW(IFREQ,2,3,1) - 2252 * BETANEW(IFREQ,2,1,3) + BETANEW(IFREQ,3,1,2) - 2253 * BETANEW(IFREQ,3,2,1) ) 2254 BETALAB = '<beta> ' 2255 WRITE(LUPRI,'(20X,A,A,A,F8.4,A,F8.4,A,F8.4,A,F8.4,A)') 2256 * '| ',BETALAB,' | ',AQRFR(IFREQ),' | ', 2257 * BQRFR(IFREQ),' | ',CQRFR(IFREQ),' | ', 2258 * MEANBETA,' |' 2259C 2260 CALL WRIPRO(MEANBETA,MOPRPC,3, 2261 * BETALAB,BETALAB,BETALAB,BETALAB, 2262 * AQRFR(IFREQ),BQRFR(IFREQ),CQRFR(IFREQ),ISYMABC, 2263 * 0,0,0) 2264C 2265 WRITE(LUPRI,'(20X,A)') 2266 * '+----------+----------+----------+----------+----------+' 2267 END DO 2268 END IF 2269 2270 WRITE(LUPRI,'(//,102("-"),///)') 2271 2272 RETURN 2273 END 2274*---------------------------------------------------------------------* 2275* END OF SUBROUTINE CCQHYPPRT * 2276*---------------------------------------------------------------------* 2277c /* deck ccqexpprt */ 2278*=====================================================================* 2279 SUBROUTINE CCQEXPPRT(EXPCOF) 2280*---------------------------------------------------------------------* 2281* 2282* Purpose: print expansion coefficients for 2283* first hyperpolarizabilities 2284* 2285* 2286* Written by Christof Haettig in October 1997. 2287* 2288*=====================================================================* 2289#if defined (IMPLICIT_NONE) 2290 IMPLICIT NONE 2291#else 2292# include "implicit.h" 2293#endif 2294#include "priunit.h" 2295#include "ccorb.h" 2296#include "ccsdinp.h" 2297#include "ccqrinf.h" 2298#include "ccroper.h" 2299 2300 2301 LOGICAL TSTISA 2302 CHARACTER*5 BLANKS 2303 CHARACTER*80 STRING 2304 INTEGER ISYMA, ISYMB, ISYMC, ISAMA, ISAMB, ISAMC, ISAPROP 2305 INTEGER IDISP, IOPER, ICTOT, ISACAU, IDISP2 2306 2307 2308#if defined (SYS_CRAY) 2309 REAL EXPCOF(NQRDISP,NQROPER) 2310#else 2311 DOUBLE PRECISION EXPCOF(NQRDISP,NQROPER) 2312#endif 2313 2314*---------------------------------------------------------------------* 2315* print header for expansion coefficients: 2316*---------------------------------------------------------------------* 2317 BLANKS = ' ' 2318 STRING = ' RESULTS FOR EXPANSION COEFFICIENTS' 2319 2320 IF (CCS) THEN 2321 CALL AROUND( BLANKS//'FINAL CCS'//STRING(1:45)//BLANKS ) 2322 ELSE IF (CC2) THEN 2323 CALL AROUND( BLANKS//'FINAL CC2'//STRING(1:45)//BLANKS ) 2324 ELSE IF (CCSD) THEN 2325 CALL AROUND( BLANKS//'FINAL CCSD'//STRING(1:45)//BLANKS ) 2326 ELSE 2327 CALL QUIT('CCQEXPPRT called for an unknown '// 2328 & 'Coupled Cluster model.') 2329 END IF 2330 2331 WRITE(LUPRI,'(/1X,3(1X,A," operator",3X),4X,A,/,68("-"))') 2332 & 'A','B','C','d_ABC' 2333 2334 2335 DO IOPER = 1, NQROPER 2336 ISYMA = ISYOPR(IAQROP(IOPER)) 2337 ISYMB = ISYOPR(IBQROP(IOPER)) 2338 ISYMC = ISYOPR(ICQROP(IOPER)) 2339 2340 ISAMA = ISYMAT(IAQROP(IOPER)) 2341 ISAMB = ISYMAT(IBQROP(IOPER)) 2342 ISAMC = ISYMAT(ICQROP(IOPER)) 2343 2344 ISAPROP = ISAMA * ISAMB * ISAMC 2345 2346 IDISP = 1 2347 2348 ICTOT = IQCAUA(IDISP) + IQCAUB(IDISP) + IQCAUC(IDISP) 2349 ISACAU = 2 * ( (ICTOT/2)*2 - ICTOT ) + 1 2350 TSTISA = ISACAU.EQ.ISAPROP .OR. ISAPROP.EQ.0 .OR. ALLDSPCF 2351 2352 IF (.NOT. TSTISA) IDISP = IDISP + 1 2353 2354 IF (MULD2H(ISYMA,ISYMB).EQ.ISYMC) THEN 2355 WRITE(LUPRI,'(/1X,3(A8,I3,3X),G16.8)') 2356 & LBLOPR(IAQROP(IOPER)),IQCAUA(IDISP), 2357 & LBLOPR(IBQROP(IOPER)),IQCAUB(IDISP), 2358 & LBLOPR(ICQROP(IOPER)),IQCAUC(IDISP), 2359 & -EXPCOF(IDISP,IOPER) 2360 ELSE 2361 IF (IPRQHYP.GT.0) THEN 2362 WRITE(LUPRI,'(/1X,3(A8,A3,3X),6X,A,7X)') 2363 & LBLOPR(IAQROP(IOPER)),' - ', 2364 & LBLOPR(IBQROP(IOPER)),' - ', 2365 & LBLOPR(ICQROP(IOPER)),' - ', 2366 & '---' 2367 END IF 2368 END IF 2369 2370 IDISP2 = IDISP + 1 2371 2372 DO IDISP = IDISP2, NQRDISP 2373 ICTOT = IQCAUA(IDISP) + IQCAUB(IDISP) + IQCAUC(IDISP) 2374 ISACAU = 2 * ( (ICTOT/2)*2 - ICTOT ) + 1 2375 TSTISA = ISACAU.EQ.ISAPROP .OR. ISAPROP.EQ.0 .OR. ALLDSPCF 2376 2377 IF (MULD2H(ISYMA,ISYMB).EQ.ISYMC .AND. TSTISA) THEN 2378 WRITE(LUPRI,'(1X,3(8X,I3,3X),G16.8)') 2379 & IQCAUA(IDISP), IQCAUB(IDISP), IQCAUC(IDISP), 2380 & -EXPCOF(IDISP,IOPER) 2381 END IF 2382 END DO 2383 2384 END DO 2385 2386 RETURN 2387 END 2388*---------------------------------------------------------------------* 2389* END OF SUBROUTINE CCQEXPPRT * 2390*---------------------------------------------------------------------* 2391c /* deck ccqdisprt */ 2392*=====================================================================* 2393 SUBROUTINE CCQDISPRT (DISPCF, AVEDISP,ABCDE, 2394 & SHGCF, ORCOF, EOPECF, 2395 & AVESHG,AVEOR,AVEEOPE) 2396*---------------------------------------------------------------------* 2397* 2398* Purpose: print dispersion coefficients first hyperpolarizabilities 2399* 2400* 2401* Written by Christof Haettig in November 1997. 2402* 2403*=====================================================================* 2404#if defined (IMPLICIT_NONE) 2405 IMPLICIT NONE 2406#else 2407# include "implicit.h" 2408#endif 2409#include "priunit.h" 2410#include "ccorb.h" 2411#include "ccsdinp.h" 2412#include "ccqrinf.h" 2413#include "ccroper.h" 2414 2415 INTEGER KA, KB, KC, KCP, KD, KDP, KE, KEP 2416 PARAMETER (KA=1,KB=2,KC=3,KCP=4,KD=5,KDP=6,KE=7,KEP=8) 2417 INTEGER KBP, KDP2, KEP2, KEP3 2418 PARAMETER (KBP=9,KDP2=10,KEP2=11,KEP3=12) 2419 2420 CHARACTER*5 BLANKS 2421 CHARACTER*80 STRING 2422 INTEGER ISYMA, ISYMB, ISYMC, ISAMA, ISAMB, ISAMC, ISAPROP 2423 INTEGER IOPER, MN, N, M, MN0, MNS 2424 2425 2426#if defined (SYS_CRAY) 2427 REAL DISPCF((NQRDSPE+2)*(NQRDSPE+1)/2,NQROPER) 2428 REAL AVEDISP((NQRDSPE+2)*(NQRDSPE+1)/2,4) 2429 REAL ABCDE(12,2) 2430 REAL SHGCF(NQRDSPE+1,NQROPER) 2431 REAL ORCOF(NQRDSPE+1,NQROPER) 2432 REAL EOPECF(NQRDSPE+1,NQROPER) 2433 REAL AVESHG(NQRDSPE+1,4) 2434 REAL AVEOR(NQRDSPE+1,4) 2435 REAL AVEEOPE(NQRDSPE+1,4) 2436#else 2437 DOUBLE PRECISION DISPCF((NQRDSPE+2)*(NQRDSPE+1)/2,NQROPER) 2438 DOUBLE PRECISION AVEDISP((NQRDSPE+2)*(NQRDSPE+1)/2,4) 2439 DOUBLE PRECISION ABCDE(12,2) 2440 DOUBLE PRECISION SHGCF(NQRDSPE+1,NQROPER) 2441 DOUBLE PRECISION ORCOF(NQRDSPE+1,NQROPER) 2442 DOUBLE PRECISION EOPECF(NQRDSPE+1,NQROPER) 2443 DOUBLE PRECISION AVESHG(NQRDSPE+1,4) 2444 DOUBLE PRECISION AVEOR(NQRDSPE+1,4) 2445 DOUBLE PRECISION AVEEOPE(NQRDSPE+1,4) 2446#endif 2447 2448C 2449 INTEGER ITRI, I, J 2450 ITRI(I,J) = (MAX(I,J)+1)*(MAX(I,J)-2)/2 + I + J + 2 2451C 2452 2453*---------------------------------------------------------------------* 2454* print header for hyperpolarizability dispersion coefficient section 2455*---------------------------------------------------------------------* 2456 BLANKS = ' ' 2457 STRING = ' RESULTS FOR DISPERSION COEFFICIENTS' 2458 2459 IF (CCS) THEN 2460 CALL AROUND( BLANKS//'FINAL CCS'//STRING(1:45)//BLANKS ) 2461 ELSE IF (CC2) THEN 2462 CALL AROUND( BLANKS//'FINAL CC2'//STRING(1:45)//BLANKS ) 2463 ELSE IF (CCSD) THEN 2464 CALL AROUND( BLANKS//'FINAL CCSD'//STRING(1:45)//BLANKS ) 2465 ELSE 2466 CALL QUIT('CCQDISPRT called for an unknown '// 2467 & 'Coupled Cluster model.') 2468 END IF 2469 2470 WRITE(LUPRI,'(/1X,3(A," operator",3X),3(4X,A),/,68("-"))') 2471 & 'A','B','C','N','M','D_ABC' 2472 2473 2474 DO IOPER = 1, NQROPER 2475 ISYMA = ISYOPR(IAQROP(IOPER)) 2476 ISYMB = ISYOPR(IBQROP(IOPER)) 2477 ISYMC = ISYOPR(ICQROP(IOPER)) 2478 2479 ISAMA = ISYMAT(IAQROP(IOPER)) 2480 ISAMB = ISYMAT(IBQROP(IOPER)) 2481 ISAMC = ISYMAT(ICQROP(IOPER)) 2482 2483 ISAPROP = ISAMA * ISAMB * ISAMC 2484 IF (ALLDSPCF) ISAPROP = 0 2485 2486 MN0 = 0 2487 MNS = 2 2488 2489 IF (ISAPROP.EQ.-1) MN0 = 1 2490 IF (ISAPROP.EQ. 0) MNS = 1 2491 2492 MN = MN0 2493 M = 0 2494 N = MN - M 2495 IF (MULD2H(ISYMA,ISYMB).EQ.ISYMC) THEN 2496 WRITE(LUPRI,'(/1X,3(A8,5X),2I5,1X,G16.8)') 2497 & LBLOPR(IAQROP(IOPER)), 2498 & LBLOPR(IBQROP(IOPER)), 2499 & LBLOPR(ICQROP(IOPER)),N,M, 2500 & -DISPCF(ITRI(M+N,N),IOPER) 2501 ELSE 2502 IF (IPRQHYP.GT.0) THEN 2503 WRITE(LUPRI,'(/1X,3(A8,5X),3A)') 2504 & LBLOPR(IAQROP(IOPER)), 2505 & LBLOPR(IBQROP(IOPER)), 2506 & LBLOPR(ICQROP(IOPER)),' - ', ' - ', '---' 2507 END IF 2508 END IF 2509 2510 IF (MN0.EQ.1 .AND. MULD2H(ISYMA,ISYMB).EQ.ISYMC) THEN 2511 M = 1 2512 N = MN - M 2513 WRITE(LUPRI,'(1X,3(8X,5X),2I5,1X,G16.8)') 2514 & N,M, -DISPCF(ITRI(M+N,N),IOPER) 2515 END IF 2516 2517 IF (MULD2H(ISYMA,ISYMB).EQ.ISYMC) THEN 2518 DO MN = MN0+MNS, NQRDSPE, MNS 2519 DO M = 0, MN 2520 N = MN - M 2521 WRITE(LUPRI,'(1X,3(8X,5X),2I5,1X,G16.8)') 2522 & N,M, -DISPCF(ITRI(M+N,N),IOPER) 2523 END DO 2524 END DO 2525 END IF 2526 2527 END DO 2528 2529 IF (BETA_AVERAGE) THEN 2530 MN = 0 2531 N = 0 2532 M = 0 2533 WRITE(LUPRI,'(/1X,A10,29X,2I5,1X,G16.8)') 2534 & 'beta_{||} ',N,M, -AVEDISP(ITRI(M+N,M),1) 2535 DO MN = 2, NQRDSPE, 2 2536 DO M = 0, MN 2537 N = MN - M 2538 WRITE(LUPRI,'(1X,3(8X,5X),2I5,1X,G16.8)') 2539 & N,M, -AVEDISP(ITRI(M+N,M),1) 2540 END DO 2541 END DO 2542 2543 MN = 0 2544 N = 0 2545 M = 0 2546 WRITE(LUPRI,'(/1X,A10,29X,2I5,1X,G16.8)') 2547 & 'beta^K ',N,M, -AVEDISP(ITRI(M+N,M),4) 2548 DO MN = 2, NQRDSPE, 2 2549 DO M = 0, MN 2550 N = MN - M 2551 WRITE(LUPRI,'(1X,3(8X,5X),2I5,1X,G16.8)') 2552 & N,M, -AVEDISP(ITRI(M+N,M),4) 2553 END DO 2554 END DO 2555 2556 MN = 0 2557 N = 0 2558 M = 0 2559 WRITE(LUPRI,'(/1X,A10,29X,2I5,1X,G16.8)') 2560 & 'beta_{_|_}',N,M, -AVEDISP(ITRI(M+N,M),2) 2561 DO MN = 2, NQRDSPE, 2 2562 DO M = 0, MN 2563 N = MN - M 2564 WRITE(LUPRI,'(1X,3(8X,5X),2I5,1X,G16.8)') 2565 & N,M, -AVEDISP(ITRI(M+N,M),2) 2566 END DO 2567 END DO 2568 2569 MN = 0 2570 N = 0 2571 M = 0 2572 WRITE(LUPRI,'(/1X,A10,29X,2I5,1X,G16.8)') 2573 & 'beta_ms ',N,M, -AVEDISP(ITRI(M+N,M),3) 2574 DO MN = 2, NQRDSPE, 2 2575 DO M = 0, MN 2576 N = MN - M 2577 WRITE(LUPRI,'(1X,3(8X,5X),2I5,1X,G16.8)') 2578 & N,M, -AVEDISP(ITRI(M+N,M),3) 2579 END DO 2580 END DO 2581 END IF 2582 2583 WRITE(LUPRI,'(68("-"),//)') 2584 2585 WRITE(LUPRI,'(////1X,3(A," operator",3X),2(4X,A),/,91("-"))') 2586 & 'A','B','C','N','D^SHG D^OR D^EOPE' 2587 2588 DO IOPER = 1, NQROPER 2589 ISYMA = ISYOPR(IAQROP(IOPER)) 2590 ISYMB = ISYOPR(IBQROP(IOPER)) 2591 ISYMC = ISYOPR(ICQROP(IOPER)) 2592 2593 ISAMA = ISYMAT(IAQROP(IOPER)) 2594 ISAMB = ISYMAT(IBQROP(IOPER)) 2595 ISAMC = ISYMAT(ICQROP(IOPER)) 2596 2597 ISAPROP = ISAMA * ISAMB * ISAMC 2598 IF (ALLDSPCF) ISAPROP = 0 2599 2600 MN0 = 0 2601 MNS = 2 2602 2603 IF (ISAPROP.EQ.-1) MN0 = 1 2604 IF (ISAPROP.EQ. 0) MNS = 1 2605 2606 MN = MN0 2607 IF (MULD2H(ISYMA,ISYMB).EQ.ISYMC) THEN 2608 WRITE(LUPRI,'(/1X,3(A8,5X),I5,1X,3G16.8)') 2609 & LBLOPR(IAQROP(IOPER)), 2610 & LBLOPR(IBQROP(IOPER)), 2611 & LBLOPR(ICQROP(IOPER)),MN, 2612 & -SHGCF(MN+1,IOPER), 2613 & -ORCOF(MN+1,IOPER), 2614 & -EOPECF(MN+1,IOPER) 2615 ELSE 2616 IF (IPRQHYP.GT.0) THEN 2617 WRITE(LUPRI,'(/1X,3(A8,5X),A,3(A,8X))') 2618 & LBLOPR(IAQROP(IOPER)), 2619 & LBLOPR(IBQROP(IOPER)), 2620 & LBLOPR(ICQROP(IOPER)),' - ', 2621 & '---','---','---' 2622 END IF 2623 END IF 2624 2625 IF (MULD2H(ISYMA,ISYMB).EQ.ISYMC) THEN 2626 DO MN = MN0+MNS, NQRDSPE, MNS 2627 WRITE(LUPRI,'(1X,3(8X,5X),I5,1X,3G16.8)') MN, 2628 & -SHGCF(MN+1,IOPER),-ORCOF(MN+1,IOPER),-EOPECF(MN+1,IOPER) 2629 END DO 2630 END IF 2631 2632 END DO 2633 2634 IF (BETA_AVERAGE) THEN 2635 MN = 0 2636 WRITE(LUPRI,'(/1X,A10,29X,I5,1X,3G16.8)') 'beta_{||} ',MN, 2637 & -AVESHG(MN+1,1), -AVEOR(MN+1,1), -AVEEOPE(MN+1,1) 2638 DO MN = 2, NQRDSPE, 2 2639 WRITE(LUPRI,'(1X,39X,I5,1X,3G16.8)') MN, 2640 & -AVESHG(MN+1,1), -AVEOR(MN+1,1), -AVEEOPE(MN+1,1) 2641 END DO 2642 MN = 0 2643 WRITE(LUPRI,'(/1X,A10,29X,I5,1X,3G16.8)') 'beta^K ',MN, 2644 & -AVESHG(MN+1,4), -AVEOR(MN+1,4), -AVEEOPE(MN+1,4) 2645 DO MN = 2, NQRDSPE, 2 2646 WRITE(LUPRI,'(1X,39X,I5,1X,3G16.8)') MN, 2647 & -AVESHG(MN+1,4), -AVEOR(MN+1,4), -AVEEOPE(MN+1,4) 2648 END DO 2649 MN = 0 2650 WRITE(LUPRI,'(/1X,A10,29X,I5,1X,3G16.8)') 'beta_{_|_}',MN, 2651 & -AVESHG(MN+1,2), -AVEOR(MN+1,2), -AVEEOPE(MN+1,2) 2652 DO MN = 2, NQRDSPE, 2 2653 WRITE(LUPRI,'(1X,39X,I5,1X,3G16.8)') MN, 2654 & -AVESHG(MN+1,2), -AVEOR(MN+1,2), -AVEEOPE(MN+1,2) 2655 END DO 2656 MN = 0 2657 WRITE(LUPRI,'(/1X,A10,29X,I5,1X,3G16.8)') 'beta_ms ',MN, 2658 & -AVESHG(MN+1,3), -AVEOR(MN+1,3), -AVEEOPE(MN+1,3) 2659 DO MN = 2, NQRDSPE, 2 2660 WRITE(LUPRI,'(1X,39X,I5,1X,3G16.8)') MN, 2661 & -AVESHG(MN+1,3), -AVEOR(MN+1,3), -AVEEOPE(MN+1,3) 2662 END DO 2663 END IF 2664 2665 WRITE(LUPRI,'(91("-"),//)') 2666 2667 IF (BETA_AVERAGE) THEN 2668 WRITE(LUPRI,'(//////1X,A,/,49("-"))') 2669 & 'PROCESS INDEPENDENT COEFFICIENTS FOR beta_{||}:' 2670 WRITE(LUPRI,'(1X,A,G16.8)')'beta(0)',-AVEDISP(1,1) 2671 IF (NQRDSPE.GE.2)WRITE(LUPRI,'(1X,A,5X,G16.8)')"A ",ABCDE(KA,1) 2672 IF (NQRDSPE.GE.4)WRITE(LUPRI,'(1X,A,5X,G16.8)')"B ",ABCDE(KB,1) 2673 IF (NQRDSPE.GE.6) THEN 2674 WRITE(LUPRI,'(1X,A,5X,G16.8,5X,A,5X,G16.8)') 2675 & "C ",ABCDE(KC,1), "C'",ABCDE(KCP,1) 2676 END IF 2677 IF (NQRDSPE.GE.8) THEN 2678 WRITE(LUPRI,'(1X,A,5X,G16.8,5X,A,5X,G16.8)') 2679 & "D ",ABCDE(KD,1), "D'",ABCDE(KDP,1) 2680 END IF 2681 IF (NQRDSPE.GE.10) THEN 2682 WRITE(LUPRI,'(1X,A,5X,G16.8,5X,A,5X,G16.8)') 2683 & "E ",ABCDE(KE,1), "E'",ABCDE(KEP,1) 2684 END IF 2685 WRITE(LUPRI,'(49("-"),//)') 2686 END IF 2687 2688 IF (BETA_AVERAGE) THEN 2689 WRITE(LUPRI,'(///1X,A,/,101("-"))') 2690 & 'PROCESS INDEPENDENT COEFFICIENTS FOR beta_ms:' 2691 IF (NQRDSPE.GE.2) 2692 & WRITE(LUPRI,'(1X,A,2X,G16.8)') "A ",ABCDE(KA,2) 2693 IF (NQRDSPE.GE.4) THEN 2694 WRITE(LUPRI,'(1X,A,2X,G16.8,5X,A,2X,G16.8)') 2695 & "B ",ABCDE(KB,2), "B' ",ABCDE(KBP,2) 2696 END IF 2697 IF (NQRDSPE.GE.6) THEN 2698 WRITE(LUPRI,'(1X,A,2X,G16.8,5X,A,2X,G16.8)') 2699 & "C ",ABCDE(KC,2), "C' ",ABCDE(KCP,2) 2700 END IF 2701 IF (NQRDSPE.GE.8) THEN 2702 WRITE(LUPRI,'(1X,3(A,2X,G16.8,5X))') 2703 & "D ",ABCDE(KD,2),"D' ",ABCDE(KDP,2),"D'' ",ABCDE(KDP2,2) 2704 END IF 2705 IF (NQRDSPE.GE.10) THEN 2706 WRITE(LUPRI,'(1X,3(A,2X,G16.8,5X),A,2X,G16.8)') 2707 & "E ",ABCDE(KE,2), "E' ",ABCDE(KEP,2), 2708 & "E'' ",ABCDE(KEP2,2), "E'''",ABCDE(KEP3,2) 2709 END IF 2710 WRITE(LUPRI,'(101("-"),////)') 2711 END IF 2712 2713 2714 RETURN 2715 END 2716*---------------------------------------------------------------------* 2717* END OF SUBROUTINE CCQDISPRT * 2718*---------------------------------------------------------------------* 2719C /* Deck around120 */ 2720 SUBROUTINE AROUND120(HEAD) 2721 CHARACTER HEAD*(*) 2722#include "priunit.h" 2723 LNG = LEN(HEAD) + 2 2724 IND = (120 - LNG)/2 + 1 2725 WRITE (LUPRI, '(//,120A)') (' ',I=1,IND), '+', ('-',I=1,LNG), '+' 2726 WRITE (LUPRI, '(120A)') (' ',I=1,IND), '! ', HEAD, ' !' 2727 WRITE (LUPRI, '(120A)') (' ',I=1,IND), '+', ('-',I=1,LNG), '+' 2728 WRITE (LUPRI, '()') 2729 CALL FLSHFO(LUPRI) 2730 RETURN 2731 END 2732*---------------------------------------------------------------------* 2733