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_2HYP */ 20*=====================================================================* 21 SUBROUTINE CC_2HYP(WORK,LWORK) 22*---------------------------------------------------------------------* 23* 24* Purpose: direct calculation of frequency-dependent second 25* hyperpolarizabilities for the Couple Cluster models 26* 27* CCS, CC2, CCSD 28* 29* Written by Christof Haettig, february 1997. 30* Dispersion coefficients, march 1998, Christof Haettig 31* 32*=====================================================================* 33#if defined (IMPLICIT_NONE) 34 IMPLICIT NONE 35#else 36# include "implicit.h" 37#endif 38#include "priunit.h" 39#include "ccsdinp.h" 40#include "ccorb.h" 41#include "cccrinf.h" 42#include "ccroper.h" 43#include "ccr1rsp.h" 44#include "ccrc1rsp.h" 45#include "cccperm.h" 46#include "second.h" 47#include "cclists.h" 48#include "dummy.h" 49 50* local parameters: 51 CHARACTER*(19) MSGDBG 52 PARAMETER (MSGDBG = '[debug] CC_2HYP> ') 53 LOGICAL LOCDBG 54 PARAMETER (LOCDBG = .FALSE.) 55 56#if defined (SYS_CRAY) 57 REAL ZERO 58#else 59 DOUBLE PRECISION ZERO 60#endif 61 PARAMETER (ZERO = 0.0d0) 62 63 LOGICAL LCHI2_SAVE, LXKS3_SAVE, LADD, LERROR 64 CHARACTER*10 MODEL, MOPRPC 65 INTEGER LWORK 66 INTEGER IZT(4),IAM(4),IOP(4),IR2(6),IX2(6),IO3(4),IL2(6),IO2(6) 67 INTEGER ISYML, ISYM1, ISYM2, ISYM3, ISYMA, ISYMB, ISYMC, ISYMD 68 INTEGER KHYPPOL, NBHYPPOL, IKAP(4) 69 INTEGER ITRAN, IVEC, IDX, I, ISIGN, IFREQ, IOPER, IOPT 70 INTEGER K0HTRAN, K0HDOTS, K0HCONS, KXCONS, KXDOTS, KXTRAN 71 INTEGER K0GTRAN, K0GDOTS, K0GCONS, KAGTRAN, KAGDOTS, KAGCONS 72 INTEGER K0FTRAN, K0FDOTS, K0FCONS, KAFTRAN, KAFDOTS, KAFCONS 73 INTEGER K0FATRAN,K0FADOTS,K0FACONS,KAFATRAN,KAFADOTS,KAFACONS 74 INTEGER K0EATRAN,K0EADOTS,K0EACONS,KAEATRAN,KAEADOTS,KAEACONS 75 INTEGER KOCONS, KODOTS, KOTRAN, KLCONS, KLDOTS, KLTRAN, KEXPCOF 76 INTEGER N0HTRAN, N0GTRAN, NAGTRAN, N0FTRAN, NAFTRAN, NBEXPCOF 77 INTEGER N0FATRAN, NAFATRAN, NAEATRAN, NXTRAN, NOTRAN, NLTRAN 78 INTEGER MX0HTRAN, MX0GTRAN, MXAGTRAN, MX0FTRAN, MXAFTRAN 79 INTEGER MX0FATRAN, MXAFATRAN, MXAEATRAN, MXXTRAN, MXOTRAN,MXLTRAN 80 INTEGER MX0HDOTS, MX0GDOTS, MXAGDOTS, MX0FDOTS, MXAFDOTS 81 INTEGER MX0FADOTS, MXAFADOTS, MXAEADOTS, MXXDOTS, MXODOTS,MXLDOTS 82 INTEGER MXTRAN3, MXTRAN2, MXVEC1, MXVEC2, IORDER 83 INTEGER KEND0, LEND0 84 INTEGER P 85 INTEGER NTRINOM, KTRINOM, KEND0A, KEND1, LEND1 86 INTEGER KDSPCF, KTHGCF, KESHGCF, KKERRCF, KDFWMCF 87 INTEGER KDSPCFA, KTHGCFA, KESHGCA, KKERRCA, KDFWMCA, KABCDE 88 89 90#if defined (SYS_CRAY) 91 REAL WORK(LWORK) 92 REAL TIM0, TIM1, TIMH, TIMG, TIMF, TIMFA, TIMEA, TIMX2, TIMO2 93 REAL TIM2, TIML2 94 REAL H0CON, G0CON(6), GACON(4), F0CON(3), FACON(12) 95 REAL F0ACON(12), FAACON(12), EAACON(12), XCON(6) 96 REAL OCON(4), LCON(3) 97 REAL SIGN 98#else 99 DOUBLE PRECISION WORK(LWORK) 100 DOUBLE PRECISION TIM0,TIM1,TIMH,TIMG,TIMF,TIMFA,TIMEA,TIMX2,TIMO2 101 DOUBLE PRECISION TIM2,TIML2 102 DOUBLE PRECISION H0CON, G0CON(6), GACON(4), F0CON(3), FACON(12) 103 DOUBLE PRECISION F0ACON(12), FAACON(12), EAACON(12), XCON(6) 104 DOUBLE PRECISION OCON(4), LCON(3) 105 DOUBLE PRECISION SIGN 106#endif 107 108* external functions 109 INTEGER IR3TAMP 110 INTEGER IR2TAMP 111 INTEGER IR1TAMP 112 INTEGER IL1ZETA 113 INTEGER IL2ZETA 114 INTEGER ICHI2 115 INTEGER IRHSR2 116 INTEGER IRHSR3 117 118 119*---------------------------------------------------------------------* 120* print header for hyperpolarizability section 121*---------------------------------------------------------------------* 122 WRITE (LUPRI,'(7(/1X,2A),/)') 123 & '************************************', 124 & '*******************************', 125 & '* ', 126 & ' *', 127 & '*---------- OUTPUT FROM COUPLED CLUST', 128 & 'ER CUBIC RESPONSE ---------*', 129 & '* ', 130 & ' *', 131 & '*-------- CALCULATION OF CC HYPER', 132 & 'POLARIZABILITIES ---------*', 133 & '* ', 134 & ' *', 135 & '************************************', 136 & '*******************************' 137 138*---------------------------------------------------------------------* 139 IF (.NOT. (CCS .OR. CC2 .OR. CCSD .OR. CC3) ) THEN 140 CALL QUIT('CC_2HYP called for unknown Coupled Cluster.') 141 END IF 142 143* print some debug/info output 144 IF (IPRCHYP .GT. 10) WRITE(LUPRI,*) 'CC_2HYP Workspace:',LWORK 145 146 TIM0 = SECOND() 147 148 IF (NCRFREQ.EQ.0 .AND. NCRDISP.EQ.0) THEN 149 WRITE(LUPRI,*) 'Nothing to be done in this section...' 150 RETURN 151 END IF 152 153 IF (NCRFREQ.GT.0) THEN 154*---------------------------------------------------------------------* 155* allocate & initialize work space for hyperpolarizabilities 156*---------------------------------------------------------------------* 157 NBHYPPOL = NCROPER * NCRFREQ 158 159 KHYPPOL = 1 160 KEND0 = KHYPPOL + 2*NBHYPPOL 161 162 163 MXTRAN2 = 2 * NLRTLBL * NLRTLBL 164 MXVEC2 = (NLRTLBL+1) * NLRTLBL /2 165 MXTRAN3 = 2 * NLRTLBL * NLRTLBL * NLRTLBL 166 MXVEC1 = NLRTLBL 167 MXTRAN2 = MIN(MXTRAN2,2*12*NBHYPPOL) ! 12 = max. # of permut. 168 MXTRAN3 = MIN(MXTRAN3,2*12*NBHYPPOL) ! 12 = max. # of permut. 169 MXVEC2 = MIN(MXVEC2,2*12*NBHYPPOL) ! 12 = max. # of permut. 170 MXVEC1 = MIN(MXVEC1,2*12*NBHYPPOL) ! 12 = max. # of permut. 171 172 MX0HTRAN = MXDIM_HTRAN * MXTRAN3 173 MX0GTRAN = MXDIM_GTRAN * MXTRAN2 174 MXAGTRAN = MXDIM_GTRAN * MXTRAN3 175 MX0FTRAN = MXDIM_FTRAN * MXTRAN2 176 MXAFTRAN = MXDIM_FTRAN * MXTRAN2 177 MX0FATRAN = MXDIM_FATRAN * MXTRAN2 178 MXAFATRAN = MXDIM_FATRAN * MXTRAN3 179 MXAEATRAN = MXDIM_XEVEC * MXTRAN2 180 MXXTRAN = 1 * MXTRAN2 181 MXOTRAN = 1 * MXTRAN3 182 MXLTRAN = 1 * MXTRAN2 183 184 MX0HDOTS = MXTRAN3 * MXVEC1 185 MX0GDOTS = MXTRAN2 * MXVEC2 186 MXAGDOTS = MXTRAN3 * MXVEC1 187 MX0FDOTS = MXTRAN2 * MXVEC2 188 MXAFDOTS = MXTRAN2 * MXVEC2 189 MX0FADOTS = MXTRAN2 * MXVEC2 190 MXAFADOTS = MXTRAN3 * MXVEC1 191 MXAEADOTS = MXTRAN2 * MXVEC2 192 MXXDOTS = MXTRAN2 * MXVEC2 193 MXODOTS = MXTRAN3 * MXVEC1 194 MXLDOTS = MXTRAN2 * MXVEC2 195 196 K0HTRAN = KEND0 197 K0HDOTS = K0HTRAN + MX0HTRAN 198 K0HCONS = K0HDOTS + MX0HDOTS 199 KEND0 = K0HCONS + MX0HDOTS 200 201 K0GTRAN = KEND0 202 K0GDOTS = K0GTRAN + MX0GTRAN 203 K0GCONS = K0GDOTS + MX0GDOTS 204 KEND0 = K0GCONS + MX0GDOTS 205 206 KAGTRAN = KEND0 207 KAGDOTS = KAGTRAN + MXAGTRAN 208 KAGCONS = KAGDOTS + MXAGDOTS 209 KEND0 = KAGCONS + MXAGDOTS 210 211 K0FTRAN = KEND0 212 K0FDOTS = K0FTRAN + MX0FTRAN 213 K0FCONS = K0FDOTS + MX0FDOTS 214 KEND0 = K0FCONS + MX0FDOTS 215 216 KAFTRAN = KEND0 217 KAFDOTS = KAFTRAN + MXAFTRAN 218 KAFCONS = KAFDOTS + MXAFDOTS 219 KEND0 = KAFCONS + MXAFDOTS 220 221 K0FATRAN = KEND0 222 K0FADOTS = K0FATRAN + MX0FATRAN 223 K0FACONS = K0FADOTS + MX0FADOTS 224 KEND0 = K0FACONS + MX0FADOTS 225 226 KAFATRAN = KEND0 227 KAFADOTS = KAFATRAN + MXAFATRAN 228 KAFACONS = KAFADOTS + MXAFADOTS 229 KEND0 = KAFACONS + MXAFADOTS 230 231 KAEATRAN = KEND0 232 KAEADOTS = KAEATRAN + MXAEATRAN 233 KAEACONS = KAEADOTS + MXAEADOTS 234 KEND0 = KAEACONS + MXAEADOTS 235 236 KXTRAN = KEND0 237 KXDOTS = KXTRAN + MXXTRAN 238 KXCONS = KXDOTS + MXXDOTS 239 KEND0 = KXCONS + MXXDOTS 240 241 KOTRAN = KEND0 242 KODOTS = KOTRAN + MXOTRAN 243 KOCONS = KODOTS + MXODOTS 244 KEND0 = KOCONS + MXODOTS 245 246 KLTRAN = KEND0 247 KLDOTS = KLTRAN + MXLTRAN 248 KLCONS = KLDOTS + MXLDOTS 249 KEND0 = KLCONS + MXLDOTS 250 251 LEND0 = LWORK - KEND0 252 253 IF (LEND0 .LT. 0) THEN 254 CALL QUIT('Insufficient memory in CC_2HYP. (1)') 255 END IF 256 257* initialize hyperpolarizabilities and the lists of contributions: 258 CALL DZERO(WORK(KHYPPOL),2*NBHYPPOL) 259 CALL DZERO(WORK(K0HTRAN),MX0HTRAN+2*MX0HDOTS) 260 CALL DZERO(WORK(K0GTRAN),MX0GTRAN+2*MX0GDOTS) 261 CALL DZERO(WORK(KAGTRAN),MXAGTRAN+2*MXAGDOTS) 262 CALL DZERO(WORK(K0FTRAN),MX0FTRAN+2*MX0FDOTS) 263 CALL DZERO(WORK(KAFTRAN),MXAFTRAN+2*MXAFDOTS) 264 CALL DZERO(WORK(K0FATRAN),MX0FATRAN+2*MX0FADOTS) 265 CALL DZERO(WORK(KAFATRAN),MXAFATRAN+2*MXAFADOTS) 266 CALL DZERO(WORK(KAEATRAN),MXAEATRAN+2*MXAEADOTS) 267 CALL DZERO(WORK(KXTRAN),MXXTRAN+2*MXXDOTS) 268 CALL DZERO(WORK(KOTRAN),MXOTRAN+2*MXODOTS) 269 CALL DZERO(WORK(KLTRAN),MXLTRAN+2*MXLDOTS) 270 271*---------------------------------------------------------------------* 272* set up lists for H, G, F and F{O} transformations and ETA{O} vectors: 273*---------------------------------------------------------------------* 274 CALL CCCR_SETUP(MXTRAN2, MXVEC2, MXTRAN3, MXVEC1, 275 & WORK(K0HTRAN), WORK(K0HDOTS), N0HTRAN, 276 & WORK(K0GTRAN), WORK(K0GDOTS), N0GTRAN, 277 & WORK(KAGTRAN), WORK(KAGDOTS), NAGTRAN, 278 & WORK(K0FTRAN), WORK(K0FDOTS), N0FTRAN, 279 & WORK(KAFTRAN), WORK(KAFDOTS), NAFTRAN, 280 & WORK(K0FATRAN),WORK(K0FADOTS),N0FATRAN, 281 & WORK(KAFATRAN),WORK(KAFADOTS),NAFATRAN, 282 & WORK(KAEATRAN),WORK(KAEADOTS),NAEATRAN, 283 & WORK(KXTRAN), WORK(KXDOTS), NXTRAN, 284 & WORK(KOTRAN), WORK(KODOTS), NOTRAN, 285 & WORK(KLTRAN), WORK(KLDOTS), NLTRAN ) 286 287*---------------------------------------------------------------------* 288* calculate H matrix contributions: 289*---------------------------------------------------------------------* 290 TIM1 = SECOND() 291 292 IOPT = 5 293 CALL CC_HMAT('L0 ','R1 ','R1 ','R1 ','R1 ',N0HTRAN, MXVEC1, 294 & WORK(K0HTRAN),WORK(K0HDOTS),WORK(K0HCONS), 295 & WORK(KEND0), LEND0, IOPT ) 296 297 TIMH = SECOND() - TIM1 298 WRITE (LUPRI,'(/A,I5,A,F12.2," seconds.")') 299 & ' Time used for',N0HTRAN,' H matrix transformations: ',TIMH 300 CALL FLSHFO(LUPRI) 301 302*---------------------------------------------------------------------* 303* calculate G matrix contributions: 304*---------------------------------------------------------------------* 305 TIM1 = SECOND() 306 307 IF (.NOT.L_USE_CHI2) THEN 308 IOPT = 5 309 CALL CC_GMATRIX('L0 ','R1 ','R1 ','R2 ',N0GTRAN, MXVEC2, 310 & WORK(K0GTRAN),WORK(K0GDOTS),WORK(K0GCONS), 311 & WORK(KEND0), LEND0, IOPT ) 312 END IF 313 314 IF (.NOT.L_USE_XKS3) THEN 315 IOPT = 5 316 CALL CC_GMATRIX('L1 ','R1 ','R1 ','R1 ',NAGTRAN, MXVEC1, 317 & WORK(KAGTRAN),WORK(KAGDOTS),WORK(KAGCONS), 318 & WORK(KEND0), LEND0, IOPT ) 319 END IF 320 321 TIMG = SECOND() - TIM1 322 WRITE (LUPRI,'(/A,I5,A,F12.2," seconds.")') 323 & ' Time used for',N0GTRAN+NAGTRAN, 324 & ' G matrix transformations: ',TIMG 325 CALL FLSHFO(LUPRI) 326 327*---------------------------------------------------------------------* 328* calculate F matrix contributions: 329*---------------------------------------------------------------------* 330 TIM1 = SECOND() 331 332 IOPT = 5 333 CALL CC_FMATRIX(WORK(K0FTRAN),N0FTRAN,'L0 ','R2 ',IOPT,'R2 ', 334 & WORK(K0FDOTS),WORK(K0FCONS),MXVEC2, 335 & WORK(KEND0),LEND0) 336 IF (L_USE_CHI2) THEN 337 CALL DSCAL(MX0FDOTS,-1.0d0,WORK(K0FCONS),1) 338 END IF 339 340 IF (.NOT. (L_USE_CHI2 .OR. L_USE_XKS3) ) THEN 341 IOPT = 5 342 CALL CC_FMATRIX(WORK(KAFTRAN),NAFTRAN,'L1 ','R1 ',IOPT,'R2 ', 343 & WORK(KAFDOTS),WORK(KAFCONS),MXVEC2, 344 & WORK(KEND0),LEND0) 345 END IF 346 347 TIMF = SECOND() - TIM1 348 WRITE (LUPRI,'(/A,I5,A,F12.2," seconds.")') 349 & ' Time used for',N0FTRAN+NAFTRAN, 350 & ' F matrix transformations: ',TIMF 351 CALL FLSHFO(LUPRI) 352 353*---------------------------------------------------------------------* 354* calculate F{O} matrix contributions: 355*---------------------------------------------------------------------* 356 TIM1 = SECOND() 357 358 IF (.NOT.L_USE_CHI2) THEN 359 CALL CCQR_FADRV('L0 ','o1 ','R1 ','R2 ',N0FATRAN, MXVEC2, 360 & WORK(K0FATRAN),WORK(K0FADOTS), 361 & WORK(K0FACONS), 362 & WORK(KEND0), LEND0, 'DOTP' ) 363 END IF 364 365 IF (.NOT.L_USE_XKS3) THEN 366 CALL CCQR_FADRV('L1 ','o1 ','R1 ','R1 ',NAFATRAN, MXVEC1, 367 & WORK(KAFATRAN),WORK(KAFADOTS), 368 & WORK(KAFACONS), 369 & WORK(KEND0), LEND0, 'DOTP' ) 370 END IF 371 372 TIMFA = SECOND() - TIM1 373 WRITE (LUPRI,'(/A,I5,A,F12.2," seconds.")') 374 & ' Time used for',N0FATRAN+NAFATRAN, 375 & ' F{O} matrix transformations:',TIMFA 376 CALL FLSHFO(LUPRI) 377 378*---------------------------------------------------------------------* 379* calculate ETA{O} vector contributions: 380*---------------------------------------------------------------------* 381 TIM1 = SECOND() 382 383 IF (.NOT. (L_USE_CHI2 .OR. L_USE_XKS3) ) THEN 384 385 IOPT = 5 386 IORDER = 1 387 CALL CC_XIETA( WORK(KAEATRAN), NAEATRAN, IOPT, IORDER, 'L1 ', 388 & '---',IDUMMY, DUMMY, 389 & 'R2 ',WORK(KAEADOTS),WORK(KAEACONS), 390 & .FALSE.,MXVEC2, WORK(KEND0), LEND0 ) 391 392 END IF 393 394 TIMEA = SECOND() - TIM1 395 WRITE (LUPRI,'(/A,I5,A,F12.2," seconds.")') ' Time used for', 396 & NAEATRAN,' ETA{O} vector calculations :',TIMEA 397 CALL FLSHFO(LUPRI) 398 399*---------------------------------------------------------------------* 400* chi vector contributions: 401*---------------------------------------------------------------------* 402 IF (L_USE_CHI2) THEN 403 TIM1 = SECOND() 404 405 CALL CC_DOTDRV('X2B','R2 ',NXTRAN,MXVEC2, 406 & WORK(KXTRAN), WORK(KXDOTS), WORK(KXCONS), 407 & WORK(KEND0), LEND0 ) 408 409 TIMX2 = SECOND() - TIM1 410 WRITE (LUPRI,'(/A,I5,A,F12.2," seconds.")') 411 & ' Time used for',NXTRAN, ' X2 x R2 dot products:',TIMX2 412 CALL FLSHFO(LUPRI) 413 414 END IF 415 416*---------------------------------------------------------------------* 417* L2 x O2 dot products: 418*---------------------------------------------------------------------* 419 IF (USE_L2BC .OR. USE_LBCD) THEN 420 TIM1 = SECOND() 421 422 CALL CC_DOTDRV('L2 ','O2 ',NLTRAN,MXVEC2, 423 & WORK(KLTRAN), WORK(KLDOTS), WORK(KLCONS), 424 & WORK(KEND0), LEND0 ) 425 426 TIML2 = SECOND() - TIM1 427 WRITE (LUPRI,'(/A,I5,A,5X,F12.2," seconds.")') 428 & ' Time used for',NLTRAN, ' L2 x R2 dot products: ',TIML2 429 CALL FLSHFO(LUPRI) 430 431 END IF 432*---------------------------------------------------------------------* 433* xksi vector contributions: 434*---------------------------------------------------------------------* 435 IF (L_USE_XKS3) THEN 436 437 TIM1 = SECOND() 438 439 CALL CC_DOTDRV('R3 ','X1B',NOTRAN,MXVEC1, 440 & WORK(KOTRAN), WORK(KODOTS), WORK(KOCONS), 441 & WORK(KEND0), LEND0 ) 442 443 TIMO2 = SECOND() - TIM1 444 WRITE (LUPRI,'(/A,I5,A,F12.2," seconds.")') 445 & ' Time used for',NOTRAN, ' R3 x X1 dot products:',TIMO2 446 CALL FLSHFO(LUPRI) 447 448 END IF 449*=====================================================================* 450* loop over quadruples of operator labels and over frequencies and 451* collect all contributions to the final hyperpolarizabilities: 452*=====================================================================* 453 DO IOPER = 1, NCROPER 454 455 IOP(A) = IACROP(IOPER) 456 IOP(B) = IBCROP(IOPER) 457 IOP(C) = ICCROP(IOPER) 458 IOP(D) = IDCROP(IOPER) 459 460 IKAP(A)= 0 461 IKAP(B)= 0 462 IKAP(C)= 0 463 IKAP(D)= 0 464 465 ISYMB = ISYOPR(IOP(B)) 466 ISYMC = ISYOPR(IOP(C)) 467 ISYMA = ISYOPR(IOP(A)) 468 ISYMD = ISYOPR(IOP(D)) 469 470 471 IF (MULD2H(ISYMA,ISYMB) .EQ. MULD2H(ISYMC,ISYMD)) THEN 472 473 DO IFREQ = 1, NCRFREQ 474 475 DO ISIGN = +1, -1, -2 476 SIGN = DBLE(ISIGN) 477 478 IZT(A) =IL1ZETA(LBLOPR(IOP(A)),.FALSE.,SIGN*ACRFR(IFREQ),ISYML) 479 IZT(B) =IL1ZETA(LBLOPR(IOP(B)),.FALSE.,SIGN*BCRFR(IFREQ),ISYML) 480 IZT(C) =IL1ZETA(LBLOPR(IOP(C)),.FALSE.,SIGN*CCRFR(IFREQ),ISYML) 481 IZT(D) =IL1ZETA(LBLOPR(IOP(D)),.FALSE.,SIGN*DCRFR(IFREQ),ISYML) 482 483 IAM(A) =IR1TAMP(LBLOPR(IOP(A)),.FALSE.,SIGN*ACRFR(IFREQ),ISYML) 484 IAM(B) =IR1TAMP(LBLOPR(IOP(B)),.FALSE.,SIGN*BCRFR(IFREQ),ISYML) 485 IAM(C) =IR1TAMP(LBLOPR(IOP(C)),.FALSE.,SIGN*CCRFR(IFREQ),ISYML) 486 IAM(D) =IR1TAMP(LBLOPR(IOP(D)),.FALSE.,SIGN*DCRFR(IFREQ),ISYML) 487 488 IF ( USE_LBCD ) THEN 489 IL2(BC)=IL2ZETA(LBLOPR(IOP(B)),SIGN*BCRFR(IFREQ),ISYM1, 490 & LBLOPR(IOP(C)),SIGN*CCRFR(IFREQ),ISYM2) 491 IL2(BD)=IL2ZETA(LBLOPR(IOP(B)),SIGN*BCRFR(IFREQ),ISYM1, 492 & LBLOPR(IOP(D)),SIGN*DCRFR(IFREQ),ISYM2) 493 IL2(CD)=IL2ZETA(LBLOPR(IOP(C)),SIGN*CCRFR(IFREQ),ISYM1, 494 & LBLOPR(IOP(D)),SIGN*DCRFR(IFREQ),ISYM2) 495 IO2(AB)=IRHSR2( LBLOPR(IOP(A)),.FALSE.,SIGN*ACRFR(IFREQ),ISYM1, 496 & LBLOPR(IOP(B)),.FALSE.,SIGN*BCRFR(IFREQ),ISYM2) 497 IO2(AC)=IRHSR2( LBLOPR(IOP(A)),.FALSE.,SIGN*ACRFR(IFREQ),ISYM1, 498 & LBLOPR(IOP(C)),.FALSE.,SIGN*CCRFR(IFREQ),ISYM2) 499 IO2(AD)=IRHSR2( LBLOPR(IOP(A)),.FALSE.,SIGN*ACRFR(IFREQ),ISYM1, 500 & LBLOPR(IOP(D)),.FALSE.,SIGN*DCRFR(IFREQ),ISYM2) 501 ELSE IF ( USE_L2BC) THEN 502 IR2(AB)=IR2TAMP(LBLOPR(IOP(A)),.FALSE.,SIGN*ACRFR(IFREQ),ISYM1, 503 & LBLOPR(IOP(B)),.FALSE.,SIGN*BCRFR(IFREQ),ISYM2) 504 IR2(AC)=IR2TAMP(LBLOPR(IOP(A)),.FALSE.,SIGN*ACRFR(IFREQ),ISYM1, 505 & LBLOPR(IOP(C)),.FALSE.,SIGN*CCRFR(IFREQ),ISYM2) 506 IL2(BC)=IL2ZETA(LBLOPR(IOP(B)),SIGN*BCRFR(IFREQ),ISYM1, 507 & LBLOPR(IOP(C)),SIGN*CCRFR(IFREQ),ISYM2) 508 IO2(AD)=IRHSR2( LBLOPR(IOP(A)),.FALSE.,SIGN*ACRFR(IFREQ),ISYM1, 509 & LBLOPR(IOP(D)),.FALSE.,SIGN*DCRFR(IFREQ),ISYM2) 510 ELSE 511 IR2(AB)=IR2TAMP(LBLOPR(IOP(A)),.FALSE.,SIGN*ACRFR(IFREQ),ISYM1, 512 & LBLOPR(IOP(B)),.FALSE.,SIGN*BCRFR(IFREQ),ISYM2) 513 IR2(AC)=IR2TAMP(LBLOPR(IOP(A)),.FALSE.,SIGN*ACRFR(IFREQ),ISYM1, 514 & LBLOPR(IOP(C)),.FALSE.,SIGN*CCRFR(IFREQ),ISYM2) 515 IR2(AD)=IR2TAMP(LBLOPR(IOP(A)),.FALSE.,SIGN*ACRFR(IFREQ),ISYM1, 516 & LBLOPR(IOP(D)),.FALSE.,SIGN*DCRFR(IFREQ),ISYM2) 517 END IF 518 519 IR2(BC)=IR2TAMP(LBLOPR(IOP(B)),.FALSE.,SIGN*BCRFR(IFREQ),ISYM1, 520 & LBLOPR(IOP(C)),.FALSE.,SIGN*CCRFR(IFREQ),ISYM2) 521 IR2(BD)=IR2TAMP(LBLOPR(IOP(B)),.FALSE.,SIGN*BCRFR(IFREQ),ISYM1, 522 & LBLOPR(IOP(D)),.FALSE.,SIGN*DCRFR(IFREQ),ISYM2) 523 IR2(CD)=IR2TAMP(LBLOPR(IOP(C)),.FALSE.,SIGN*CCRFR(IFREQ),ISYM1, 524 & LBLOPR(IOP(D)),.FALSE.,SIGN*DCRFR(IFREQ),ISYM2) 525 526 IF (L_USE_CHI2) THEN 527c IX2(AB)=ICHI2(LBLOPR(IOP(A)),.FALSE.,SIGN*ACRFR(IFREQ),ISYM1, 528c & LBLOPR(IOP(B)),.FALSE.,SIGN*BCRFR(IFREQ),ISYM2) 529c IX2(AC)=ICHI2(LBLOPR(IOP(A)),.FALSE.,SIGN*ACRFR(IFREQ),ISYM1, 530c & LBLOPR(IOP(C)),.FALSE.,SIGN*CCRFR(IFREQ),ISYM2) 531c IX2(AD)=ICHI2(LBLOPR(IOP(A)),.FALSE.,SIGN*ACRFR(IFREQ),ISYM1, 532c & LBLOPR(IOP(D)),.FALSE.,SIGN*DCRFR(IFREQ),ISYM2) 533c IX2(BC)=ICHI2(LBLOPR(IOP(B)),.FALSE.,SIGN*BCRFR(IFREQ),ISYM1, 534c & LBLOPR(IOP(C)),.FALSE.,SIGN*CCRFR(IFREQ),ISYM2) 535c IX2(BD)=ICHI2(LBLOPR(IOP(B)),.FALSE.,SIGN*BCRFR(IFREQ),ISYM1, 536c & LBLOPR(IOP(D)),.FALSE.,SIGN*DCRFR(IFREQ),ISYM2) 537c IX2(CD)=ICHI2(LBLOPR(IOP(C)),.FALSE.,SIGN*CCRFR(IFREQ),ISYM1, 538c & LBLOPR(IOP(D)),.FALSE.,SIGN*DCRFR(IFREQ),ISYM2) 539 540 IX2(AB)=IL2ZETA(LBLOPR(IOP(A)), SIGN*ACRFR(IFREQ),ISYM1, 541 & LBLOPR(IOP(B)), SIGN*BCRFR(IFREQ),ISYM2) 542 IX2(AC)=IL2ZETA(LBLOPR(IOP(A)), SIGN*ACRFR(IFREQ),ISYM1, 543 & LBLOPR(IOP(C)), SIGN*CCRFR(IFREQ),ISYM2) 544 IX2(AD)=IL2ZETA(LBLOPR(IOP(A)), SIGN*ACRFR(IFREQ),ISYM1, 545 & LBLOPR(IOP(D)), SIGN*DCRFR(IFREQ),ISYM2) 546 IX2(BC)=IL2ZETA(LBLOPR(IOP(B)), SIGN*BCRFR(IFREQ),ISYM1, 547 & LBLOPR(IOP(C)), SIGN*CCRFR(IFREQ),ISYM2) 548 IX2(BD)=IL2ZETA(LBLOPR(IOP(B)), SIGN*BCRFR(IFREQ),ISYM1, 549 & LBLOPR(IOP(D)), SIGN*DCRFR(IFREQ),ISYM2) 550 IX2(CD)=IL2ZETA(LBLOPR(IOP(C)), SIGN*CCRFR(IFREQ),ISYM1, 551 & LBLOPR(IOP(D)), SIGN*DCRFR(IFREQ),ISYM2) 552 END IF 553 554 IF (L_USE_XKS3) THEN 555 IO3(ABC) = IR3TAMP(LBLOPR(IOP(A)),SIGN*ACRFR(IFREQ),ISYM1, 556 & LBLOPR(IOP(B)),SIGN*BCRFR(IFREQ),ISYM2, 557 & LBLOPR(IOP(C)),SIGN*CCRFR(IFREQ),ISYM3) 558 IO3(ABD) = IR3TAMP(LBLOPR(IOP(A)),SIGN*ACRFR(IFREQ),ISYM1, 559 & LBLOPR(IOP(B)),SIGN*BCRFR(IFREQ),ISYM2, 560 & LBLOPR(IOP(D)),SIGN*DCRFR(IFREQ),ISYM3) 561 IO3(ACD) = IR3TAMP(LBLOPR(IOP(A)),SIGN*ACRFR(IFREQ),ISYM1, 562 & LBLOPR(IOP(C)),SIGN*CCRFR(IFREQ),ISYM2, 563 & LBLOPR(IOP(D)),SIGN*DCRFR(IFREQ),ISYM3) 564 IO3(BCD) = IR3TAMP(LBLOPR(IOP(B)),SIGN*BCRFR(IFREQ),ISYM1, 565 & LBLOPR(IOP(C)),SIGN*CCRFR(IFREQ),ISYM2, 566 & LBLOPR(IOP(D)),SIGN*DCRFR(IFREQ),ISYM3) 567 END IF 568 569 570*---------------------------------------------------------------------* 571* get H^0 matrix transformations, 1 permutation 572*---------------------------------------------------------------------* 573 CALL CC_SETH1111(WORK(K0HTRAN),WORK(K0HDOTS),MXTRAN3,MXVEC1, 574 & 0,IAM(A),IAM(B),IAM(C),IAM(D),ITRAN,IVEC) 575 H0CON = WORK(K0HCONS-1 + (ITRAN-1)*MXVEC1 + IVEC) 576 577*---------------------------------------------------------------------* 578* get G^0 matrix transformations, 6 permutations 579*---------------------------------------------------------------------* 580 DO P = 1, 6 581 IF ((USE_L2BC.AND.P.EQ.4) .OR. (USE_LBCD.AND.P.GT.3) ) THEN 582 G0CON(P) = ZERO 583 ELSE 584 CALL CC_SETG112(WORK(K0GTRAN),WORK(K0GDOTS),MXTRAN2,MXVEC2, 585 & 0,IAM(I1(P)),IAM(I2(P)),IR2(IP(P)),ITRAN,IVEC) 586 G0CON(P) = WORK(K0GCONS-1 + (ITRAN-1)*MXVEC2 + IVEC) 587 END IF 588 END DO 589 590*---------------------------------------------------------------------* 591* get G^A matrix transformations, 4 permutations 592*---------------------------------------------------------------------* 593 CALL CCQR_SETG(WORK(KAGTRAN),WORK(KAGDOTS),MXTRAN3,MXVEC1, 594 & IZT(A),IAM(B),IAM(C),IAM(D),ITRAN,IVEC) 595 GACON(1) = WORK(KAGCONS-1 + (ITRAN-1)*MXVEC1 + IVEC) 596 597 CALL CCQR_SETG(WORK(KAGTRAN),WORK(KAGDOTS),MXTRAN3,MXVEC1, 598 & IZT(B),IAM(A),IAM(C),IAM(D),ITRAN,IVEC) 599 GACON(2) = WORK(KAGCONS-1 + (ITRAN-1)*MXVEC1 + IVEC) 600 601 CALL CCQR_SETG(WORK(KAGTRAN),WORK(KAGDOTS),MXTRAN3,MXVEC1, 602 & IZT(C),IAM(B),IAM(A),IAM(D),ITRAN,IVEC) 603 GACON(3) = WORK(KAGCONS-1 + (ITRAN-1)*MXVEC1 + IVEC) 604 605 CALL CCQR_SETG(WORK(KAGTRAN),WORK(KAGDOTS),MXTRAN3,MXVEC1, 606 & IZT(D),IAM(B),IAM(C),IAM(A),ITRAN,IVEC) 607 GACON(4) = WORK(KAGCONS-1 + (ITRAN-1)*MXVEC1 + IVEC) 608 609*---------------------------------------------------------------------* 610* get F^0 matrix transformations, 3 permutations 611*---------------------------------------------------------------------* 612 IF (.NOT. USE_LBCD) THEN 613 CALL CCQR_SETF(WORK(K0FTRAN),WORK(K0FDOTS),MXTRAN2,MXVEC2, 614 & 0,IR2(AB),IR2(CD),ITRAN,IVEC) 615 F0CON(1) = WORK(K0FCONS-1 + (ITRAN-1)*MXVEC2 + IVEC) 616 ELSE 617 F0CON(1) = ZERO 618 END IF 619 620 IF (.NOT. USE_LBCD) THEN 621 CALL CCQR_SETF(WORK(K0FTRAN),WORK(K0FDOTS),MXTRAN2,MXVEC2, 622 & 0,IR2(AC),IR2(BD),ITRAN,IVEC) 623 F0CON(2) = WORK(K0FCONS-1 + (ITRAN-1)*MXVEC2 + IVEC) 624 ELSE 625 F0CON(2) = ZERO 626 END IF 627 628 IF (.NOT. (USE_LBCD .OR. USE_L2BC)) THEN 629 CALL CCQR_SETF(WORK(K0FTRAN),WORK(K0FDOTS),MXTRAN2,MXVEC2, 630 & 0,IR2(AD),IR2(BC),ITRAN,IVEC) 631 F0CON(3) = WORK(K0FCONS-1 + (ITRAN-1)*MXVEC2 + IVEC) 632 ELSE 633 F0CON(3) = ZERO 634 END IF 635 636*---------------------------------------------------------------------* 637* get F^A matrix transformations, 12 permutations 638*---------------------------------------------------------------------* 639 DO P = 1, 6 640 IF ((USE_L2BC.AND.P.EQ.4) .OR. (USE_LBCD.AND.P.GT.3) ) THEN 641 FACON(P) = ZERO 642 FACON(6+P) = ZERO 643 ELSE 644 CALL CC_SETF12(WORK(KAFTRAN),WORK(KAFDOTS),MXTRAN2,MXVEC2, 645 & IZT(I1(P)),IAM(I2(P)),IR2(IP(P)),ITRAN,IVEC) 646 FACON(P) = WORK(KAFCONS-1 + (ITRAN-1)*MXVEC2 + IVEC) 647 648 CALL CC_SETF12(WORK(KAFTRAN),WORK(KAFDOTS),MXTRAN2,MXVEC2, 649 & IZT(I2(P)),IAM(I1(P)),IR2(IP(P)),ITRAN,IVEC) 650 FACON(6+P) = WORK(KAFCONS-1 + (ITRAN-1)*MXVEC2 + IVEC) 651 END IF 652 END DO 653 654*---------------------------------------------------------------------* 655* get F^0{O} matrix transformations, 12 permutations 656*---------------------------------------------------------------------* 657 DO P = 1, 6 658 IF ((USE_L2BC.AND.P.EQ.4) .OR. (USE_LBCD.AND.P.GT.3) ) THEN 659 F0ACON(P) = ZERO 660 F0ACON(6+P) = ZERO 661 ELSE 662 CALL CC_SETFA12(WORK(K0FATRAN),WORK(K0FADOTS),MXTRAN2,MXVEC2, 663 & 0,IOP(I1(P)),IAM(I2(P)),IR2(IP(P)),ITRAN,IVEC) 664 F0ACON(P) = WORK(K0FACONS-1 + (ITRAN-1)*MXVEC2 + IVEC) 665 666 CALL CC_SETFA12(WORK(K0FATRAN),WORK(K0FADOTS),MXTRAN2,MXVEC2, 667 & 0,IOP(I2(P)),IAM(I1(P)),IR2(IP(P)),ITRAN,IVEC) 668 F0ACON(6+P) = WORK(K0FACONS-1 + (ITRAN-1)*MXVEC2 + IVEC) 669 END IF 670 END DO 671 672*---------------------------------------------------------------------* 673* get F^A{O} matrix transformations, 12 permutations 674*---------------------------------------------------------------------* 675 DO P = 1, 6 676 CALL CCQR_SETFA(WORK(KAFATRAN),WORK(KAFADOTS),MXTRAN3,MXVEC1, 677 & IZT(I1(P)),IOP(I2(P)),IAM(I3(P)),IAM(I4(P)),ITRAN,IVEC) 678 FAACON(P) = WORK(KAFACONS-1 + (ITRAN-1)*MXVEC1 + IVEC) 679 680 CALL CCQR_SETFA(WORK(KAFATRAN),WORK(KAFADOTS),MXTRAN3,MXVEC1, 681 & IZT(I2(P)),IOP(I1(P)),IAM(I3(P)),IAM(I4(P)),ITRAN,IVEC) 682 FAACON(6+P) = WORK(KAFACONS-1 + (ITRAN-1)*MXVEC1 + IVEC) 683 END DO 684 685*---------------------------------------------------------------------* 686* get ETA{O} vector calculations, 12 permutations 687*---------------------------------------------------------------------* 688 DO P = 1, 6 689 IF ((USE_L2BC.AND.P.EQ.4) .OR. (USE_LBCD.AND.P.GT.3) ) THEN 690 EAACON(P) = ZERO 691 EAACON(6+P) = ZERO 692 ELSE 693 CALL CC_SETXE('Eta',WORK(KAEATRAN),WORK(KAEADOTS), !1x2x3,4 694 & MXTRAN2,MXVEC2, 695 & IZT(I1(P)),IOP(I2(P)),IKAP(I2(P)),0,0,0, 696 & IR2(IP(P)),ITRAN,IVEC) 697 EAACON(P) = WORK(KAEACONS-1 + (ITRAN-1)*MXVEC2 + IVEC) 698 699 CALL CC_SETXE('Eta',WORK(KAEATRAN),WORK(KAEADOTS), !2x1x3,4 700 & MXTRAN2,MXVEC2, 701 & IZT(I2(P)),IOP(I1(P)),IKAP(I1(P)),0,0,0, 702 & IR2(IP(P)),ITRAN,IVEC) 703 EAACON(6+P) = WORK(KAEACONS-1 + (ITRAN-1)*MXVEC2 + IVEC) 704 END IF 705 END DO 706 707*---------------------------------------------------------------------* 708* get L2 x O2 vector dot products, max. 3 permutations 709*---------------------------------------------------------------------* 710 IF (USE_LBCD .OR. USE_L2BC) THEN 711 CALL CC_SETDOT(WORK(KLTRAN),WORK(KLDOTS),MXTRAN2,MXVEC2, 712 & IL2(BC),IO2(AD),ITRAN,IVEC) 713 LCON(1) = WORK(KLCONS-1 + (ITRAN-1)*MXVEC2 + IVEC) 714 ELSE 715 LCON(1) = ZERO 716 END IF 717 718 IF (USE_LBCD) THEN 719 CALL CC_SETDOT(WORK(KLTRAN),WORK(KLDOTS),MXTRAN2,MXVEC2, 720 & IL2(BD),IO2(AC),ITRAN,IVEC) 721 LCON(2) = WORK(KLCONS-1 + (ITRAN-1)*MXVEC2 + IVEC) 722 723 CALL CC_SETDOT(WORK(KLTRAN),WORK(KLDOTS),MXTRAN2,MXVEC2, 724 & IL2(CD),IO2(AB),ITRAN,IVEC) 725 LCON(3) = WORK(KLCONS-1 + (ITRAN-1)*MXVEC2 + IVEC) 726 ELSE 727 LCON(2) = ZERO 728 LCON(3) = ZERO 729 END IF 730*---------------------------------------------------------------------* 731* get chi vector dot products, 6 permutations 732*---------------------------------------------------------------------* 733 IF (L_USE_CHI2) THEN 734 CALL CC_SETDOT(WORK(KXTRAN),WORK(KXDOTS),MXTRAN2,MXVEC2, 735 & IX2(AB),IR2(CD),ITRAN,IVEC) 736 XCON(1) = WORK(KXCONS-1 + (ITRAN-1)*MXVEC2 + IVEC) 737 738 CALL CC_SETDOT(WORK(KXTRAN),WORK(KXDOTS),MXTRAN2,MXVEC2, 739 & IX2(AC),IR2(BD),ITRAN,IVEC) 740 XCON(2) = WORK(KXCONS-1 + (ITRAN-1)*MXVEC2 + IVEC) 741 742 CALL CC_SETDOT(WORK(KXTRAN),WORK(KXDOTS),MXTRAN2,MXVEC2, 743 & IX2(AD),IR2(BC),ITRAN,IVEC) 744 XCON(3) = WORK(KXCONS-1 + (ITRAN-1)*MXVEC2 + IVEC) 745 746 CALL CC_SETDOT(WORK(KXTRAN),WORK(KXDOTS),MXTRAN2,MXVEC2, 747 & IX2(CD),IR2(AB),ITRAN,IVEC) 748 XCON(4) = WORK(KXCONS-1 + (ITRAN-1)*MXVEC2 + IVEC) 749 750 CALL CC_SETDOT(WORK(KXTRAN),WORK(KXDOTS),MXTRAN2,MXVEC2, 751 & IX2(BD),IR2(AC),ITRAN,IVEC) 752 XCON(5) = WORK(KXCONS-1 + (ITRAN-1)*MXVEC2 + IVEC) 753 754 CALL CC_SETDOT(WORK(KXTRAN),WORK(KXDOTS),MXTRAN2,MXVEC2, 755 & IX2(BC),IR2(AD),ITRAN,IVEC) 756 XCON(6) = WORK(KXCONS-1 + (ITRAN-1)*MXVEC2 + IVEC) 757 ELSE 758 XCON(1) = ZERO 759 XCON(2) = ZERO 760 XCON(3) = ZERO 761 XCON(4) = ZERO 762 XCON(5) = ZERO 763 XCON(6) = ZERO 764 END IF 765 766 767*---------------------------------------------------------------------* 768* get xksi vector dot products, 4 permutations 769*---------------------------------------------------------------------* 770 IF (L_USE_XKS3) THEN 771 CALL CC_SETDOT(WORK(KOTRAN),WORK(KODOTS),MXTRAN3,MXVEC1, 772 & IO3(ABC),IZT(D),ITRAN,IVEC) 773 OCON(1) = WORK(KOCONS-1 + (ITRAN-1)*MXVEC1 + IVEC) 774 775 CALL CC_SETDOT(WORK(KOTRAN),WORK(KODOTS),MXTRAN3,MXVEC1, 776 & IO3(ABD),IZT(C),ITRAN,IVEC) 777 OCON(2) = WORK(KOCONS-1 + (ITRAN-1)*MXVEC1 + IVEC) 778 779 CALL CC_SETDOT(WORK(KOTRAN),WORK(KODOTS),MXTRAN3,MXVEC1, 780 & IO3(ACD),IZT(B),ITRAN,IVEC) 781 OCON(3) = WORK(KOCONS-1 + (ITRAN-1)*MXVEC1 + IVEC) 782 783 CALL CC_SETDOT(WORK(KOTRAN),WORK(KODOTS),MXTRAN3,MXVEC1, 784 & IO3(BCD),IZT(A),ITRAN,IVEC) 785 OCON(4) = WORK(KOCONS-1 + (ITRAN-1)*MXVEC1 + IVEC) 786 ELSE 787 OCON(1) = ZERO 788 OCON(2) = ZERO 789 OCON(3) = ZERO 790 OCON(4) = ZERO 791 END IF 792 793*---------------------------------------------------------------------* 794* add contributions to hyperpolarizabilities: 795*---------------------------------------------------------------------* 796 797 IDX = NCRFREQ*(IOPER-1)+ IFREQ 798 799 IF (ISIGN.EQ.-1) IDX = IDX + NCRFREQ*NCROPER 800 801 WORK(KHYPPOL-1+IDX) = WORK(KHYPPOL-1+IDX) + 802 & H0CON + 803 & G0CON(1) +G0CON(2) +G0CON(3) +G0CON(4) +G0CON(5) +G0CON(6) + 804 & GACON(1) +GACON(2) +GACON(3) +GACON(4) + 805 & F0CON(1) +F0CON(2) +F0CON(3) + 806 & FACON(1) +FACON(2) +FACON(3) +FACON(4) +FACON(5) +FACON(6) + 807 & FACON(7) +FACON(8) +FACON(9) +FACON(10) +FACON(11) +FACON(12) + 808 & F0ACON(1)+F0ACON(2)+F0ACON(3)+F0ACON(4) +F0ACON(5) +F0ACON(6) + 809 & F0ACON(7)+F0ACON(8)+F0ACON(9)+F0ACON(10)+F0ACON(11)+F0ACON(12)+ 810 & FAACON(1)+FAACON(2)+FAACON(3)+FAACON(4) +FAACON(5) +FAACON(6) + 811 & FAACON(7)+FAACON(8)+FAACON(9)+FAACON(10)+FAACON(11)+FAACON(12)+ 812 & EAACON(1)+EAACON(2)+EAACON(3)+EAACON(4) +EAACON(5) +EAACON(6) + 813 & EAACON(7)+EAACON(8)+EAACON(9)+EAACON(10)+EAACON(11)+EAACON(12)+ 814 & XCON(1) +XCON(2) +XCON(3) +XCON(4) +XCON(5) +XCON(6) + 815 & OCON(1) +OCON(2) +OCON(3) +OCON(4) + 816 & LCON(1) +LCON(2) +LCON(3) 817 818 IF (LOCDBG) THEN 819 WRITE(LUPRI,'(A,3I5)') 'IOPER, IFREQ, ISIGN:',IOPER,IFREQ,ISIGN 820 WRITE(LUPRI,'(A,4I5)') 'IZT:',(IZT(I),I=1,4) 821 WRITE(LUPRI,'(A,4I5)') 'IAM:',(IAM(I),I=1,4) 822 WRITE(LUPRI,'(A,6I5)') 'IR2:',(IR2(I),I=1,6) 823 IF (L_USE_CHI2) WRITE(LUPRI,'(A,6I5)') 'IX2:',(IX2(I),I=1,6) 824 IF (L_USE_XKS3) WRITE(LUPRI,'(A,6I5)') 'IO3:',(IO3(I),I=1,4) 825 WRITE(LUPRI,*) 'H0CON: ',H0CON 826 WRITE(LUPRI,*) 'G0CON: ',(G0CON(I),I=1,6) 827 WRITE(LUPRI,*) 'GACON: ',(GACON(I),I=1,4) 828 WRITE(LUPRI,*) 'F0CON: ',(F0CON(I),I=1,3) 829 WRITE(LUPRI,*) 'FACON: ',(FACON(I),I=1,12) 830 WRITE(LUPRI,*) 'F0ACON:',(F0ACON(I),I=1,12) 831 WRITE(LUPRI,*) 'FAACON:',(FAACON(I),I=1,12) 832 WRITE(LUPRI,*) 'EAACON:',(EAACON(I),I=1,12) 833 WRITE(LUPRI,*) 'XCON:',(XCON(I),I=1,6) 834 WRITE(LUPRI,*) 'OCON:',(OCON(I),I=1,4) 835 WRITE(LUPRI,*) 'LCON:',(LCON(I),I=1,3) 836 END IF 837 838 END DO 839 END DO 840 END IF 841 END DO 842 843*---------------------------------------------------------------------* 844* print frequency-dependent hyperpolarizabilities 845*---------------------------------------------------------------------* 846 WRITE (LUPRI,'(/A,I4,A,F12.2," seconds.")') ' Total time for', 847 & NBHYPPOL,' cubic response functions: ', SECOND() - TIM0 848 849 CALL CCCHYPPRT(WORK(KHYPPOL)) 850 851*---------------------------------------------------------------------* 852* now we are finished with frequency-dependent hyperpolarizabilities... 853*---------------------------------------------------------------------* 854 END IF ! (NCRFREQ.GT.0) 855 856 IF (NCRDISP.EQ.0) THEN 857 RETURN 858 END IF 859 860*---------------------------------------------------------------------* 861* allocate & initialize work space for expansion coefficients 862*---------------------------------------------------------------------* 863 NBEXPCOF = NCROPER * NCRDISP 864 865 KEXPCOF = 1 866 KEND0 = KEXPCOF + NBEXPCOF 867 KEND0A = KEND0 868 869 870 MXTRAN2 = 2 * NLRCLBL * NLRCLBL 871 MXVEC2 = (NLRCLBL+1) * NLRCLBL /2 872 MXTRAN3 = 2 * NLRCLBL * NLRCLBL * NLRCLBL 873 MXVEC1 = NLRCLBL 874 MXTRAN2 = MIN(MXTRAN2,2*12*NBEXPCOF) ! 12 = max. # of permut. 875 MXTRAN3 = MIN(MXTRAN3,2*12*NBEXPCOF) ! 12 = max. # of permut. 876 MXVEC2 = MIN(MXVEC2,2*12*NBEXPCOF) ! 12 = max. # of permut. 877 MXVEC1 = MIN(MXVEC1,2*12*NBEXPCOF) ! 12 = max. # of permut. 878 879 MX0HTRAN = 5 * MXTRAN3 880 MX0GTRAN = 4 * MXTRAN2 881 MXAGTRAN = 4 * MXTRAN3 882 MX0FTRAN = 3 * MXTRAN2 883 MXAFTRAN = 3 * MXTRAN2 884 MX0FATRAN = 5 * MXTRAN2 885 MXAFATRAN = 5 * MXTRAN3 886 MXAEATRAN = 3 * MXTRAN2 887 MXXTRAN = 1 * MXTRAN2 888 MXOTRAN = 1 * MXTRAN2 889 MXLTRAN = 1 * MXTRAN2 890 891 MX0HDOTS = MXTRAN3 * MXVEC1 892 MX0GDOTS = MXTRAN2 * MXVEC2 893 MXAGDOTS = MXTRAN3 * MXVEC1 894 MX0FDOTS = MXTRAN2 * MXVEC2 895 MXAFDOTS = MXTRAN2 * MXVEC2 896 MX0FADOTS = MXTRAN2 * MXVEC2 897 MXAFADOTS = MXTRAN3 * MXVEC1 898 MXAEADOTS = MXTRAN2 * MXVEC2 899 MXXDOTS = MXTRAN2 * MXVEC2 900 MXODOTS = MXTRAN2 * MXVEC2 901 MXLDOTS = MXTRAN2 * MXVEC2 902 903 K0HTRAN = KEND0 904 K0HDOTS = K0HTRAN + MX0HTRAN 905 K0HCONS = K0HDOTS + MX0HDOTS 906 KEND0 = K0HCONS + MX0HDOTS 907 908 K0GTRAN = KEND0 909 K0GDOTS = K0GTRAN + MX0GTRAN 910 K0GCONS = K0GDOTS + MX0GDOTS 911 KEND0 = K0GCONS + MX0GDOTS 912 913 KAGTRAN = KEND0 914 KAGDOTS = KAGTRAN + MXAGTRAN 915 KAGCONS = KAGDOTS + MXAGDOTS 916 KEND0 = KAGCONS + MXAGDOTS 917 918 K0FTRAN = KEND0 919 K0FDOTS = K0FTRAN + MX0FTRAN 920 K0FCONS = K0FDOTS + MX0FDOTS 921 KEND0 = K0FCONS + MX0FDOTS 922 923 KAFTRAN = KEND0 924 KAFDOTS = KAFTRAN + MXAFTRAN 925 KAFCONS = KAFDOTS + MXAFDOTS 926 KEND0 = KAFCONS + MXAFDOTS 927 928 K0FATRAN = KEND0 929 K0FADOTS = K0FATRAN + MX0FATRAN 930 K0FACONS = K0FADOTS + MX0FADOTS 931 KEND0 = K0FACONS + MX0FADOTS 932 933 KAFATRAN = KEND0 934 KAFADOTS = KAFATRAN + MXAFATRAN 935 KAFACONS = KAFADOTS + MXAFADOTS 936 KEND0 = KAFACONS + MXAFADOTS 937 938 KAEATRAN = KEND0 939 KAEADOTS = KAEATRAN + MXAEATRAN 940 KAEACONS = KAEADOTS + MXAEADOTS 941 KEND0 = KAEACONS + MXAEADOTS 942 943 KXTRAN = KEND0 944 KXDOTS = KXTRAN + MXXTRAN 945 KXCONS = KXDOTS + MXXDOTS 946 KEND0 = KXCONS + MXXDOTS 947 948 KOTRAN = KEND0 949 KODOTS = KOTRAN + MXOTRAN 950 KOCONS = KODOTS + MXODOTS 951 KEND0 = KOCONS + MXODOTS 952 953 KLTRAN = KEND0 954 KLDOTS = KLTRAN + MXLTRAN 955 KLCONS = KLDOTS + MXLDOTS 956 KEND0 = KLCONS + MXLDOTS 957 958 LEND0 = LWORK - KEND0 959 960 IF (LEND0 .LT. 0) THEN 961 WRITE (LUPRI,*) 'Insufficient memory in CC_2HYP. (2)' 962 WRITE (LUPRI,*) 'available:',LWORK 963 WRITE (LUPRI,*) ' needed:',KEND0 964 WRITE (LUPRI,*) ' MXTRAN2:',MXTRAN2 965 WRITE (LUPRI,*) ' MXTRAN3:',MXTRAN3 966 WRITE (LUPRI,*) ' MXVEC1 :',MXVEC1 967 WRITE (LUPRI,*) ' MXVEC2 :',MXVEC2 968 CALL QUIT('Insufficient memory in CC_2HYP. (2)') 969 END IF 970 971 972* initialize expansion coefficients and the lists of contributions: 973 CALL DZERO(WORK(KEXPCOF),NBEXPCOF) 974 CALL DZERO(WORK(K0HTRAN),MX0HTRAN+2*MX0HDOTS) 975 CALL DZERO(WORK(K0GTRAN),MX0GTRAN+2*MX0GDOTS) 976 CALL DZERO(WORK(KAGTRAN),MXAGTRAN+2*MXAGDOTS) 977 CALL DZERO(WORK(K0FTRAN),MX0FTRAN+2*MX0FDOTS) 978 CALL DZERO(WORK(KAFTRAN),MXAFTRAN+2*MXAFDOTS) 979 CALL DZERO(WORK(K0FATRAN),MX0FATRAN+2*MX0FADOTS) 980 CALL DZERO(WORK(KAFATRAN),MXAFATRAN+2*MXAFADOTS) 981 CALL DZERO(WORK(KAEATRAN),MXAEATRAN+2*MXAEADOTS) 982 CALL DZERO(WORK(KXTRAN),MXXTRAN+2*MXXDOTS) 983 CALL DZERO(WORK(KOTRAN),MXOTRAN+2*MXODOTS) 984 CALL DZERO(WORK(KLTRAN),MXLTRAN+2*MXLDOTS) 985 986* initialize timing: 987 TIM2 = SECOND() 988 989* save options L_USE_CHI2 and L_USE_XKS3, which are not implemented 990* for the dispersion coefficients: 991 LCHI2_SAVE = L_USE_CHI2 992 LXKS3_SAVE = L_USE_XKS3 993 L_USE_CHI2 = .FALSE. 994 L_USE_XKS3 = .FALSE. 995 996*---------------------------------------------------------------------* 997* set up lists for H, G, F and F{O} transformations and ETA{O} vectors: 998*---------------------------------------------------------------------* 999 LADD = .FALSE. 1000 1001 CALL CCCR_DISP_SETUP(MXTRAN2, MXVEC2, MXTRAN3, MXVEC1, 1002 & WORK(K0HTRAN), WORK(K0HDOTS), WORK(K0HCONS), N0HTRAN, 1003 & WORK(K0GTRAN), WORK(K0GDOTS), WORK(K0GCONS), N0GTRAN, 1004 & WORK(KAGTRAN), WORK(KAGDOTS), WORK(KAGCONS), NAGTRAN, 1005 & WORK(K0FTRAN), WORK(K0FDOTS), WORK(K0FCONS), N0FTRAN, 1006 & WORK(KAFTRAN), WORK(KAFDOTS), WORK(KAFCONS), NAFTRAN, 1007 & WORK(K0FATRAN),WORK(K0FADOTS),WORK(K0FACONS),N0FATRAN, 1008 & WORK(KAFATRAN),WORK(KAFADOTS),WORK(KAFACONS),NAFATRAN, 1009 & WORK(KAEATRAN),WORK(KAEADOTS),WORK(KAEACONS),NAEATRAN, 1010 & WORK(KXTRAN), WORK(KXDOTS), WORK(KXCONS), NXTRAN, 1011 & WORK(KOTRAN), WORK(KODOTS), WORK(KOCONS), NOTRAN, 1012 & WORK(KLTRAN), WORK(KLDOTS), WORK(KLCONS), NLTRAN, 1013 & WORK(KEXPCOF), NBEXPCOF, LADD ) 1014 1015*---------------------------------------------------------------------* 1016* calculate H matrix contributions: 1017*---------------------------------------------------------------------* 1018 TIM1 = SECOND() 1019 1020 IOPT = 5 1021 CALL CC_HMAT('L0 ','RC ','RC ','RC ','RC ',N0HTRAN, MXVEC1, 1022 & WORK(K0HTRAN),WORK(K0HDOTS),WORK(K0HCONS), 1023 & WORK(KEND0), LEND0, IOPT ) 1024 1025 TIMH = SECOND() - TIM1 1026 WRITE (LUPRI,'(/A,I5,A,F12.2," seconds.")') 1027 & ' Time used for',N0HTRAN,' H matrix transformations: ',TIMH 1028 CALL FLSHFO(LUPRI) 1029 1030*---------------------------------------------------------------------* 1031* calculate G matrix contributions: 1032*---------------------------------------------------------------------* 1033 TIM1 = SECOND() 1034 1035 IF ((.NOT.L_USE_CHI2) .AND. NO_2NP1_RULE) THEN 1036 IOPT = 5 1037 CALL CC_GMATRIX('L0 ','RC ','RC ','CR2',N0GTRAN, MXVEC2, 1038 & WORK(K0GTRAN),WORK(K0GDOTS),WORK(K0GCONS), 1039 & WORK(KEND0), LEND0, IOPT ) 1040 END IF 1041 1042 IF (.NOT.L_USE_XKS3) THEN 1043 IOPT = 5 1044 CALL CC_GMATRIX('LC ','RC ','RC ','RC ',NAGTRAN, MXVEC1, 1045 & WORK(KAGTRAN),WORK(KAGDOTS),WORK(KAGCONS), 1046 & WORK(KEND0), LEND0, IOPT ) 1047 END IF 1048 1049 TIMG = SECOND() - TIM1 1050 WRITE (LUPRI,'(/A,I5,A,F12.2," seconds.")') 1051 & ' Time used for',N0GTRAN+NAGTRAN, 1052 & ' G matrix transformations: ',TIMG 1053 CALL FLSHFO(LUPRI) 1054 1055*---------------------------------------------------------------------* 1056* calculate F matrix contributions: 1057*---------------------------------------------------------------------* 1058 TIM1 = SECOND() 1059 1060 IOPT = 5 1061 CALL CC_FMATRIX(WORK(K0FTRAN),N0FTRAN,'L0 ','CR2',IOPT,'CR2', 1062 & WORK(K0FDOTS),WORK(K0FCONS),MXVEC2, 1063 & WORK(KEND0),LEND0) 1064 1065 IF (L_USE_CHI2) THEN 1066 CALL DSCAL(MX0FDOTS,-1.0d0,WORK(K0FCONS),1) 1067 END IF 1068 1069 IF ((.NOT. (L_USE_CHI2 .OR. L_USE_XKS3)) .AND. NO_2NP1_RULE) THEN 1070 IOPT = 5 1071 CALL CC_FMATRIX(WORK(KAFTRAN),NAFTRAN,'LC ','RC ',IOPT,'CR2', 1072 & WORK(KAFDOTS),WORK(KAFCONS),MXVEC2, 1073 & WORK(KEND0),LEND0) 1074 END IF 1075 1076 TIMF = SECOND() - TIM1 1077 WRITE (LUPRI,'(/A,I5,A,F12.2," seconds.")') 1078 & ' Time used for',N0FTRAN+NAFTRAN, 1079 & ' F matrix transformations: ',TIMF 1080 CALL FLSHFO(LUPRI) 1081 1082*---------------------------------------------------------------------* 1083* calculate F{O} matrix contributions: 1084*---------------------------------------------------------------------* 1085 TIM1 = SECOND() 1086 1087 IF ((.NOT.L_USE_CHI2) .AND. NO_2NP1_RULE) THEN 1088 CALL CCQR_FADRV('L0 ','o1 ','RC ','CR2',N0FATRAN, MXVEC2, 1089 & WORK(K0FATRAN),WORK(K0FADOTS), 1090 & WORK(K0FACONS), 1091 & WORK(KEND0), LEND0, 'DOTP' ) 1092 END IF 1093 1094 IF (.NOT.L_USE_XKS3) THEN 1095 CALL CCQR_FADRV('LC ','o1 ','RC ','RC ',NAFATRAN, MXVEC1, 1096 & WORK(KAFATRAN),WORK(KAFADOTS), 1097 & WORK(KAFACONS), 1098 & WORK(KEND0), LEND0, 'DOTP' ) 1099 END IF 1100 1101 TIMFA = SECOND() - TIM1 1102 WRITE (LUPRI,'(/A,I5,A,F12.2," seconds.")') 1103 & ' Time used for',N0FATRAN+NAFATRAN, 1104 & ' F{O} matrix transformations:',TIMFA 1105 CALL FLSHFO(LUPRI) 1106 1107*---------------------------------------------------------------------* 1108* calculate ETA{O} vector contributions: 1109*---------------------------------------------------------------------* 1110 TIM1 = SECOND() 1111 1112 IF ((.NOT.(L_USE_CHI2 .OR. L_USE_XKS3)) .AND. NO_2NP1_RULE ) THEN 1113 CALL CCQR_EADRV('LC ','o1 ','CR2',NAEATRAN, MXVEC2, 1114 & WORK(KAEATRAN),WORK(KAEADOTS),WORK(KAEACONS), 1115 & WORK(KEND0), LEND0, 'DOTP' ) 1116 END IF 1117 1118 TIMEA = SECOND() - TIM1 1119 WRITE (LUPRI,'(/A,I5,A,F12.2," seconds.")') ' Time used for', 1120 & NAEATRAN,' ETA{O} vector calculations :',TIMEA 1121 CALL FLSHFO(LUPRI) 1122 1123*---------------------------------------------------------------------* 1124* chi vector contributions: 1125*---------------------------------------------------------------------* 1126 IF (.NOT. NO_2NP1_RULE) THEN 1127 TIM1 = SECOND() 1128 1129 CALL CC_DOTDRV('CX2','CR2',NXTRAN,MXVEC2, 1130 & WORK(KXTRAN), WORK(KXDOTS), WORK(KXCONS), 1131 & WORK(KEND0), LEND0 ) 1132 1133 TIMX2 = SECOND() - TIM1 1134 WRITE (LUPRI,'(/A,I5,A,5X,F12.2," seconds.")') 1135 & ' Time used for',NXTRAN, ' CX2 x CR2 dot products:',TIMX2 1136 CALL FLSHFO(LUPRI) 1137 1138 END IF 1139*---------------------------------------------------------------------* 1140* CL2 x CR2 dot products: 1141*---------------------------------------------------------------------* 1142 IF (.NOT. NO_2NP1_RULE) THEN 1143 TIM1 = SECOND() 1144 1145 CALL CC_DOTDRV('CL2','CR2',NLTRAN,MXVEC2, 1146 & WORK(KLTRAN), WORK(KLDOTS), WORK(KLCONS), 1147 & WORK(KEND0), LEND0 ) 1148 1149 TIML2 = SECOND() - TIM1 1150 WRITE (LUPRI,'(/A,I5,A,5X,F12.2," seconds.")') 1151 & ' Time used for',NLTRAN, ' CL2 x CR2 dot products:',TIML2 1152 CALL FLSHFO(LUPRI) 1153 1154 END IF 1155*---------------------------------------------------------------------* 1156* xksi vector contributions: 1157*---------------------------------------------------------------------* 1158 IF (.NOT. NO_2NP1_RULE) THEN 1159 1160 TIM1 = SECOND() 1161 1162 CALL CC_DOTDRV('CO2','CL2',NOTRAN,MXVEC2, 1163 & WORK(KOTRAN), WORK(KODOTS), WORK(KOCONS), 1164 & WORK(KEND0), LEND0 ) 1165 1166 TIMO2 = SECOND() - TIM1 1167 WRITE (LUPRI,'(/A,I5,A,5X,F12.2," seconds.")') 1168 & ' Time used for',NOTRAN, ' CL2 x CO2 dot products:',TIMO2 1169 CALL FLSHFO(LUPRI) 1170 1171 END IF 1172 1173*---------------------------------------------------------------------* 1174* collect contributions and add them to hyperpolarizabilities 1175*---------------------------------------------------------------------* 1176 LADD = .TRUE. 1177 1178 CALL CCCR_DISP_SETUP(MXTRAN2, MXVEC2, MXTRAN3, MXVEC1, 1179 & WORK(K0HTRAN), WORK(K0HDOTS), WORK(K0HCONS), N0HTRAN, 1180 & WORK(K0GTRAN), WORK(K0GDOTS), WORK(K0GCONS), N0GTRAN, 1181 & WORK(KAGTRAN), WORK(KAGDOTS), WORK(KAGCONS), NAGTRAN, 1182 & WORK(K0FTRAN), WORK(K0FDOTS), WORK(K0FCONS), N0FTRAN, 1183 & WORK(KAFTRAN), WORK(KAFDOTS), WORK(KAFCONS), NAFTRAN, 1184 & WORK(K0FATRAN),WORK(K0FADOTS),WORK(K0FACONS),N0FATRAN, 1185 & WORK(KAFATRAN),WORK(KAFADOTS),WORK(KAFACONS),NAFATRAN, 1186 & WORK(KAEATRAN),WORK(KAEADOTS),WORK(KAEACONS),NAEATRAN, 1187 & WORK(KXTRAN), WORK(KXDOTS), WORK(KXCONS), NXTRAN, 1188 & WORK(KOTRAN), WORK(KODOTS), WORK(KOCONS), NOTRAN, 1189 & WORK(KLTRAN), WORK(KLDOTS), WORK(KLCONS), NLTRAN, 1190 & WORK(KEXPCOF), NBEXPCOF, LADD ) 1191 1192 1193* recover options L_USE_CHI2 and L_USE_XKS3, which are not implemented 1194* for the dispersion coefficients: 1195 L_USE_CHI2 = LCHI2_SAVE 1196 L_USE_XKS3 = LXKS3_SAVE 1197 1198*---------------------------------------------------------------------* 1199* print timing: 1200*---------------------------------------------------------------------* 1201 WRITE (LUPRI,'(/A,I5,A,F12.2," seconds.")') ' Total time for', 1202 & NBEXPCOF,' expansion coefficients: ', SECOND() - TIM2 1203 1204*---------------------------------------------------------------------* 1205* print expansion coefficients for second hyperpolarizabilities: 1206*---------------------------------------------------------------------* 1207 IF (IPRCHYP .GE. 15) THEN 1208 CALL CCCEXPPRT(WORK(KEXPCOF)) 1209 END IF 1210 1211 CALL FLSHFO(LUPRI) 1212 1213*---------------------------------------------------------------------* 1214* calculate from the d_ABCD expansion coeff. the D_ABCD disp coeff.: 1215*---------------------------------------------------------------------* 1216 IF (NCRDSPE.GE.0) THEN 1217 NTRINOM = (NCRDSPE+3)*(NCRDSPE+2)*(NCRDSPE+1)/6 1218 1219 KTRINOM = KEND0A 1220 KDSPCF = KTRINOM + NTRINOM 1221 KTHGCF = KDSPCF + NTRINOM * NCROPER 1222 KESHGCF = KTHGCF + (NCRDSPE+1) * NCROPER 1223 KKERRCF = KESHGCF + (NCRDSPE+1) * NCROPER 1224 KDFWMCF = KKERRCF + (NCRDSPE+1) * NCROPER 1225 KEND1 = KDFWMCF + (NCRDSPE+1) * NCROPER 1226 LEND1 = LWORK - KEND1 1227 1228 IF (GAMMA_PAR .OR. GAMMA_ORT) THEN 1229 KDSPCFA = KEND1 1230 KTHGCFA = KDSPCFA + NTRINOM * 4 1231 KESHGCA = KTHGCFA + (NCRDSPE+1) * 4 1232 KKERRCA = KESHGCA + (NCRDSPE+1) * 4 1233 KDFWMCA = KKERRCA + (NCRDSPE+1) * 4 1234 KABCDE = KDFWMCA + (NCRDSPE+1) * 4 1235 KEND1 = KABCDE + 15 * 2 1236 LEND1 = LWORK - KEND1 1237 END IF 1238 1239 IF (LEND1 .LT. 0) THEN 1240 CALL QUIT('Insufficient memory in CC_2HYP. (3)') 1241 END IF 1242 1243 CALL DZERO(WORK(KTRINOM),KEND1-KTRINOM) 1244 1245 CALL CCCRDISP(WORK(KDSPCF),WORK(KDSPCFA),WORK(KABCDE), 1246 & WORK(KTHGCF),WORK(KESHGCF), 1247 & WORK(KKERRCF),WORK(KDFWMCF), 1248 & WORK(KTHGCFA),WORK(KESHGCA), 1249 & WORK(KKERRCA),WORK(KDFWMCA), 1250 & WORK(KEXPCOF),WORK(KTRINOM), NTRINOM,LERROR) 1251 1252 CALL CCCDISPRT(WORK(KDSPCF),WORK(KDSPCFA),WORK(KABCDE), 1253 & WORK(KTHGCF),WORK(KESHGCF), 1254 & WORK(KKERRCF),WORK(KDFWMCF), 1255 & WORK(KTHGCFA),WORK(KESHGCA), 1256 & WORK(KKERRCA),WORK(KDFWMCA),LERROR) 1257 1258 END IF 1259 1260 RETURN 1261 END 1262 1263*=====================================================================* 1264* END OF SUBROUTINE CC_2HYP * 1265*=====================================================================* 1266 SUBROUTINE CCCHYPPRT(HYPERPOL) 1267*---------------------------------------------------------------------* 1268* 1269* Purpose: print second hyperpolarizabilities 1270* 1271* 1272* Written by Christof Haettig in February 1997. 1273* 1274*=====================================================================* 1275#if defined (IMPLICIT_NONE) 1276 IMPLICIT NONE 1277#else 1278# include "implicit.h" 1279#endif 1280#include "priunit.h" 1281#include "ccorb.h" 1282#include "ccsdinp.h" 1283#include "cccrinf.h" 1284#include "ccroper.h" 1285 1286 1287 CHARACTER*10 BLANKS 1288 CHARACTER*80 STRING 1289 INTEGER ISYMA, ISYMB, ISYMC, ISYMD 1290 INTEGER IFREQ, IOPER, ISYOLD, ICOMP, IADD 1291 1292 1293#if defined (SYS_CRAY) 1294 REAL HYPERPOL(NCRFREQ,NCROPER,2) 1295 REAL HALF, GAMMA, PROP 1296#else 1297 DOUBLE PRECISION HYPERPOL(NCRFREQ,NCROPER,2) 1298 DOUBLE PRECISION HALF, GAMMA, PROP 1299#endif 1300 PARAMETER (HALF = 0.5D0) 1301 1302 CHARACTER MODEL*10,MOPRPC 1303 INTEGER ISYMABCD 1304*---------------------------------------------------------------------* 1305* print header for hyperpolarizability section 1306*---------------------------------------------------------------------* 1307 BLANKS = ' ' 1308 STRING = ' RESULTS FOR THE SECOND HYPERPOLARIZABILITIES ' 1309 1310 IF (CCS) THEN 1311 CALL AROUND( BLANKS//'FINAL CCS'//STRING(1:45)//BLANKS ) 1312 ELSE IF (CC2) THEN 1313 CALL AROUND( BLANKS//'FINAL CC2'//STRING(1:45)//BLANKS ) 1314 ELSE IF (CCSD) THEN 1315 CALL AROUND( BLANKS//'FINAL CCSD'//STRING(1:45)//BLANKS ) 1316 ELSE IF (CC3) THEN 1317 CALL AROUND( BLANKS//'FINAL CC3'//STRING(1:45)//BLANKS ) 1318 ELSE 1319 CALL QUIT('CCCHYPPRT called for an unknown Coupled '// 1320 & 'Cluster model.') 1321 END IF 1322 IF (CCSD) THEN 1323 MODEL = 'CCSD' 1324 MOPRPC = 'CCSD ' 1325 ELSE IF (CIS) THEN 1326 MODEL = 'CIS' 1327 MOPRPC = 'CIS ' 1328 ELSE IF (CCS) THEN 1329 MODEL = 'CCS' 1330 MOPRPC = 'CCS ' 1331 ELSE IF (CC2) THEN 1332 MODEL = 'CC2' 1333 MOPRPC = 'CC2 ' 1334 ELSE IF (CC3) THEN 1335 MODEL = 'CC3' 1336 MOPRPC = 'CC3 ' 1337 END IF 1338 1339 IF (IPRCHYP.GT.5) THEN 1340 WRITE(LUPRI,'(/1X,4(1X,A," operator",7X),4X,A,8X,A,/,107("-"))') 1341 & 'A','B','C','D','gamma','(asy. Resp.)' 1342 ELSE 1343 WRITE(LUPRI,'(/1X,4(1X,A," operator",7X),4X,A,/,86("-"))') 1344 & 'A','B','C','D','gamma' 1345 END IF 1346 1347 ISYOLD = 1 1348 1349 DO IOPER = 1, NCROPER 1350 ISYMA = ISYOPR(IACROP(IOPER)) 1351 ISYMB = ISYOPR(IBCROP(IOPER)) 1352 ISYMC = ISYOPR(ICCROP(IOPER)) 1353 ISYMD = ISYOPR(IDCROP(IOPER)) 1354 ISYMABCD = MULD2H(MULD2H(ISYMA,ISYMB),MULD2H(ISYMC,ISYMD)) 1355 1356 1357 IFREQ = 1 1358cmbh initialize 1359 PROP=0.0d0 1360cmbh end 1361 IF (MULD2H(ISYMA,ISYMB).EQ.MULD2H(ISYMC,ISYMD)) THEN 1362 PROP = HALF*(HYPERPOL(IFREQ,IOPER,1)+HYPERPOL(IFREQ,IOPER,2)) 1363 IF (IPRCHYP.GT.5) THEN 1364 WRITE(LUPRI,'(/1X,4(A8,F7.4,3X),G18.10," (",G18.10,")")') 1365 & LBLOPR(IACROP(IOPER)),ACRFR(IFREQ), 1366 & LBLOPR(IBCROP(IOPER)),BCRFR(IFREQ), 1367 & LBLOPR(ICCROP(IOPER)),CCRFR(IFREQ), 1368 & LBLOPR(IDCROP(IOPER)),DCRFR(IFREQ), 1369 & HALF*(HYPERPOL(IFREQ,IOPER,1)+HYPERPOL(IFREQ,IOPER,2)), 1370 & HALF*(HYPERPOL(IFREQ,IOPER,1)-HYPERPOL(IFREQ,IOPER,2)) 1371 ELSE 1372 WRITE(LUPRI,'(/1X,4(A8,F7.4,3X),G16.8)') 1373 & LBLOPR(IACROP(IOPER)),ACRFR(IFREQ), 1374 & LBLOPR(IBCROP(IOPER)),BCRFR(IFREQ), 1375 & LBLOPR(ICCROP(IOPER)),CCRFR(IFREQ), 1376 & LBLOPR(IDCROP(IOPER)),DCRFR(IFREQ), 1377 & HALF*(HYPERPOL(IFREQ,IOPER,1)+HYPERPOL(IFREQ,IOPER,2)) 1378 END IF 1379 ISYOLD = 1 1380 ELSE 1381 IF (IPRCHYP.GT.5) THEN 1382 IF (ISYOLD.EQ.1) WRITE(LUPRI,*) 1383 WRITE(LUPRI,'(1X,4(A8,A7,3X),7X,A,9X," (",5X,A,6X,")")') 1384 & LBLOPR(IACROP(IOPER)),' -.- ', 1385 & LBLOPR(IBCROP(IOPER)),' -.- ', 1386 & LBLOPR(ICCROP(IOPER)),' -.- ', 1387 & LBLOPR(IDCROP(IOPER)),' -.- ', 1388 & '---', 1389 & '---' 1390 ELSE IF (IPRCHYP.GT.0) THEN 1391 IF (ISYOLD.EQ.1) WRITE(LUPRI,*) 1392 WRITE(LUPRI,'(1X,4(A8,A7,3X),6X,A,7X)') 1393 & LBLOPR(IACROP(IOPER)),' -.- ', 1394 & LBLOPR(IBCROP(IOPER)),' -.- ', 1395 & LBLOPR(ICCROP(IOPER)),' -.- ', 1396 & LBLOPR(IDCROP(IOPER)),' -.- ', 1397 & '---' 1398 END IF 1399 ISYOLD = 0 1400 END IF 1401 CALL WRIPRO(PROP,model,4, 1402 & LBLOPR(IACROP(IOPER)),LBLOPR(IBCROP(IOPER)), 1403 & LBLOPR(ICCROP(IOPER)),LBLOPR(IDCROP(IOPER)), 1404 & BCRFR(IFREQ),CCRFR(IFREQ),DCRFR(IFREQ),ISYMABCD, 1405 & 0,0,0) 1406 1407 1408 DO IFREQ = 2, NCRFREQ 1409 PROP = 0.0D0 1410 IF (MULD2H(ISYMA,ISYMB).EQ.MULD2H(ISYMC,ISYMD)) THEN 1411 PROP = HALF*(HYPERPOL(IFREQ,IOPER,1)+HYPERPOL(IFREQ,IOPER,2)) 1412 IF (IPRCHYP.GT.5) THEN 1413 WRITE(LUPRI,'(1X,4(8X,F7.4,3X),G18.10," (",G18.10,")")') 1414 & ACRFR(IFREQ), BCRFR(IFREQ), CCRFR(IFREQ), DCRFR(IFREQ), 1415 & HALF*(HYPERPOL(IFREQ,IOPER,1)+HYPERPOL(IFREQ,IOPER,2)), 1416 & HALF*(HYPERPOL(IFREQ,IOPER,1)-HYPERPOL(IFREQ,IOPER,2)) 1417 ELSE 1418 WRITE(LUPRI,'(1X,4(8X,F7.4,3X),G16.8)') 1419 & ACRFR(IFREQ), BCRFR(IFREQ), CCRFR(IFREQ), DCRFR(IFREQ), 1420 & HALF*(HYPERPOL(IFREQ,IOPER,1)+HYPERPOL(IFREQ,IOPER,2)) 1421 END IF 1422 END IF 1423cmbh 1424 CALL WRIPRO(PROP,model,4, 1425 & LBLOPR(IACROP(IOPER)),LBLOPR(IBCROP(IOPER)), 1426 & LBLOPR(ICCROP(IOPER)),LBLOPR(IDCROP(IOPER)), 1427 & BCRFR(IFREQ),CCRFR(IFREQ),DCRFR(IFREQ),ISYMABCD, 1428 & 0,0,0) 1429cmbh end 1430 END DO 1431 1432 END DO 1433 1434 WRITE(LUPRI,'(86("-"))') 1435 1436 1437 IF (GAMMA_PAR .OR. GAMMA_ORT) THEN 1438 WRITE(LUPRI,'(/////A,40X,A,/,75("-"))') 1439 & ' average frequencies','value' 1440 1441 IF (GAMMA_PAR) THEN 1442 1443* calculate & print gamma_{||} for all requested frequencies: 1444 DO IFREQ = 1, NCRFREQ 1445 1446 IF (CSYM(1:6).EQ.'ATOMIC') THEN 1447 GAMMA = HALF*(HYPERPOL(IFREQ,1,1)+HYPERPOL(IFREQ,1,2)) 1448 ELSE IF (CSYM(1:6).EQ.'SPHTOP') THEN 1449 GAMMA = 0.3d0 * (HYPERPOL(IFREQ,1,1)+HYPERPOL(IFREQ,1,2)) 1450 DO ICOMP = 2, 4 1451 GAMMA = GAMMA + 0.2d0 * 1452 & (HYPERPOL(IFREQ,ICOMP,1)+HYPERPOL(IFREQ,ICOMP,2)) 1453 END DO 1454 ELSE IF (CSYM(1:6).EQ.'LINEAR') THEN 1455 GAMMA= 1.5d0*(HYPERPOL(IFREQ,1,1)+HYPERPOL(IFREQ,1,2)) 1456 GAMMA=GAMMA+4.0d0*(HYPERPOL(IFREQ,8,1)+HYPERPOL(IFREQ,8,2)) 1457 DO ICOMP = 2, 7 1458 GAMMA = GAMMA + 1459 & (HYPERPOL(IFREQ,ICOMP,1)+HYPERPOL(IFREQ,ICOMP,2)) 1460 END DO 1461 GAMMA = GAMMA / 15.0d0 1462 ELSE 1463 GAMMA = HYPERPOL(IFREQ,1,1) + HYPERPOL(IFREQ,1,2) 1464 GAMMA = GAMMA + HYPERPOL(IFREQ,8,1) + HYPERPOL(IFREQ,8,2) 1465 GAMMA = GAMMA + HYPERPOL(IFREQ,15,1) + HYPERPOL(IFREQ,15,2) 1466 DO ICOMP = 1, 21 1467 GAMMA = GAMMA + 1468 & HALF*(HYPERPOL(IFREQ,ICOMP,1)+HYPERPOL(IFREQ,ICOMP,2)) 1469 END DO 1470 GAMMA = GAMMA / 15.0d0 1471 END IF 1472 1473 IF (IFREQ.EQ.1) THEN 1474 WRITE(LUPRI,'(1X,A8,4X,F7.4,2X,3(3X,F7.4,2X),G16.8)') 1475 & 'gamma_||', ACRFR(IFREQ), BCRFR(IFREQ), 1476 & CCRFR(IFREQ), DCRFR(IFREQ), GAMMA 1477 ELSE 1478 WRITE(LUPRI,'(1X,8X,4X,F7.4,2X,3(3X,F7.4,2X),G16.8)') 1479 & ACRFR(IFREQ), BCRFR(IFREQ), 1480 & CCRFR(IFREQ), DCRFR(IFREQ), GAMMA 1481 END IF 1482 END DO 1483 1484 END IF 1485 1486 IF (GAMMA_ORT) THEN 1487 1488* calculate & print gamma_{_|_} for all requested frequencies: 1489 DO IFREQ = 1, NCRFREQ 1490 1491 IF (CSYM(1:6).EQ.'ATOMIC') THEN 1492 GAMMA = (HYPERPOL(IFREQ,2,1)+HYPERPOL(IFREQ,2,2)) 1493 & - HALF*(HYPERPOL(IFREQ,1,1)+HYPERPOL(IFREQ,1,2)) 1494 & + (HYPERPOL(IFREQ,3,1)+HYPERPOL(IFREQ,3,2)) 1495 ELSE IF (CSYM(1:6).EQ.'SPHTOP') THEN 1496 GAMMA = 0.1d0 * (HYPERPOL(IFREQ,1,1)+HYPERPOL(IFREQ,1,2)) 1497 & + 0.4d0 * (HYPERPOL(IFREQ,2,1)+HYPERPOL(IFREQ,2,2)) 1498 & - 0.1d0 * (HYPERPOL(IFREQ,3,1)+HYPERPOL(IFREQ,3,2)) 1499 & - 0.1d0 * (HYPERPOL(IFREQ,4,1)+HYPERPOL(IFREQ,4,2)) 1500 ELSE IF (CSYM(1:6).EQ.'LINEAR') THEN 1501 GAMMA= 1.0d0 * (HYPERPOL(IFREQ,1,1)+HYPERPOL(IFREQ,1,2)) 1502 & + 4.0d0 * (HYPERPOL(IFREQ,2,1)+HYPERPOL(IFREQ,2,2)) 1503 & - 1.0d0 * (HYPERPOL(IFREQ,3,1)+HYPERPOL(IFREQ,3,2)) 1504 & - 1.0d0 * (HYPERPOL(IFREQ,4,1)+HYPERPOL(IFREQ,4,2)) 1505 & + 4.0d0 * (HYPERPOL(IFREQ,5,1)+HYPERPOL(IFREQ,5,2)) 1506 & - 1.0d0 * (HYPERPOL(IFREQ,6,1)+HYPERPOL(IFREQ,6,2)) 1507 & - 1.0d0 * (HYPERPOL(IFREQ,7,1)+HYPERPOL(IFREQ,7,2)) 1508 & + 1.0d0 * (HYPERPOL(IFREQ,8,1)+HYPERPOL(IFREQ,8,2)) 1509 & + 5.0d0 * (HYPERPOL(IFREQ,9,1)+HYPERPOL(IFREQ,9,2)) 1510 GAMMA = GAMMA / 30.0d0 1511 ELSE 1512 GAMMA = 0.0d0 1513 DO IADD = 0, 14, 7 1514 GAMMA = GAMMA + 1515 & 1.0d0*(HYPERPOL(IFREQ,1+IADD,1)+HYPERPOL(IFREQ,1+IADD,2)) 1516 & +2.0d0*(HYPERPOL(IFREQ,2+IADD,1)+HYPERPOL(IFREQ,2+IADD,2)) 1517 & -0.5d0*(HYPERPOL(IFREQ,3+IADD,1)+HYPERPOL(IFREQ,3+IADD,2)) 1518 & -0.5d0*(HYPERPOL(IFREQ,4+IADD,1)+HYPERPOL(IFREQ,4+IADD,2)) 1519 & +2.0d0*(HYPERPOL(IFREQ,5+IADD,1)+HYPERPOL(IFREQ,5+IADD,2)) 1520 & -0.5d0*(HYPERPOL(IFREQ,6+IADD,1)+HYPERPOL(IFREQ,6+IADD,2)) 1521 & -0.5d0*(HYPERPOL(IFREQ,7+IADD,1)+HYPERPOL(IFREQ,7+IADD,2)) 1522 END DO 1523 GAMMA = GAMMA / 30.0d0 1524 END IF 1525 1526 IF (IFREQ.EQ.1) THEN 1527 WRITE(LUPRI,'(/1X,A8,4X,F7.4,2X,3(3X,F7.4,2X),G16.8)') 1528 & 'gamma_|_', ACRFR(IFREQ), BCRFR(IFREQ), 1529 & CCRFR(IFREQ), DCRFR(IFREQ), GAMMA 1530 ELSE 1531 WRITE(LUPRI,'(1X,8X,4X,F7.4,2X,3(3X,F7.4,2X),G16.8)') 1532 & ACRFR(IFREQ), BCRFR(IFREQ), 1533 & CCRFR(IFREQ), DCRFR(IFREQ), GAMMA 1534 END IF 1535 END DO 1536 1537 END IF 1538 1539 IF (GAMMA_ORT) THEN 1540 1541* calculate & print gamma_{ms} for all requested frequencies: 1542 DO IFREQ = 1, NCRFREQ 1543 1544 IF (CSYM(1:6).EQ.'ATOMIC') THEN 1545 GAMMA = 3.0d0*HALF*(HYPERPOL(IFREQ,2,1)+HYPERPOL(IFREQ,2,2)) 1546 & - 2.0d0*HALF*(HYPERPOL(IFREQ,1,1)+HYPERPOL(IFREQ,1,2)) 1547 & + 3.0d0*HALF*(HYPERPOL(IFREQ,3,1)+HYPERPOL(IFREQ,3,2)) 1548 ELSE IF (CSYM(1:6).EQ.'SPHTOP') THEN 1549 GAMMA = 0.5d0 * (HYPERPOL(IFREQ,2,1)+HYPERPOL(IFREQ,2,2)) 1550 & + 0.5d0 * (HYPERPOL(IFREQ,3,1)+HYPERPOL(IFREQ,3,2)) 1551 & - 1.0d0 * (HYPERPOL(IFREQ,4,1)+HYPERPOL(IFREQ,4,2)) 1552 ELSE IF (CSYM(1:6).EQ.'LINEAR') THEN 1553 GAMMA = 0.5d0 * (HYPERPOL(IFREQ,2,1)+HYPERPOL(IFREQ,2,2)) 1554 & + 0.5d0 * (HYPERPOL(IFREQ,3,1)+HYPERPOL(IFREQ,3,2)) 1555 & - 1.0d0 * (HYPERPOL(IFREQ,4,1)+HYPERPOL(IFREQ,4,2)) 1556 & + 0.5d0 * (HYPERPOL(IFREQ,5,1)+HYPERPOL(IFREQ,5,2)) 1557 & + 0.5d0 * (HYPERPOL(IFREQ,6,1)+HYPERPOL(IFREQ,6,2)) 1558 & - 1.0d0 * (HYPERPOL(IFREQ,7,1)+HYPERPOL(IFREQ,7,2)) 1559 & + 1.5d0 * (HYPERPOL(IFREQ,9,1)+HYPERPOL(IFREQ,9,2)) 1560 & + 1.5d0 * (HYPERPOL(IFREQ,10,1)+HYPERPOL(IFREQ,10,2)) 1561 & - 1.0d0 * (HYPERPOL(IFREQ,8,1)+HYPERPOL(IFREQ,8,2)) 1562 GAMMA = GAMMA / 3.0d0 1563 ELSE 1564 GAMMA = 0.0d0 1565 DO IADD = 0, 14, 7 1566 GAMMA = GAMMA + 1567 & 0.5d0*(HYPERPOL(IFREQ,2+IADD,1)+HYPERPOL(IFREQ,2+IADD,2)) 1568 & +0.5d0*(HYPERPOL(IFREQ,3+IADD,1)+HYPERPOL(IFREQ,3+IADD,2)) 1569 & -1.0d0*(HYPERPOL(IFREQ,4+IADD,1)+HYPERPOL(IFREQ,4+IADD,2)) 1570 & +0.5d0*(HYPERPOL(IFREQ,5+IADD,1)+HYPERPOL(IFREQ,5+IADD,2)) 1571 & +0.5d0*(HYPERPOL(IFREQ,6+IADD,1)+HYPERPOL(IFREQ,6+IADD,2)) 1572 & -1.0d0*(HYPERPOL(IFREQ,7+IADD,1)+HYPERPOL(IFREQ,7+IADD,2)) 1573 END DO 1574 GAMMA = GAMMA / 6.0d0 1575 END IF 1576 1577 IF (IFREQ.EQ.1) THEN 1578 WRITE(LUPRI,'(/1X,A8,4X,F7.4,2X,3(3X,F7.4,2X),F16.7)') 1579 & 'gamma_ms', ACRFR(IFREQ), BCRFR(IFREQ), 1580 & CCRFR(IFREQ), DCRFR(IFREQ), GAMMA 1581 ELSE 1582 WRITE(LUPRI,'(1X,8X,4X,F7.4,2X,3(3X,F7.4,2X),G16.8)') 1583 & ACRFR(IFREQ), BCRFR(IFREQ), 1584 & CCRFR(IFREQ), DCRFR(IFREQ), GAMMA 1585 END IF 1586 END DO 1587 1588 END IF 1589 1590 END IF 1591 1592 WRITE(LUPRI,'(75("-"))') 1593 1594 CALL FLSHFO(LUPRI) 1595 1596 RETURN 1597 END 1598*---------------------------------------------------------------------* 1599* END OF SUBROUTINE CCCHYPPRT * 1600*---------------------------------------------------------------------* 1601 SUBROUTINE CCCEXPPRT(EXPCOF) 1602*---------------------------------------------------------------------* 1603* 1604* Purpose: print expansion coefficients for 1605* second hyperpolarizabilities 1606* 1607* 1608* Written by Christof Haettig in March 1998. 1609* 1610*=====================================================================* 1611#if defined (IMPLICIT_NONE) 1612 IMPLICIT NONE 1613#else 1614# include "implicit.h" 1615#endif 1616#include "priunit.h" 1617#include "ccorb.h" 1618#include "ccsdinp.h" 1619#include "cccrinf.h" 1620#include "ccroper.h" 1621 1622 1623 CHARACTER*10 BLANKS 1624 CHARACTER*80 STRING 1625 INTEGER ISYMA, ISYMB, ISYMC, ISYMD 1626 INTEGER IDISP, IOPER, ISYOLD, ICOMP, IADD 1627 1628 1629#if defined (SYS_CRAY) 1630 REAL EXPCOF(NCRDISP,NCROPER) 1631 REAL GAMMA 1632#else 1633 DOUBLE PRECISION EXPCOF(NCRDISP,NCROPER) 1634 DOUBLE PRECISION GAMMA 1635#endif 1636 1637*---------------------------------------------------------------------* 1638* print header for hyperpolarizability section 1639*---------------------------------------------------------------------* 1640 BLANKS = ' ' 1641 STRING = ' RESULTS FOR EXPANSION COEFFICIENTS' 1642 1643 IF (CCS) THEN 1644 CALL AROUND( BLANKS//'FINAL CCS'//STRING(1:45)//BLANKS ) 1645 ELSE IF (CC2) THEN 1646 CALL AROUND( BLANKS//'FINAL CC2'//STRING(1:45)//BLANKS ) 1647 ELSE IF (CCSD) THEN 1648 CALL AROUND( BLANKS//'FINAL CCSD'//STRING(1:45)//BLANKS ) 1649 ELSE IF (CC3) THEN 1650 CALL AROUND( BLANKS//'FINAL CC3'//STRING(1:45)//BLANKS ) 1651 ELSE 1652 CALL QUIT('CCCEXPPRT called for an unknown Coupled '// 1653 & 'Cluster model.') 1654 END IF 1655 1656 WRITE(LUPRI,'(/1X,4(1X,A," operator",3X),4X,A,/,72("-"))') 1657 & 'A','B','C','D','d_ABCD' 1658 1659 ISYOLD = 1 1660 1661 DO IOPER = 1, NCROPER 1662 ISYMA = ISYOPR(IACROP(IOPER)) 1663 ISYMB = ISYOPR(IBCROP(IOPER)) 1664 ISYMC = ISYOPR(ICCROP(IOPER)) 1665 ISYMD = ISYOPR(IDCROP(IOPER)) 1666 1667 1668 IDISP = 1 1669 IF (MULD2H(ISYMA,ISYMB).EQ.MULD2H(ISYMC,ISYMD)) THEN 1670 WRITE(LUPRI,'(/1X,4(A8,I3,3X),G16.8)') 1671 & LBLOPR(IACROP(IOPER)),ICCAUA(IDISP), 1672 & LBLOPR(IBCROP(IOPER)),ICCAUB(IDISP), 1673 & LBLOPR(ICCROP(IOPER)),ICCAUC(IDISP), 1674 & LBLOPR(IDCROP(IOPER)),ICCAUD(IDISP), 1675 & -EXPCOF(IDISP,IOPER) 1676 ELSE 1677 IF (IPRCHYP.GT.0) THEN 1678 WRITE(LUPRI,'(1X,4(A8,A3,3X),6X,A,7X)') 1679 & LBLOPR(IACROP(IOPER)),' - ', 1680 & LBLOPR(IBCROP(IOPER)),' - ', 1681 & LBLOPR(ICCROP(IOPER)),' - ', 1682 & LBLOPR(IDCROP(IOPER)),' - ', 1683 & '---' 1684 END IF 1685 END IF 1686 1687 DO IDISP = 2, NCRDISP 1688 IF (MULD2H(ISYMA,ISYMB).EQ.MULD2H(ISYMC,ISYMD)) THEN 1689 WRITE(LUPRI,'(1X,4(8X,I3,3X),G16.8)') 1690 & ICCAUA(IDISP),ICCAUB(IDISP),ICCAUC(IDISP),ICCAUD(IDISP), 1691 & -EXPCOF(IDISP,IOPER) 1692 END IF 1693 END DO 1694 1695 END DO 1696 1697 IF (GAMMA_PAR .OR. GAMMA_ORT) THEN 1698 WRITE(LUPRI,'(/////A,/,86("-"))')' average ' 1699 END IF 1700 1701 IF (GAMMA_PAR) THEN 1702* calculate & print gamma_{||} for all requested orders: 1703 DO IDISP = 1, NCRDISP 1704 1705 IF (CSYM(1:6).EQ.'ATOMIC') THEN 1706 GAMMA = EXPCOF(IDISP,1) 1707 ELSE IF (CSYM(1:6).EQ.'SPHTOP') THEN 1708 GAMMA = 0.6d0 * EXPCOF(IDISP,1) 1709 DO ICOMP = 2, 4 1710 GAMMA = GAMMA + 0.4d0 * EXPCOF(IDISP,ICOMP) 1711 END DO 1712 ELSE IF (CSYM(1:6).EQ.'LINEAR') THEN 1713 GAMMA = 3.0d0 * EXPCOF(IDISP,1) 1714 GAMMA = GAMMA + 8.0d0 * EXPCOF(IDISP,8) 1715 DO ICOMP = 2, 7 1716 GAMMA = GAMMA + 2.0d0 * EXPCOF(IDISP,ICOMP) 1717 END DO 1718 GAMMA = GAMMA / 15.0d0 1719 ELSE 1720 GAMMA = 2.0d0 * EXPCOF(IDISP,1) 1721 GAMMA = GAMMA + 2.0d0 * EXPCOF(IDISP,8) 1722 GAMMA = GAMMA + 2.0d0 * EXPCOF(IDISP,15) 1723 DO ICOMP = 1, 21 1724 GAMMA = GAMMA + EXPCOF(IDISP,ICOMP) 1725 END DO 1726 GAMMA = GAMMA / 15.0d0 1727 END IF 1728 1729 IF (IDISP.EQ.1) THEN 1730 WRITE(LUPRI,'(1X,A8,4X,I3,2X,3(8X,I3,2X),G16.8)') 1731 & 'gamma_||', ICCAUA(IDISP), ICCAUB(IDISP), 1732 & ICCAUC(IDISP), ICCAUD(IDISP), GAMMA 1733 ELSE 1734 WRITE(LUPRI,'(1X,8X,4X,I3,2X,3(8X,I3,2X),G16.8)') 1735 & ICCAUA(IDISP), ICCAUB(IDISP), 1736 & ICCAUC(IDISP), ICCAUD(IDISP), GAMMA 1737 END IF 1738 END DO 1739 END IF ! (GAMMA_PAR) 1740 1741 IF (GAMMA_ORT) THEN 1742* calculate & print gamma_{_|_} for all requested orders: 1743 DO IDISP = 1, NCRDISP 1744 1745 IF (CSYM(1:6).EQ.'ATOMIC') THEN 1746 GAMMA = 2.0d0*EXPCOF(IDISP,2) 1747 & - 1.0d0*EXPCOF(IDISP,1) 1748 & + 2.0d0*EXPCOF(IDISP,3) 1749 ELSE IF (CSYM(1:6).EQ.'SPHTOP') THEN 1750 GAMMA = 0.2d0 * EXPCOF(IDISP,1) 1751 & + 0.8d0 * EXPCOF(IDISP,2) 1752 & - 0.2d0 * EXPCOF(IDISP,3) 1753 & - 0.2d0 * EXPCOF(IDISP,4) 1754 ELSE IF (CSYM(1:6).EQ.'LINEAR') THEN 1755 GAMMA= 1.0d0 * EXPCOF(IDISP,1) 1756 & + 4.0d0 * EXPCOF(IDISP,2) 1757 & - 1.0d0 * EXPCOF(IDISP,3) 1758 & - 1.0d0 * EXPCOF(IDISP,4) 1759 & + 4.0d0 * EXPCOF(IDISP,5) 1760 & - 1.0d0 * EXPCOF(IDISP,6) 1761 & - 1.0d0 * EXPCOF(IDISP,7) 1762 & + 1.0d0 * EXPCOF(IDISP,8) 1763 & + 5.0d0 * EXPCOF(IDISP,9) 1764 GAMMA = GAMMA / 15.0d0 1765 ELSE 1766 GAMMA = 0.0d0 1767 DO IADD = 0, 14, 7 1768 GAMMA = GAMMA + 1769 & 1.0d0 * EXPCOF(IDISP,1+IADD) 1770 & +2.0d0 * EXPCOF(IDISP,2+IADD) 1771 & -0.5d0 * EXPCOF(IDISP,3+IADD) 1772 & -0.5d0 * EXPCOF(IDISP,4+IADD) 1773 & +2.0d0 * EXPCOF(IDISP,5+IADD) 1774 & -0.5d0 * EXPCOF(IDISP,6+IADD) 1775 & -0.5d0 * EXPCOF(IDISP,7+IADD) 1776 END DO 1777 GAMMA = GAMMA / 15.0d0 1778 END IF 1779 1780 IF (IDISP.EQ.1) THEN 1781 WRITE(LUPRI,'(/1X,A8,4X,I3,2X,3(8X,I3,2X),G16.8)') 1782 & 'gamma_|_', ICCAUA(IDISP), ICCAUB(IDISP), 1783 & ICCAUC(IDISP), ICCAUD(IDISP), GAMMA 1784 ELSE 1785 WRITE(LUPRI,'(1X,8X,4X,I3,2X,3(8X,I3,2X),G16.8)') 1786 & ICCAUA(IDISP), ICCAUB(IDISP), 1787 & ICCAUC(IDISP), ICCAUD(IDISP), GAMMA 1788 END IF 1789 END DO 1790 END IF ! (GAMMA_ORT) 1791 1792 IF (GAMMA_ORT) THEN 1793* calculate & print gamma_{ms} for all requested orders: 1794 DO IDISP = 1, NCRDISP 1795 1796 IF (CSYM(1:6).EQ.'ATOMIC') THEN 1797 GAMMA = 3.0d0*EXPCOF(IDISP,2) 1798 & - 2.0d0*EXPCOF(IDISP,1) 1799 & + 3.0d0*EXPCOF(IDISP,3) 1800 ELSE IF (CSYM(1:6).EQ.'SPHTOP') THEN 1801 GAMMA = 1.0d0 * EXPCOF(IDISP,2) 1802 & + 1.0d0 * EXPCOF(IDISP,3) 1803 & - 2.0d0 * EXPCOF(IDISP,4) 1804 ELSE IF (CSYM(1:6).EQ.'LINEAR') THEN 1805 GAMMA = 0.5d0 * EXPCOF(IDISP,2) 1806 & + 0.5d0 * EXPCOF(IDISP,3) 1807 & - 1.0d0 * EXPCOF(IDISP,4) 1808 & + 0.5d0 * EXPCOF(IDISP,5) 1809 & + 0.5d0 * EXPCOF(IDISP,6) 1810 & - 1.0d0 * EXPCOF(IDISP,7) 1811 & + 1.5d0 * EXPCOF(IDISP,9) 1812 & + 1.5d0 * EXPCOF(IDISP,10) 1813 & - 1.0d0 * EXPCOF(IDISP,8) 1814 GAMMA = GAMMA / 1.5d0 1815 ELSE 1816 GAMMA = 0.0d0 1817 DO IADD = 0, 14, 7 1818 GAMMA = GAMMA + 1819 & 0.5d0 * EXPCOF(IDISP,2+IADD) 1820 & + 0.5d0 * EXPCOF(IDISP,3+IADD) 1821 & - 1.0d0 * EXPCOF(IDISP,4+IADD) 1822 & + 0.5d0 * EXPCOF(IDISP,5+IADD) 1823 & + 0.5d0 * EXPCOF(IDISP,6+IADD) 1824 & - 1.0d0 * EXPCOF(IDISP,7+IADD) 1825 END DO 1826 GAMMA = GAMMA / 3.0d0 1827 END IF 1828 1829 IF (IDISP.EQ.1) THEN 1830 WRITE(LUPRI,'(/1X,A8,4X,I3,2X,3(8X,I3,2X),G16.8)') 1831 & 'gamma_ms', ICCAUA(IDISP), ICCAUB(IDISP), 1832 & ICCAUC(IDISP), ICCAUD(IDISP), GAMMA 1833 ELSE 1834 WRITE(LUPRI,'(1X,8X,4X,I3,2X,3(8X,I3,2X),G16.8)') 1835 & ICCAUA(IDISP), ICCAUB(IDISP), 1836 & ICCAUC(IDISP), ICCAUD(IDISP), GAMMA 1837 END IF 1838 END DO 1839 END IF ! (GAMMA_ORT) 1840 1841 1842 CALL FLSHFO(LUPRI) 1843 1844 RETURN 1845 END 1846*---------------------------------------------------------------------* 1847* END OF SUBROUTINE CCCEXPPRT * 1848*---------------------------------------------------------------------* 1849c /* deck CCCRDISP */ 1850*=====================================================================* 1851 SUBROUTINE CCCRDISP(DISPCF,AVEDISP,ABCDE, 1852 & THGCF, ESHGCF, KERRCF,DFWMCF, 1853 & AVETHG,AVESHG, AVEKERR, AVEDFWM, 1854 & EXPCOF,TRINOM, NDSPCOF,LERROR) 1855*---------------------------------------------------------------------* 1856* 1857* Purpose: calculate from the expansion coefficients defined by 1858* gamma_ABCD=sum_{klmn} w_A^k w_B^l w_C^m w_D^n d_ABCD(klmn) 1859* the physical more relevant dispersion coefficients defined 1860* by gamma_ABCD = sum_{lmn} w_B^l w_C^m w_D^l D_ABCD(lmn), 1861* special versions thereof for THG, ESHG, Kerr & DFWM, 1862* as well as coefficients for the isotropic averages and 1863* the `universal' dispersion coefficients A, B, ... 1864* 1865* 1866* Written by Christof Haettig in March 1998. 1867* 1868*=====================================================================* 1869#if defined (IMPLICIT_NONE) 1870 IMPLICIT NONE 1871#else 1872# include "implicit.h" 1873#endif 1874#include "priunit.h" 1875#include "cccrinf.h" 1876 1877 LOGICAL LOCDBG 1878 PARAMETER (LOCDBG = .FALSE.) 1879 1880 INTEGER KA, KB, KBP, KC, KCP, KCPP, KD, KDP, KDPP, KDPPP 1881 INTEGER KE, KEP, KEPP, KEPPP, KEPPPP 1882 PARAMETER (KA=1, KB=2, KBP=3, KC=4, KCP=5, KCPP=6) 1883 PARAMETER (KD=7, KDP=8, KDPP=9, KDPPP=10) 1884 PARAMETER (KE=11, KEP=12, KEPP=13, KEPPP=14, KEPPPP=15) 1885 1886 INTEGER NDSPCOF ! leading dimension of DISPCF, AVEDISP, TRINOM 1887 LOGICAL LERROR 1888 1889#if defined (SYS_CRAY) 1890 REAL DISPCF(NDSPCOF,NCROPER) 1891 REAL AVEDISP(NDSPCOF,4) 1892 REAL ABCDE(15,2) 1893 REAL THGCF(NCRDSPE+1,NCROPER) 1894 REAL ESHGCF(NCRDSPE+1,NCROPER) 1895 REAL KERRCF(NCRDSPE+1,NCROPER) 1896 REAL DFWMCF(NCRDSPE+1,NCROPER) 1897 REAL AVESHG(NCRDSPE+1,4) 1898 REAL AVETHG(NCRDSPE+1,4) 1899 REAL AVEKERR(NCRDSPE+1,4) 1900 REAL AVEDFWM(NCRDSPE+1,4) 1901 REAL EXPCOF(NCRDISP,NCROPER) 1902 REAL TRINOM(NDSPCOF) 1903 REAL GAMMA0,A,B,BP,C,CP,CPP 1904 REAL D, DP, DPP, DPPP, E, EP, EPP, EPPP, EPPPP 1905 REAL ERROR, ERR1, ERR2, ERR3, ERR4, ERR5, ERR6, ERR7 1906 REAL ERR0, ERR8, ERR9, ERR10, ERR11, ERR12, ERR13, ERR14 1907 REAL ERR15, ERR16 1908 REAL Y, X1, X2, X3, X4, Z, W 1909 REAL THRSDISP 1910#else 1911 DOUBLE PRECISION DISPCF(NDSPCOF,NCROPER) 1912 DOUBLE PRECISION AVEDISP(NDSPCOF,4) 1913 DOUBLE PRECISION ABCDE(15,2) 1914 DOUBLE PRECISION THGCF(NCRDSPE+1,NCROPER) 1915 DOUBLE PRECISION ESHGCF(NCRDSPE+1,NCROPER) 1916 DOUBLE PRECISION KERRCF(NCRDSPE+1,NCROPER) 1917 DOUBLE PRECISION DFWMCF(NCRDSPE+1,NCROPER) 1918 DOUBLE PRECISION AVESHG(NCRDSPE+1,4) 1919 DOUBLE PRECISION AVETHG(NCRDSPE+1,4) 1920 DOUBLE PRECISION AVEKERR(NCRDSPE+1,4) 1921 DOUBLE PRECISION AVEDFWM(NCRDSPE+1,4) 1922 DOUBLE PRECISION EXPCOF(NCRDISP,NCROPER) 1923 DOUBLE PRECISION TRINOM(NDSPCOF) 1924 DOUBLE PRECISION GAMMA0,A,B,BP,C,CP,CPP 1925 DOUBLE PRECISION D, DP, DPP, DPPP, E, EP, EPP, EPPP, EPPPP 1926 DOUBLE PRECISION ERROR, ERR1, ERR2, ERR3, ERR4, ERR5, ERR6, ERR7 1927 DOUBLE PRECISION ERR0,ERR8,ERR9,ERR10,ERR11,ERR12,ERR13,ERR14 1928 DOUBLE PRECISION ERR15, ERR16 1929 DOUBLE PRECISION Y, X1, X2, X3, X4, Z, W 1930 DOUBLE PRECISION THRSDISP 1931#endif 1932 PARAMETER (THRSDISP = 1.0d-5) 1933 1934 INTEGER IDISP, ICOMP, IOPER,OFFEX, IEXPCF 1935 INTEGER L, N, M, LMN, MN, P, Q, R, I, J, IADD 1936 1937 INTEGER ITRI 1938#if defined (SYS_CRAY) 1939 REAL SIGNEO 1940#else 1941 DOUBLE PRECISION SIGNEO 1942#endif 1943 1944* index for 3-dim array with L >= M >= N and L,M,N >= 0 : 1945 ITRI(L,M,N) = L*(L+1)*(L+2)/6 + M*(M+1)/2 + N + 1 1946 1947* this gives (-1)^I: 1948 SIGNEO(I) = DBLE(1-MOD(I,2)*2) 1949 1950C 1951*---------------------------------------------------------------------* 1952* initialze flag for numerical errors in A,B... coefficients: 1953*---------------------------------------------------------------------* 1954 LERROR = .FALSE. 1955 1956*---------------------------------------------------------------------* 1957* precompute trinomial coefficients: L! / ( (L-M)! (M-N)! N! ) 1958*---------------------------------------------------------------------* 1959 DO L = 0, NCRDSPE 1960 TRINOM(ITRI(L,0,0)) = 1.0d0 1961 TRINOM(ITRI(L,L,0)) = 1.0d0 1962 TRINOM(ITRI(L,L,L)) = 1.0d0 1963 DO M = 1, L-1 1964 TRINOM(ITRI(L,M,0)) = 1965 & TRINOM(ITRI(L-1,M,0)) + TRINOM(ITRI(L-1,M-1,0)) 1966 TRINOM(ITRI(L,M,M)) = 1967 & TRINOM(ITRI(L-1,M,M)) + TRINOM(ITRI(L-1,M-1,M-1)) 1968 DO N = 1, M-1 1969 TRINOM(ITRI(L,M,N)) = TRINOM(ITRI(L-1,M,N)) + 1970 & TRINOM(ITRI(L-1,M-1,N)) + TRINOM(ITRI(L-1,M-1,N-1)) 1971 END DO 1972 END DO 1973 DO N = 1, L-1 1974 TRINOM(ITRI(L,L,N)) = 1975 & TRINOM(ITRI(L-1,L-1,N)) + TRINOM(ITRI(L-1,L-1,N-1)) 1976 END DO 1977 END DO 1978 1979 IF (LOCDBG) THEN 1980 WRITE (LUPRI,'(///)') 1981 WRITE (LUPRI,*) 'DEBUG_CCCRDISP> trinomial coefficients:' 1982 DO L = 0, NCRDSPE 1983 DO M = 0, L 1984 WRITE (LUPRI,'(2I3,4X,15F8.2)') 1985 & L,M,(TRINOM(ITRI(L,M,N)),N=0,M) 1986 END DO 1987 END DO 1988 WRITE (LUPRI,*) 1989 END IF 1990 1991 DO IOPER = 1, NCROPER 1992 OFFEX = 0 1993 DO LMN = 0, NCRDSPE, 2 1994 DO MN = 0, LMN 1995 DO N = 0, MN 1996 L = LMN - MN 1997 M = MN - N 1998 1999 IF (LOCDBG) THEN 2000 WRITE (LUPRI,'(//A,A)') 2001 & ' L M N P Q R SIGN TRINOM', 2002 & ' I J K L EXPCOF DISPCF' 2003 END IF 2004 2005 DISPCF(ITRI(LMN,MN,N),IOPER) = 0.0d0 2006 2007 DO P = 0, L 2008 DO Q = 0, M 2009 DO R = 0, N 2010 IEXPCF = OFFEX + ITRI(LMN-P-Q-R,MN-Q-R,N-R) 2011 2012 DISPCF(ITRI(LMN,MN,N),IOPER) = 2013 & DISPCF(ITRI(LMN,MN,N),IOPER) + SIGNEO(P+Q+R) * 2014 & TRINOM(ITRI(P+Q+R,Q+R,R)) * EXPCOF(IEXPCF,IOPER) 2015 2016 IF (LOCDBG) THEN 2017 WRITE (LUPRI,'(6I3,2F8.2,3X,4I2,2E16.8,5I5)') 2018 & L,M,N,P,Q,R,SIGNEO(P+Q+R),TRINOM(ITRI(P+Q+R,Q+R,R)), 2019 & P+Q+R,L-P,M-Q,N-R, EXPCOF(IEXPCF,IOPER), 2020 & DISPCF(ITRI(LMN,MN,N),IOPER), 2021 & ITRI(P+Q+R,Q+R,R),ITRI(LMN,MN,N),IEXPCF,IOPER 2022 END IF 2023 2024 END DO ! R 2025 END DO ! Q 2026 END DO ! P 2027 2028 END DO ! N 2029 END DO ! MN 2030 2031 OFFEX = OFFEX + (LMN+3)*(LMN+2)*(LMN+1)/6 2032 2033 END DO ! LMN 2034 END DO ! IOPER 2035 2036 2037 DO IOPER = 1, NCROPER 2038 DO LMN = 0, NCRDSPE, 2 2039 IF (LOCDBG) THEN 2040 WRITE (LUPRI,'(//A,5(5X,A,5X))') ' L M N ', 2041 & 'DISPCF', 'THGCF ', 'ESHGCF', 'KERRCF', 'DFWMCF' 2042 END IF 2043 THGCF(LMN+1,IOPER) = 0.0d0 2044 ESHGCF(LMN+1,IOPER) = 0.0d0 2045 KERRCF(LMN+1,IOPER) = 0.0d0 2046 DFWMCF(LMN+1,IOPER) = 0.0d0 2047 2048 KERRCF(LMN+1,IOPER) = KERRCF(LMN+1,IOPER) + 2049 & DISPCF(ITRI(LMN,LMN,LMN),IOPER) 2050 DO MN = 0, LMN 2051 ESHGCF(LMN+1,IOPER) = ESHGCF(LMN+1,IOPER) + 2052 & DISPCF(ITRI(LMN,MN,0),IOPER) 2053 DO N = 0, MN 2054 THGCF(LMN+1,IOPER) = THGCF(LMN+1,IOPER) + 2055 & DISPCF(ITRI(LMN,MN,N),IOPER) 2056 2057 DFWMCF(LMN+1,IOPER) = DFWMCF(LMN+1,IOPER) + 2058 & SIGNEO(N) * DISPCF(ITRI(LMN,MN,N),IOPER) 2059 2060 IF (LOCDBG) THEN 2061 WRITE (LUPRI,'(3I3,5E16.8)') LMN-MN,MN-N,N, 2062 & DISPCF(ITRI(LMN,MN,N),IOPER), 2063 & THGCF(LMN+1,IOPER), ESHGCF(LMN+1,IOPER), 2064 & KERRCF(LMN+1,IOPER),DFWMCF(LMN+1,IOPER) 2065 END IF 2066 2067 END DO 2068 END DO 2069 2070 END DO 2071 END DO 2072 2073 IF (GAMMA_PAR .OR. GAMMA_ORT) THEN 2074 2075* calculate averaged vector components parallel/orthogonal to z axis: 2076 DO IDISP = 1, (NCRDSPE+3)*(NCRDSPE+2)*(NCRDSPE+1)/6 2077 IF (CSYM(1:6).EQ.'ATOMIC') THEN 2078 IF (GAMMA_PAR) THEN 2079 AVEDISP(IDISP,1) = 1.0d0*DISPCF(IDISP,1) 2080 END IF 2081 IF (GAMMA_ORT) THEN 2082 AVEDISP(IDISP,2) = 2.0d0*DISPCF(IDISP,2) 2083 & - 1.0d0*DISPCF(IDISP,1) 2084 & + 2.0d0*DISPCF(IDISP,3) 2085 AVEDISP(IDISP,3) = 3.0d0*DISPCF(IDISP,2) 2086 & - 2.0d0*DISPCF(IDISP,1) 2087 & + 3.0d0*DISPCF(IDISP,3) 2088 END IF 2089 ELSE IF (CSYM(1:6).EQ.'SPHTOP') THEN 2090 IF (GAMMA_PAR) THEN 2091 AVEDISP(IDISP,1) = 0.6d0 * DISPCF(IDISP,1) 2092 & + 0.4d0 * DISPCF(IDISP,2) 2093 & + 0.4d0 * DISPCF(IDISP,3) 2094 & + 0.4d0 * DISPCF(IDISP,4) 2095 END IF 2096 IF (GAMMA_ORT) THEN 2097 AVEDISP(IDISP,2) = 0.2d0 * DISPCF(IDISP,1) 2098 & + 0.8d0 * DISPCF(IDISP,2) 2099 & + 0.8d0 * DISPCF(IDISP,3) 2100 & - 1.2d0 * DISPCF(IDISP,4) 2101 AVEDISP(IDISP,3) = 1.0d0 * DISPCF(IDISP,2) 2102 & + 1.0d0 * DISPCF(IDISP,3) 2103 & - 2.0d0 * DISPCF(IDISP,4) 2104 END IF 2105 ELSE IF (CSYM(1:6).EQ.'LINEAR') THEN 2106 IF (GAMMA_PAR) THEN 2107 AVEDISP(IDISP,1) = 3.0d0 * DISPCF(IDISP,1) 2108 AVEDISP(IDISP,1) = AVEDISP(IDISP,1) + 2109 & 8.0d0 * DISPCF(IDISP,8) 2110 DO ICOMP = 2, 7 2111 AVEDISP(IDISP,1) = AVEDISP(IDISP,1) + 2112 & 2.0d0 * DISPCF(IDISP,ICOMP) 2113 END DO 2114 AVEDISP(IDISP,1) = AVEDISP(IDISP,1) / 15.0d0 2115 END IF 2116 IF (GAMMA_ORT) THEN 2117 AVEDISP(IDISP,2) = 1.0d0 * DISPCF(IDISP,1) 2118 & + 4.0d0 * DISPCF(IDISP,2) 2119 & + 4.0d0 * DISPCF(IDISP,3) 2120 & - 6.0d0 * DISPCF(IDISP,4) 2121 & + 4.0d0 * DISPCF(IDISP,5) 2122 & + 4.0d0 * DISPCF(IDISP,6) 2123 & - 6.0d0 * DISPCF(IDISP,7) 2124 & + 10.0d0 * DISPCF(IDISP,9) 2125 & + 10.0d0 * DISPCF(IDISP,10) 2126 & - 4.0d0 * DISPCF(IDISP,8) 2127 AVEDISP(IDISP,2) = AVEDISP(IDISP,2) / 15.0d0 2128 AVEDISP(IDISP,3) = 0.5d0 * DISPCF(IDISP,2) 2129 & + 0.5d0 * DISPCF(IDISP,3) 2130 & - 1.0d0 * DISPCF(IDISP,4) 2131 & + 0.5d0 * DISPCF(IDISP,5) 2132 & + 0.5d0 * DISPCF(IDISP,6) 2133 & - 1.0d0 * DISPCF(IDISP,7) 2134 & + 1.5d0 * DISPCF(IDISP,9) 2135 & + 1.5d0 * DISPCF(IDISP,10) 2136 & - 1.0d0 * DISPCF(IDISP,8) 2137 AVEDISP(IDISP,3) = AVEDISP(IDISP,3) / 1.5d0 2138 END IF 2139 ELSE 2140 IF (GAMMA_PAR) THEN 2141 AVEDISP(IDISP,1) = 2.0d0*DISPCF(IDISP,1) 2142 AVEDISP(IDISP,1) = AVEDISP(IDISP,1)+2.0d0*DISPCF(IDISP,8) 2143 AVEDISP(IDISP,1) =AVEDISP(IDISP,1)+2.0d0*DISPCF(IDISP,15) 2144 DO ICOMP = 1, 21 2145 AVEDISP(IDISP,1)=AVEDISP(IDISP,1) + DISPCF(IDISP,ICOMP) 2146 END DO 2147 AVEDISP(IDISP,1) = AVEDISP(IDISP,1) / 15.0d0 2148 END IF 2149 IF (GAMMA_ORT) THEN 2150 AVEDISP(IDISP,2) = 0.0d0 2151 DO IADD = 0, 14, 7 2152 AVEDISP(IDISP,2) = AVEDISP(IDISP,2) + 2153 & 1.0d0 * DISPCF(IDISP,1+IADD) 2154 & +2.0d0 * DISPCF(IDISP,2+IADD) 2155 & +2.0d0 * DISPCF(IDISP,3+IADD) 2156 & -3.0d0 * DISPCF(IDISP,4+IADD) 2157 & +2.0d0 * DISPCF(IDISP,5+IADD) 2158 & -2.0d0 * DISPCF(IDISP,6+IADD) 2159 & -3.0d0 * DISPCF(IDISP,7+IADD) 2160 END DO 2161 AVEDISP(IDISP,2) = AVEDISP(IDISP,2) / 15.0d0 2162 AVEDISP(IDISP,3) = 0.0d0 2163 DO IADD = 0, 14, 7 2164 AVEDISP(IDISP,3) = AVEDISP(IDISP,3) + 2165 & 0.5d0 * DISPCF(IDISP,2+IADD) 2166 & + 0.5d0 * DISPCF(IDISP,3+IADD) 2167 & - 1.0d0 * DISPCF(IDISP,4+IADD) 2168 & + 0.5d0 * DISPCF(IDISP,5+IADD) 2169 & + 0.5d0 * DISPCF(IDISP,6+IADD) 2170 & - 1.0d0 * DISPCF(IDISP,7+IADD) 2171 END DO 2172 AVEDISP(IDISP,3) = AVEDISP(IDISP,3) / 3.0d0 2173 END IF 2174 END IF 2175 END DO 2176 2177 DO IDISP = 1, NCRDSPE+1 2178 IF (CSYM(1:6).EQ.'ATOMIC') THEN 2179 IF (GAMMA_PAR) THEN 2180 AVETHG(IDISP,1) = THGCF(IDISP,1) 2181 AVESHG(IDISP,1) = ESHGCF(IDISP,1) 2182 AVEKERR(IDISP,1) = KERRCF(IDISP,1) 2183 AVEDFWM(IDISP,1) = DFWMCF(IDISP,1) 2184 END IF 2185 IF (GAMMA_ORT) THEN 2186 AVETHG(IDISP,2) = 2.0d0*THGCF(IDISP,2) 2187 & - 1.0d0*THGCF(IDISP,1) 2188 & + 2.0d0*THGCF(IDISP,3) 2189 AVETHG(IDISP,3) = 3.0d0*THGCF(IDISP,2) 2190 & - 2.0d0*THGCF(IDISP,1) 2191 & + 3.0d0*THGCF(IDISP,3) 2192 AVESHG(IDISP,2) = 2.0d0*ESHGCF(IDISP,2) 2193 & - 1.0d0*ESHGCF(IDISP,1) 2194 & + 2.0d0*ESHGCF(IDISP,3) 2195 AVESHG(IDISP,3) = 3.0d0*ESHGCF(IDISP,2) 2196 & - 2.0d0*ESHGCF(IDISP,1) 2197 & + 3.0d0*ESHGCF(IDISP,3) 2198 AVEKERR(IDISP,2) = 2.0d0*KERRCF(IDISP,2) 2199 & - 1.0d0*KERRCF(IDISP,1) 2200 & + 2.0d0*KERRCF(IDISP,3) 2201 AVEKERR(IDISP,3) = 3.0d0*KERRCF(IDISP,2) 2202 & - 2.0d0*KERRCF(IDISP,1) 2203 & + 3.0d0*KERRCF(IDISP,3) 2204 AVEDFWM(IDISP,2) = 2.0d0*DFWMCF(IDISP,2) 2205 & - 1.0d0*DFWMCF(IDISP,1) 2206 & + 2.0d0*DFWMCF(IDISP,3) 2207 AVEDFWM(IDISP,3) = 3.0d0*DFWMCF(IDISP,2) 2208 & - 2.0d0*DFWMCF(IDISP,1) 2209 & + 3.0d0*DFWMCF(IDISP,3) 2210 END IF 2211 ELSE IF (CSYM(1:6).EQ.'SPHTOP') THEN 2212 IF (GAMMA_PAR) THEN 2213 AVETHG(IDISP,1) = 0.6d0 * THGCF(IDISP,1) 2214 AVESHG(IDISP,1) = 0.6d0 * ESHGCF(IDISP,1) 2215 AVEKERR(IDISP,1) = 0.6d0 * KERRCF(IDISP,1) 2216 AVEDFWM(IDISP,1) = 0.6d0 * DFWMCF(IDISP,1) 2217 DO ICOMP = 2, 4 2218 AVETHG(IDISP,1) = AVETHG(IDISP,1) + 2219 & 0.4d0 * THGCF(IDISP,ICOMP) 2220 AVESHG(IDISP,1) = AVESHG(IDISP,1) + 2221 & 0.4d0 * ESHGCF(IDISP,ICOMP) 2222 AVEKERR(IDISP,1) = AVEKERR(IDISP,1) + 2223 & 0.4d0 * KERRCF(IDISP,ICOMP) 2224 AVEDFWM(IDISP,1) = AVEDFWM(IDISP,1) + 2225 & 0.4d0 * DFWMCF(IDISP,ICOMP) 2226 END DO 2227 END IF 2228 IF (GAMMA_ORT) THEN 2229 AVETHG(IDISP,2) = 0.2d0 * THGCF(IDISP,1) 2230 & + 0.8d0 * THGCF(IDISP,2) 2231 & - 0.2d0 * THGCF(IDISP,3) 2232 & - 0.2d0 * THGCF(IDISP,4) 2233 AVETHG(IDISP,3) = 1.0d0 * THGCF(IDISP,2) 2234 & + 1.0d0 * THGCF(IDISP,3) 2235 & - 2.0d0 * THGCF(IDISP,4) 2236 AVESHG(IDISP,2) = 0.2d0 * ESHGCF(IDISP,1) 2237 & + 0.8d0 * ESHGCF(IDISP,2) 2238 & - 0.2d0 * ESHGCF(IDISP,3) 2239 & - 0.2d0 * ESHGCF(IDISP,4) 2240 AVESHG(IDISP,3) = 1.0d0 * ESHGCF(IDISP,2) 2241 & + 1.0d0 * ESHGCF(IDISP,3) 2242 & - 2.0d0 * ESHGCF(IDISP,4) 2243 AVEKERR(IDISP,2) = 0.2d0 * KERRCF(IDISP,1) 2244 & + 0.8d0 * KERRCF(IDISP,2) 2245 & - 0.2d0 * KERRCF(IDISP,3) 2246 & - 0.2d0 * KERRCF(IDISP,4) 2247 AVEKERR(IDISP,3) = 1.0d0 * KERRCF(IDISP,2) 2248 & + 1.0d0 * KERRCF(IDISP,3) 2249 & - 2.0d0 * KERRCF(IDISP,4) 2250 AVEDFWM(IDISP,2) = 0.2d0 * DFWMCF(IDISP,1) 2251 & + 0.8d0 * DFWMCF(IDISP,2) 2252 & - 0.2d0 * DFWMCF(IDISP,3) 2253 & - 0.2d0 * DFWMCF(IDISP,4) 2254 AVEDFWM(IDISP,3) = 1.0d0 * DFWMCF(IDISP,2) 2255 & + 1.0d0 * DFWMCF(IDISP,3) 2256 & - 2.0d0 * DFWMCF(IDISP,4) 2257 END IF 2258 ELSE IF (CSYM(1:6).EQ.'LINEAR') THEN 2259 IF (GAMMA_PAR) THEN 2260 AVETHG(IDISP,1) = 3.0d0 * THGCF(IDISP,1) 2261 AVESHG(IDISP,1) = 3.0d0 * ESHGCF(IDISP,1) 2262 AVEKERR(IDISP,1) = 3.0d0 * KERRCF(IDISP,1) 2263 AVEDFWM(IDISP,1) = 3.0d0 * DFWMCF(IDISP,1) 2264 2265 AVETHG(IDISP,1) = AVETHG(IDISP,1) + 2266 & 8.0d0 * THGCF(IDISP,8) 2267 AVESHG(IDISP,1) = AVESHG(IDISP,1) + 2268 & 8.0d0 * ESHGCF(IDISP,8) 2269 AVEKERR(IDISP,1) = AVEKERR(IDISP,1) + 2270 & 8.0d0 * KERRCF(IDISP,8) 2271 AVEDFWM(IDISP,1) = AVEDFWM(IDISP,1) + 2272 & 8.0d0 * DFWMCF(IDISP,8) 2273 2274 DO ICOMP = 2, 7 2275 AVETHG(IDISP,1) = AVETHG(IDISP,1) + 2276 & 2.0d0 * THGCF(IDISP,ICOMP) 2277 AVESHG(IDISP,1) = AVESHG(IDISP,1) + 2278 & 2.0d0 * ESHGCF(IDISP,ICOMP) 2279 AVEKERR(IDISP,1) = AVEKERR(IDISP,1) + 2280 & 2.0d0 * KERRCF(IDISP,ICOMP) 2281 AVEDFWM(IDISP,1) = AVEDFWM(IDISP,1) + 2282 & 2.0d0 * DFWMCF(IDISP,ICOMP) 2283 END DO 2284 2285 AVETHG(IDISP,1) = AVETHG(IDISP,1) / 15.0d0 2286 AVESHG(IDISP,1) = AVESHG(IDISP,1) / 15.0d0 2287 AVEKERR(IDISP,1) = AVEKERR(IDISP,1) / 15.0d0 2288 AVEDFWM(IDISP,1) = AVEDFWM(IDISP,1) / 15.0d0 2289 END IF 2290 IF (GAMMA_ORT) THEN 2291 AVETHG(IDISP,2) = 1.0d0 * THGCF(IDISP,1) 2292 & + 4.0d0 * THGCF(IDISP,2) 2293 & - 1.0d0 * THGCF(IDISP,3) 2294 & - 1.0d0 * THGCF(IDISP,4) 2295 & + 4.0d0 * THGCF(IDISP,5) 2296 & - 1.0d0 * THGCF(IDISP,6) 2297 & - 1.0d0 * THGCF(IDISP,7) 2298 & + 1.0d0 * THGCF(IDISP,8) 2299 & + 5.0d0 * THGCF(IDISP,9) 2300 AVETHG(IDISP,2) = AVETHG(IDISP,2) / 15.0d0 2301 AVETHG(IDISP,3) = 0.5d0 * THGCF(IDISP,2) 2302 & + 0.5d0 * THGCF(IDISP,3) 2303 & - 1.0d0 * THGCF(IDISP,4) 2304 & + 0.5d0 * THGCF(IDISP,5) 2305 & + 0.5d0 * THGCF(IDISP,6) 2306 & - 1.0d0 * THGCF(IDISP,7) 2307 & + 1.5d0 * THGCF(IDISP,9) 2308 & + 1.5d0 * THGCF(IDISP,10) 2309 & - 1.0d0 * THGCF(IDISP,8) 2310 AVETHG(IDISP,3) = AVETHG(IDISP,3) / 1.5d0 2311 AVESHG(IDISP,2) = 1.0d0 * ESHGCF(IDISP,1) 2312 & + 4.0d0 * ESHGCF(IDISP,2) 2313 & - 1.0d0 * ESHGCF(IDISP,3) 2314 & - 1.0d0 * ESHGCF(IDISP,4) 2315 & + 4.0d0 * ESHGCF(IDISP,5) 2316 & - 1.0d0 * ESHGCF(IDISP,6) 2317 & - 1.0d0 * ESHGCF(IDISP,7) 2318 & + 1.0d0 * ESHGCF(IDISP,8) 2319 & + 5.0d0 * ESHGCF(IDISP,9) 2320 AVESHG(IDISP,2) = AVESHG(IDISP,2) / 15.0d0 2321 AVESHG(IDISP,3) = 0.5d0 * ESHGCF(IDISP,2) 2322 & + 0.5d0 * ESHGCF(IDISP,3) 2323 & - 1.0d0 * ESHGCF(IDISP,4) 2324 & + 0.5d0 * ESHGCF(IDISP,5) 2325 & + 0.5d0 * ESHGCF(IDISP,6) 2326 & - 1.0d0 * ESHGCF(IDISP,7) 2327 & + 1.5d0 * ESHGCF(IDISP,9) 2328 & + 1.5d0 * ESHGCF(IDISP,10) 2329 & - 1.0d0 * ESHGCF(IDISP,8) 2330 AVESHG(IDISP,3) = AVESHG(IDISP,3) / 1.5d0 2331 AVEKERR(IDISP,2) = 1.0d0 * KERRCF(IDISP,1) 2332 & + 4.0d0 * KERRCF(IDISP,2) 2333 & - 1.0d0 * KERRCF(IDISP,3) 2334 & - 1.0d0 * KERRCF(IDISP,4) 2335 & + 4.0d0 * KERRCF(IDISP,5) 2336 & - 1.0d0 * KERRCF(IDISP,6) 2337 & - 1.0d0 * KERRCF(IDISP,7) 2338 & + 1.0d0 * KERRCF(IDISP,8) 2339 & + 5.0d0 * KERRCF(IDISP,9) 2340 AVEKERR(IDISP,2) = AVEKERR(IDISP,2) / 15.0d0 2341 AVEKERR(IDISP,3) = 0.5d0 * KERRCF(IDISP,2) 2342 & + 0.5d0 * KERRCF(IDISP,3) 2343 & - 1.0d0 * KERRCF(IDISP,4) 2344 & + 0.5d0 * KERRCF(IDISP,5) 2345 & + 0.5d0 * KERRCF(IDISP,6) 2346 & - 1.0d0 * KERRCF(IDISP,7) 2347 & + 1.5d0 * KERRCF(IDISP,9) 2348 & + 1.5d0 * KERRCF(IDISP,10) 2349 & - 1.0d0 * KERRCF(IDISP,8) 2350 AVEKERR(IDISP,3) = AVEKERR(IDISP,3) / 1.5d0 2351 AVEDFWM(IDISP,2) = 1.0d0 * DFWMCF(IDISP,1) 2352 & + 4.0d0 * DFWMCF(IDISP,2) 2353 & - 1.0d0 * DFWMCF(IDISP,3) 2354 & - 1.0d0 * DFWMCF(IDISP,4) 2355 & + 4.0d0 * DFWMCF(IDISP,5) 2356 & - 1.0d0 * DFWMCF(IDISP,6) 2357 & - 1.0d0 * DFWMCF(IDISP,7) 2358 & + 1.0d0 * DFWMCF(IDISP,8) 2359 & + 5.0d0 * DFWMCF(IDISP,9) 2360 AVEDFWM(IDISP,2) = AVEDFWM(IDISP,2) / 15.0d0 2361 AVEDFWM(IDISP,3) = 0.5d0 * DFWMCF(IDISP,2) 2362 & + 0.5d0 * DFWMCF(IDISP,3) 2363 & - 1.0d0 * DFWMCF(IDISP,4) 2364 & + 0.5d0 * DFWMCF(IDISP,5) 2365 & + 0.5d0 * DFWMCF(IDISP,6) 2366 & - 1.0d0 * DFWMCF(IDISP,7) 2367 & + 1.5d0 * DFWMCF(IDISP,9) 2368 & + 1.5d0 * DFWMCF(IDISP,10) 2369 & - 1.0d0 * DFWMCF(IDISP,8) 2370 AVEDFWM(IDISP,3) = AVEDFWM(IDISP,3) / 1.5d0 2371 END IF 2372 ELSE 2373 IF (GAMMA_PAR) THEN 2374 AVETHG(IDISP,1) = 2.0d0*THGCF(IDISP,1) 2375 AVESHG(IDISP,1) = 2.0d0*ESHGCF(IDISP,1) 2376 AVEKERR(IDISP,1)= 2.0d0*KERRCF(IDISP,1) 2377 AVEDFWM(IDISP,1)= 2.0d0*DFWMCF(IDISP,1) 2378 2379 AVETHG(IDISP,1) =AVETHG(IDISP,1) +2.0d0*THGCF(IDISP,8) 2380 AVESHG(IDISP,1) =AVESHG(IDISP,1) +2.0d0*ESHGCF(IDISP,8) 2381 AVEKERR(IDISP,1)=AVEKERR(IDISP,1)+2.0d0*KERRCF(IDISP,8) 2382 AVEDFWM(IDISP,1)=AVEDFWM(IDISP,1)+2.0d0*DFWMCF(IDISP,8) 2383 2384 AVETHG(IDISP,1) =AVETHG(IDISP,1) +2.0d0*THGCF(IDISP,15) 2385 AVESHG(IDISP,1) =AVESHG(IDISP,1) +2.0d0*ESHGCF(IDISP,15) 2386 AVEKERR(IDISP,1)=AVEKERR(IDISP,1)+2.0d0*KERRCF(IDISP,15) 2387 AVEDFWM(IDISP,1)=AVEDFWM(IDISP,1)+2.0d0*DFWMCF(IDISP,15) 2388 2389 DO ICOMP = 1, 21 2390 AVETHG(IDISP,1) =AVETHG(IDISP,1) + THGCF(IDISP,ICOMP) 2391 AVESHG(IDISP,1) =AVESHG(IDISP,1) + ESHGCF(IDISP,ICOMP) 2392 AVEKERR(IDISP,1)=AVEKERR(IDISP,1) + KERRCF(IDISP,ICOMP) 2393 AVEDFWM(IDISP,1)=AVEDFWM(IDISP,1) + DFWMCF(IDISP,ICOMP) 2394 END DO 2395 2396 AVETHG(IDISP,1) = AVETHG(IDISP,1) / 15.0d0 2397 AVESHG(IDISP,1) = AVESHG(IDISP,1) / 15.0d0 2398 AVEKERR(IDISP,1) = AVEKERR(IDISP,1) / 15.0d0 2399 AVEDFWM(IDISP,1) = AVEDFWM(IDISP,1) / 15.0d0 2400 END IF 2401 IF (GAMMA_ORT) THEN 2402 AVETHG(IDISP,2) = 0.0d0 2403 DO IADD = 0, 14, 7 2404 AVETHG(IDISP,2) = AVETHG(IDISP,2) + 2405 & 1.0d0 * THGCF(IDISP,1+IADD) 2406 & +2.0d0 * THGCF(IDISP,2+IADD) 2407 & -0.5d0 * THGCF(IDISP,3+IADD) 2408 & -0.5d0 * THGCF(IDISP,4+IADD) 2409 & +2.0d0 * THGCF(IDISP,5+IADD) 2410 & -0.5d0 * THGCF(IDISP,6+IADD) 2411 & -0.5d0 * THGCF(IDISP,7+IADD) 2412 END DO 2413 AVETHG(IDISP,2) = AVETHG(IDISP,2) / 15.0d0 2414 AVETHG(IDISP,3) = 0.0d0 2415 DO IADD = 0, 14, 7 2416 AVETHG(IDISP,3) = AVETHG(IDISP,3) + 2417 & 0.5d0 * THGCF(IDISP,2+IADD) 2418 & + 0.5d0 * THGCF(IDISP,3+IADD) 2419 & - 1.0d0 * THGCF(IDISP,4+IADD) 2420 & + 0.5d0 * THGCF(IDISP,5+IADD) 2421 & + 0.5d0 * THGCF(IDISP,6+IADD) 2422 & - 1.0d0 * THGCF(IDISP,7+IADD) 2423 END DO 2424 AVETHG(IDISP,3) = AVETHG(IDISP,3) / 3.0d0 2425 2426 AVESHG(IDISP,2) = 0.0d0 2427 DO IADD = 0, 14, 7 2428 AVESHG(IDISP,2) = AVESHG(IDISP,2) + 2429 & 1.0d0 * ESHGCF(IDISP,1+IADD) 2430 & +2.0d0 * ESHGCF(IDISP,2+IADD) 2431 & -0.5d0 * ESHGCF(IDISP,3+IADD) 2432 & -0.5d0 * ESHGCF(IDISP,4+IADD) 2433 & +2.0d0 * ESHGCF(IDISP,5+IADD) 2434 & -0.5d0 * ESHGCF(IDISP,6+IADD) 2435 & -0.5d0 * ESHGCF(IDISP,7+IADD) 2436 END DO 2437 AVESHG(IDISP,2) = AVESHG(IDISP,2) / 15.0d0 2438 AVESHG(IDISP,3) = 0.0d0 2439 DO IADD = 0, 14, 7 2440 AVESHG(IDISP,3) = AVESHG(IDISP,3) + 2441 & 0.5d0 * ESHGCF(IDISP,2+IADD) 2442 & + 0.5d0 * ESHGCF(IDISP,3+IADD) 2443 & - 1.0d0 * ESHGCF(IDISP,4+IADD) 2444 & + 0.5d0 * ESHGCF(IDISP,5+IADD) 2445 & + 0.5d0 * ESHGCF(IDISP,6+IADD) 2446 & - 1.0d0 * ESHGCF(IDISP,7+IADD) 2447 END DO 2448 AVESHG(IDISP,3) = AVESHG(IDISP,3) / 3.0d0 2449 2450 AVEKERR(IDISP,2) = 0.0d0 2451 DO IADD = 0, 14, 7 2452 AVEKERR(IDISP,2) = AVEKERR(IDISP,2) + 2453 & 1.0d0 * KERRCF(IDISP,1+IADD) 2454 & +2.0d0 * KERRCF(IDISP,2+IADD) 2455 & -0.5d0 * KERRCF(IDISP,3+IADD) 2456 & -0.5d0 * KERRCF(IDISP,4+IADD) 2457 & +2.0d0 * KERRCF(IDISP,5+IADD) 2458 & -0.5d0 * KERRCF(IDISP,6+IADD) 2459 & -0.5d0 * KERRCF(IDISP,7+IADD) 2460 END DO 2461 AVEKERR(IDISP,2) = AVEKERR(IDISP,2) / 15.0d0 2462 AVEKERR(IDISP,3) = 0.0d0 2463 DO IADD = 0, 14, 7 2464 AVEKERR(IDISP,3) = AVEKERR(IDISP,3) + 2465 & 0.5d0 * KERRCF(IDISP,2+IADD) 2466 & + 0.5d0 * KERRCF(IDISP,3+IADD) 2467 & - 1.0d0 * KERRCF(IDISP,4+IADD) 2468 & + 0.5d0 * KERRCF(IDISP,5+IADD) 2469 & + 0.5d0 * KERRCF(IDISP,6+IADD) 2470 & - 1.0d0 * KERRCF(IDISP,7+IADD) 2471 END DO 2472 AVEKERR(IDISP,3) = AVEKERR(IDISP,3) / 3.0d0 2473 2474 AVEDFWM(IDISP,2) = 0.0d0 2475 DO IADD = 0, 14, 7 2476 AVEDFWM(IDISP,2) = AVEDFWM(IDISP,2) + 2477 & 1.0d0 * DFWMCF(IDISP,1+IADD) 2478 & +2.0d0 * DFWMCF(IDISP,2+IADD) 2479 & -0.5d0 * DFWMCF(IDISP,3+IADD) 2480 & -0.5d0 * DFWMCF(IDISP,4+IADD) 2481 & +2.0d0 * DFWMCF(IDISP,5+IADD) 2482 & -0.5d0 * DFWMCF(IDISP,6+IADD) 2483 & -0.5d0 * DFWMCF(IDISP,7+IADD) 2484 END DO 2485 AVEDFWM(IDISP,2) = AVEDFWM(IDISP,2) / 15.0d0 2486 AVEDFWM(IDISP,3) = 0.0d0 2487 DO IADD = 0, 14, 7 2488 AVEDFWM(IDISP,3) = AVEDFWM(IDISP,3) + 2489 & 0.5d0 * DFWMCF(IDISP,2+IADD) 2490 & + 0.5d0 * DFWMCF(IDISP,3+IADD) 2491 & - 1.0d0 * DFWMCF(IDISP,4+IADD) 2492 & + 0.5d0 * DFWMCF(IDISP,5+IADD) 2493 & + 0.5d0 * DFWMCF(IDISP,6+IADD) 2494 & - 1.0d0 * DFWMCF(IDISP,7+IADD) 2495 END DO 2496 AVEDFWM(IDISP,3) = AVEDFWM(IDISP,3) / 3.0d0 2497 END IF 2498 END IF 2499 END DO 2500 2501 IF (LOCDBG) THEN 2502 WRITE (LUPRI,'(//,A)') 2503 & 'DISPERSION COEFF. FOR PARALLEL/ORTHOGONAL/MS AVERAGE:' 2504 DO LMN = 0, NCRDSPE 2505 DO MN = 0, LMN 2506 DO N = 0, MN 2507 L = LMN - MN 2508 M = MN - N 2509 WRITE (LUPRI,'(3I5,4G16.8)') 2510 & L,M,N,(AVEDISP(ITRI(LMN,MN,N),I),I=1,3) 2511 END DO 2512 END DO 2513 END DO 2514 2515 WRITE (LUPRI,'(//,A)') 2516 & 'THG COEFFICIENTS FOR PARALLEL/ORTHOGONAL/MS AVERAGE:' 2517 DO IDISP = 1, NCRDSPE+1 2518 WRITE (LUPRI,'(I5,4G16.8)') IDISP-1,(AVETHG(IDISP,I),I=1,3) 2519 END DO 2520 2521 WRITE (LUPRI,'(//,A)') 2522 & 'ESHG COEFFICIENTS FOR PARALLEL/ORTHOGONAL/MS AVERAGE:' 2523 DO IDISP = 1, NCRDSPE+1 2524 WRITE (LUPRI,'(I5,4G16.8)') IDISP-1,(AVESHG(IDISP,I),I=1,3) 2525 END DO 2526 2527 WRITE (LUPRI,'(//,A)') 2528 & 'KERR COEFFICIENTS FOR PARALLEL/ORTHOGONAL/MS AVERAGE:' 2529 DO IDISP = 1, NCRDSPE+1 2530 WRITE (LUPRI,'(I5,4G16.8)') IDISP-1,(AVEKERR(IDISP,I),I=1,3) 2531 END DO 2532 2533 WRITE (LUPRI,'(//,A)') 2534 & 'DFWM COEFFICIENTS FOR PARALLEL/ORTHOGONAL/MS AVERAGE:' 2535 DO IDISP = 1, NCRDSPE+1 2536 WRITE (LUPRI,'(I5,4G16.8)') IDISP-1,(AVEDFWM(IDISP,I),I=1,3) 2537 END DO 2538 END IF 2539 2540* calculate A, B, C, etc. coefficients: 2541 GAMMA0 = AVEDISP(ITRI(0,0,0),1) 2542 2543 IF (NCRDSPE.GE.2) THEN 2544 2545 A = AVEDISP(ITRI(2,0,0),1) / (2.0d0 * GAMMA0) 2546 ERR1 = AVEDISP(ITRI(2,1,0),1) / (2.0d0 * GAMMA0) - A 2547 ERROR = DABS(ERR1) / DABS(A) 2548 IF (ERROR .GT. THRSDISP) THEN 2549 WRITE (LUPRI,*) 2550 & 'WARNING: LARGE NUMERICAL (?) ERRORS DURING', 2551 & ' CALCULATION OF A COEFFICIENT:' 2552 WRITE (LUPRI,*) 'TEST #1 GAVE:',ERR1 2553 LERROR = .TRUE. 2554 END IF 2555 ABCDE(KA,1) = A 2556 IF (LOCDBG) THEN 2557 WRITE (LUPRI,'(//,A)') 2558 & 'A,B,C,D,E COEFF. FOR PARALLEL AVERAGE:' 2559 WRITE (LUPRI,*) " A coefficient :",A 2560 END IF 2561 2562 IF (GAMMA_ORT) THEN 2563 A = AVEDISP(ITRI(2,0,0),3) / GAMMA0 2564 ERR1 = AVEDISP(ITRI(2,1,0),3) / (-2.0d0 * GAMMA0) - A 2565 ERROR = DABS(ERR1) / DABS(A) 2566 IF (ERROR .GT. THRSDISP) THEN 2567 WRITE (LUPRI,*) 2568 & 'WARNING: LARGE NUMERICAL (?) ERRORS DURING', 2569 & ' CALCULATION OF A_ms COEFFICIENT:' 2570 WRITE (LUPRI,*) 'TEST #1 GAVE:',ERR1 2571 LERROR = .TRUE. 2572 END IF 2573 ABCDE(KA,2) = A 2574 IF (LOCDBG) THEN 2575 WRITE (LUPRI,*) " A_ms coefficient :",A 2576 END IF 2577 END IF 2578 END IF 2579 2580 IF (NCRDSPE.GE.4) THEN 2581 B = 0.25d0*AVEDISP(ITRI(4,0,0),1) 2582 BP=0.25d0*AVEDISP(ITRI(4,2,1),1)-0.5d0*AVEDISP(ITRI(4,1,0),1) 2583 ERR1=AVEDISP(ITRI(4,1,0),1) - 8.0d0 * B 2584 ERR2=AVEDISP(ITRI(4,2,0),1) - 12.0d0 * B 2585 ERROR = (DABS(ERR1)+DABS(ERR2)) / DABS(2.0d0*B) 2586 B = B / GAMMA0 2587 BP = BP / GAMMA0 2588 ABCDE(KB,1) = B 2589 ABCDE(KBP,1) = BP 2590 IF (ERROR .GT. THRSDISP) THEN 2591 WRITE (LUPRI,*) 2592 & 'WARNING: LARGE NUMERICAL (?) ERRORS DURING', 2593 & ' CALCULATION OF B COEFFICIENTS:' 2594 WRITE (LUPRI,*) 'TEST #1 GAVE:',ERR1 2595 WRITE (LUPRI,*) 'TEST #2 GAVE:',ERR2 2596 LERROR = .TRUE. 2597 END IF 2598 IF (LOCDBG) THEN 2599 WRITE(LUPRI,*)" B, B' coefficients:",B,BP 2600 END IF 2601 2602 2603 IF (GAMMA_ORT) THEN 2604C B = AVEDISP(ITRI(4,2,2),3) / 24.0d0 2605C BP = 2.0d0 * B - 0.5d0 * AVEDISP(ITRI(4,0,0),3) 2606C ERR1 = AVEDISP(ITRI(4,0,0),3) - ( 4.0d0*B - 2.0d0*BP) 2607C ERR2 = AVEDISP(ITRI(4,1,0),3) - ( -4.0d0*B - 4.0d0*BP) 2608C ERR3 = AVEDISP(ITRI(4,1,1),3) - ( 20.0d0*B - 4.0d0*BP) 2609C ERR4 = AVEDISP(ITRI(4,2,0),3) - (-12.0d0*B - 0.0d0*BP) 2610C ERR5 = AVEDISP(ITRI(4,2,1),3) - ( 4.0d0*B - 4.0d0*BP) 2611C ERR6 = AVEDISP(ITRI(4,2,2),3) - ( 24.0d0*B - 0.0d0*BP) 2612C ERR7 = AVEDISP(ITRI(4,3,0),3) - (-16.0d0*B + 8.0d0*BP) 2613C ERR8 = AVEDISP(ITRI(4,3,1),3) - ( -8.0d0*B + 8.0d0*BP) 2614C ERR9 = AVEDISP(ITRI(4,4,0),3) - ( -8.0d0*B + 4.0d0*BP) 2615 B = AVEDISP(ITRI(4,0,0),3) / 2.0d0 2616 BP =-(AVEDISP(ITRI(4,1,0),3)+AVEDISP(ITRI(4,0,0),3))/9.0d0 2617 ERR1 = AVEDISP(ITRI(4,0,0),3) - ( 2.0d0*B + 0.0d0*BP) 2618 ERR2 = AVEDISP(ITRI(4,1,0),3) - ( -2.0d0*B - 9.0d0*BP) 2619 ERR3 = AVEDISP(ITRI(4,1,1),3) - ( 10.0d0*B + 9.0d0*BP) 2620 ERR4 = AVEDISP(ITRI(4,2,0),3) - ( -6.0d0*B - 9.0d0*BP) 2621 ERR5 = AVEDISP(ITRI(4,2,1),3) - ( 2.0d0*B - 3.0d0*BP) 2622 ERR6 = AVEDISP(ITRI(4,2,2),3) - ( 12.0d0*B +18.0d0*BP) 2623 ERR7 = AVEDISP(ITRI(4,3,0),3) - ( -8.0d0*B + 0.0d0*BP) 2624 ERR8 = AVEDISP(ITRI(4,3,1),3) - ( -4.0d0*B + 6.0d0*BP) 2625 ERR9 = AVEDISP(ITRI(4,4,0),3) - ( -4.0d0*B + 0.0d0*BP) 2626 ERROR = (DABS(ERR1) + DABS(ERR2) + DABS(ERR3) + DABS(ERR4)+ 2627 & DABS(ERR5) + DABS(ERR6) + DABS(ERR7) + DABS(ERR8)+ 2628 & DABS(ERR9) ) / DABS(9.0d0*B) 2629 IF (ERROR .GT. THRSDISP) THEN 2630 WRITE(LUPRI,*) 2631 & 'WARNING: LARGE NUMERICAL (?) ERRORS DURING', 2632 & ' CALCULATION OF B_ms COEFFICIENTS:' 2633 WRITE (LUPRI,*) 'TEST #1 GAVE:',ERR1 2634 WRITE (LUPRI,*) 'TEST #2 GAVE:',ERR2 2635 WRITE (LUPRI,*) 'TEST #3 GAVE:',ERR3 2636 WRITE (LUPRI,*) 'TEST #4 GAVE:',ERR4 2637 WRITE (LUPRI,*) 'TEST #5 GAVE:',ERR5 2638 WRITE (LUPRI,*) 'TEST #6 GAVE:',ERR6 2639 WRITE (LUPRI,*) 'TEST #7 GAVE:',ERR7 2640 WRITE (LUPRI,*) 'TEST #8 GAVE:',ERR8 2641 WRITE (LUPRI,*) 'TEST #9 GAVE:',ERR9 2642 LERROR = .TRUE. 2643 END IF 2644 B = B / GAMMA0 2645 BP = BP / GAMMA0 2646 ABCDE(KB,2) = B 2647 ABCDE(KBP,2) = BP 2648 IF (LOCDBG) THEN 2649 WRITE (LUPRI,*) 2650 & " B_ms, B_ms' coefficients:",B,BP 2651 END IF 2652 END IF 2653 END IF 2654 2655 IF (NCRDSPE.GE.6) THEN 2656 C = AVEDISP(ITRI(6,0,0),1) / 8.0d0 2657 CP =(AVEDISP(ITRI(6,2,0),1) 2658 & - 6.0d0*AVEDISP(ITRI(6,0,0),1) ) / 9.0d0 2659 CPP =(AVEDISP(ITRI(6,2,1),1) -2.0d0*AVEDISP(ITRI(6,2,0),1) 2660 & + AVEDISP(ITRI(6,1,0),1) ) / 8.0d0 2661 ERR1 = AVEDISP(ITRI(6,1,0),1)-24.0d0*C 2662 ERR2 = AVEDISP(ITRI(6,3,0),1)-56.0d0*C-18.0d0*CP 2663 ERR3 = AVEDISP(ITRI(6,3,1),1)-120.0d0*C-16.0d0*CPP-54.0d0*CP 2664 ERR4 = AVEDISP(ITRI(6,4,2),1)-168.0d0*C-24.0d0*CPP-90.0d0*CP 2665 C = C / GAMMA0 2666 CP = CP / GAMMA0 2667 CPP = CPP / GAMMA0 2668 ERROR = ( DABS(ERR1) + DABS(ERR2) + DABS(ERR3) 2669 & + DABS(ERR4) ) / DABS(4.0d0*C) 2670 IF (ERROR .GT. THRSDISP) THEN 2671 WRITE (LUPRI,*) 2672 & 'WARNING: LARGE NUMERICAL (?) ERRORS DURING', 2673 & ' CALCULATION OF C COEFFICIENTS:' 2674 WRITE (LUPRI,*) 'TEST #1 GAVE:',ERR1 2675 WRITE (LUPRI,*) 'TEST #2 GAVE:',ERR2 2676 WRITE (LUPRI,*) 'TEST #3 GAVE:',ERR3 2677 WRITE (LUPRI,*) 'TEST #4 GAVE:',ERR4 2678 LERROR = .TRUE. 2679 END IF 2680 ABCDE(KC,1) = C 2681 ABCDE(KCP,1) = CP 2682 ABCDE(KCPP,1) = CPP 2683 IF (LOCDBG) THEN 2684 WRITE (LUPRI,*) " C, C', C'' coefficients:", 2685 & C, CP, CPP 2686 END IF 2687 2688 2689 IF (GAMMA_ORT) THEN 2690C CPP = - AVEDISP(ITRI(6,1,0),3) / 12.0d0 2691C CP = AVEDISP(ITRI(6,3,1),3)/16.0d0 + 0.50d0 * CPP 2692C C = (AVEDISP(ITRI(6,0,0),3)-AVEDISP(ITRI(6,1,0),3)/3.0d0 2693C & -4.0d0 * CP ) / 8.0d0 2694C ERR1 = ( 8.d0*C+ 4.d0*CP- 4.d0*CPP)-AVEDISP(ITRI(6,0,0),3) 2695C ERR2 = ( 0.d0*C+ 0.d0*CP-12.d0*CPP)-AVEDISP(ITRI(6,1,0),3) 2696C ERR3 = ( 48.d0*C+24.d0*CP-12.d0*CPP)-AVEDISP(ITRI(6,1,1),3) 2697C ERR4 = (-24.d0*C-12.d0*CP-12.d0*CPP)-AVEDISP(ITRI(6,2,0),3) 2698C ERR5 = ( 48.d0*C+ 8.d0*CP-28.d0*CPP)-AVEDISP(ITRI(6,2,1),3) 2699C ERR6 = ( 96.d0*C+48.d0*CP-12.d0*CPP)-AVEDISP(ITRI(6,2,2),3) 2700C ERR7 = (-64.d0*C-32.d0*CP+ 8.d0*CPP)-AVEDISP(ITRI(6,3,0),3) 2701C ERR8 = ( 0.d0*C+16.d0*CP- 8.d0*CPP)-AVEDISP(ITRI(6,3,1),3) 2702C ERR9 = ( 96.d0*C-32.d0*CP-32.d0*CPP)-AVEDISP(ITRI(6,3,2),3) 2703C ERR10 =(128.d0*C+64.d0*CP-16.d0*CPP)-AVEDISP(ITRI(6,3,3),3) 2704C ERR11 =(-72.d0*C-36.d0*CP+24.d0*CPP)-AVEDISP(ITRI(6,4,0),3) 2705C ERR12 =(-96.d0*C+16.d0*CP+40.d0*CPP)-AVEDISP(ITRI(6,4,1),3) 2706C ERR13 =( 0.d0*C+ 0.d0*CP+ 0.d0*CPP)-AVEDISP(ITRI(6,4,2),3) 2707C ERR14 =(-48.d0*C-24.d0*CP+24.d0*CPP)-AVEDISP(ITRI(6,5,0),3) 2708C ERR15 =(-96.d0*C-16.d0*CP+56.d0*CPP)-AVEDISP(ITRI(6,5,1),3) 2709C ERR16 =(-16.d0*C- 8.d0*CP+ 8.d0*CPP)-AVEDISP(ITRI(6,6,0),3) 2710 C = AVEDISP(ITRI(6,0,0),3) / 4.0d0 2711 CP = -AVEDISP(ITRI(6,1,0),3) / 18.0d0 2712 CPP = ( 3.0d0*AVEDISP(ITRI(6,3,1),3) - 2713 & 2.0d0*AVEDISP(ITRI(6,1,0),3) ) / 24.0d0 2714 ERR1 = ( 4.d0*C+ 0.d0*CPP+ 0.d0*CP)-AVEDISP(ITRI(6,0,0),3) 2715 ERR2 = ( 0.d0*C+ 0.d0*CPP-18.d0*CP)-AVEDISP(ITRI(6,1,0),3) 2716 ERR3 = ( 24.d0*C+ 0.d0*CPP+18.d0*CP)-AVEDISP(ITRI(6,1,1),3) 2717 ERR4 = (-12.d0*C+ 0.d0*CPP-36.d0*CP)-AVEDISP(ITRI(6,2,0),3) 2718 ERR5 = ( 24.d0*C- 8.d0*CPP- 6.d0*CP)-AVEDISP(ITRI(6,2,1),3) 2719 ERR6 = ( 48.d0*C+ 0.d0*CPP+54.d0*CP)-AVEDISP(ITRI(6,2,2),3) 2720 ERR7 = (-32.d0*C+ 0.d0*CPP-36.d0*CP)-AVEDISP(ITRI(6,3,0),3) 2721 ERR8 = ( 0.d0*C+ 8.d0*CPP-12.d0*CP)-AVEDISP(ITRI(6,3,1),3) 2722 ERR9 = ( 48.d0*C-40.d0*CPP+24.d0*CP)-AVEDISP(ITRI(6,3,2),3) 2723 ERR10 =( 64.d0*C+ 0.d0*CPP+72.d0*CP)-AVEDISP(ITRI(6,3,3),3) 2724 ERR11 =(-36.d0*C+ 0.d0*CPP-18.d0*CP)-AVEDISP(ITRI(6,4,0),3) 2725 ERR12 =(-48.d0*C+32.d0*CPP-12.d0*CP)-AVEDISP(ITRI(6,4,1),3) 2726 ERR13 =( 0.d0*C+ 0.d0*CPP+ 0.d0*CP)-AVEDISP(ITRI(6,4,2),3) 2727 ERR14 =(-24.d0*C+ 0.d0*CPP+ 0.d0*CP)-AVEDISP(ITRI(6,5,0),3) 2728 ERR15 =(-48.d0*C+16.d0*CPP+12.d0*CP)-AVEDISP(ITRI(6,5,1),3) 2729 ERR16 =( -8.d0*C+ 0.d0*CPP+ 0.d0*CP)-AVEDISP(ITRI(6,6,0),3) 2730 ERROR = ( DABS(ERR1) + DABS(ERR2) + DABS(ERR3) + 2731 & DABS(ERR4) + DABS(ERR5) + DABS(ERR6) + 2732 & DABS(ERR7) + DABS(ERR8) + DABS(ERR9) + 2733 & DABS(ERR10)+ DABS(ERR11)+ DABS(ERR12)+ 2734 & DABS(ERR13)+ DABS(ERR14)+ DABS(ERR15)+ 2735 & DABS(ERR16) ) / DABS(16.0d0*C) 2736 IF ( ERROR .GT. THRSDISP) THEN 2737 WRITE (LUPRI,*) 2738 & 'WARNING: LARGE NUMERICAL (?) ERRORS DURING', 2739 & ' CALCULATION OF C_ms COEFFICIENTS:' 2740 WRITE (LUPRI,*) 'TEST #1 GAVE:',ERR1 2741 WRITE (LUPRI,*) 'TEST #2 GAVE:',ERR2 2742 WRITE (LUPRI,*) 'TEST #3 GAVE:',ERR3 2743 WRITE (LUPRI,*) 'TEST #4 GAVE:',ERR4 2744 WRITE (LUPRI,*) 'TEST #5 GAVE:',ERR5 2745 WRITE (LUPRI,*) 'TEST #6 GAVE:',ERR6 2746 WRITE (LUPRI,*) 'TEST #7 GAVE:',ERR7 2747 WRITE (LUPRI,*) 'TEST #8 GAVE:',ERR8 2748 WRITE (LUPRI,*) 'TEST #9 GAVE:',ERR9 2749 WRITE (LUPRI,*) 'TEST #10 GAVE:',ERR10 2750 WRITE (LUPRI,*) 'TEST #11 GAVE:',ERR11 2751 WRITE (LUPRI,*) 'TEST #12 GAVE:',ERR12 2752 WRITE (LUPRI,*) 'TEST #13 GAVE:',ERR13 2753 WRITE (LUPRI,*) 'TEST #14 GAVE:',ERR14 2754 WRITE (LUPRI,*) 'TEST #15 GAVE:',ERR15 2755 WRITE (LUPRI,*) 'TEST #16 GAVE:',ERR16 2756 LERROR = .TRUE. 2757 END IF 2758 C = C / GAMMA0 2759 CP = CP / GAMMA0 2760 CPP = CPP / GAMMA0 2761 ABCDE(KC,2) = C 2762 ABCDE(KCP,2) = CP 2763 ABCDE(KCPP,2) = CPP 2764 IF (LOCDBG) THEN 2765 WRITE (LUPRI,*) 2766 & " C_ms, C_ms', Cms'' coefficients:",C,CP,CPP 2767 END IF 2768 END IF 2769 2770 END IF 2771 2772 IF (NCRDSPE.GE.8) THEN 2773 D = AVEDISP(ITRI(8,0,0),1) / 16.0d0 2774 DP = (AVEDISP(ITRI(8,2,1),1) -2.d0*AVEDISP(ITRI(8,2,0),1) 2775 & + AVEDISP(ITRI(8,1,0),1) ) / 16.0d0 2776 DPP = (AVEDISP(ITRI(8,4,2),1) - AVEDISP(ITRI(8,4,1),1) 2777 & -AVEDISP(ITRI(8,3,1),1)+4.d0*AVEDISP(ITRI(8,1,0),1))/16.d0 2778 DPPP = (AVEDISP(ITRI(8,2,0),1) 2779 & - 10.0d0 * AVEDISP(ITRI(8,0,0),1) ) / 18.0d0 2780 ERR1 = AVEDISP(ITRI(8,1,0),1) - 64.0d0*D 2781 ERR2 = AVEDISP(ITRI(8,3,0),1) - (256.0d0*D + 54.0d0*DPPP) 2782 ERR3 = AVEDISP(ITRI(8,3,1),1) - 2783 & ( 576.0d0*D + 48.0d0*DP + 162.0d0*DPPP) 2784 ERR4 = AVEDISP(ITRI(8,4,0),1) - (304.0d0*D + 72.0d0*DPPP) 2785 ERR5 = AVEDISP(ITRI(8,4,1),1) - 2786 & ( 832.0d0*D + 80.0d0*DP + 306.0d0*DPPP) 2787 ERR6 = AVEDISP(ITRI(8,5,2),1) - 2788 & (1408.0d0*D + 176.0d0*DP + 32.0d0*DPP + 648.0d0*DPPP) 2789 2790 DPPP = DPPP / GAMMA0 2791 DPP = DPP / GAMMA0 2792 DP = DP / GAMMA0 2793 D = D / GAMMA0 2794 2795 ERROR = ( DABS(ERR1) + DABS(ERR2) + DABS(ERR3) + DABS(ERR4) + 2796 & DABS(ERR5) + DABS(ERR6) ) / DABS(6.0d0*D) 2797 IF ( ERROR .GT. THRSDISP) THEN 2798 WRITE (LUPRI,*) 2799 & 'WARNING: LARGE NUMERICAL (?) ERRORS DURING', 2800 & ' CALCULATION OF D COEFFICIENTS:' 2801 WRITE (LUPRI,*) 'TEST #1 GAVE:',ERR1 2802 WRITE (LUPRI,*) 'TEST #2 GAVE:',ERR2 2803 WRITE (LUPRI,*) 'TEST #3 GAVE:',ERR3 2804 WRITE (LUPRI,*) 'TEST #4 GAVE:',ERR4 2805 WRITE (LUPRI,*) 'TEST #5 GAVE:',ERR5 2806 WRITE (LUPRI,*) 'TEST #6 GAVE:',ERR6 2807 LERROR = .TRUE. 2808 END IF 2809 ABCDE(KD,1) = D 2810 ABCDE(KDP,1) = DP 2811 ABCDE(KDPP,1) = DPP 2812 ABCDE(KDPPP,1) = DPPP 2813 IF (LOCDBG) THEN 2814 WRITE (LUPRI,*) " D, D', D'', D''' coefficients:", 2815 & D, DP, DPP, DPPP 2816 END IF 2817 END IF 2818 2819 IF (NCRDSPE.GE.10) THEN 2820 E = AVEDISP(ITRI(10,0,0),1) / 32.0d0 2821 EP =(AVEDISP(ITRI(10,2,1),1)-2.d0*AVEDISP(ITRI(10,2,0),1) 2822 & + AVEDISP(ITRI(10,1,0),1) ) / 32.0d0 2823 EPP = (AVEDISP(ITRI(10,4,2),1)-2.0d0*AVEDISP(ITRI(10,4,1),1) 2824 & +3.0d0*AVEDISP(ITRI(10,2,1),1)+10.0d0*AVEDISP(ITRI(10,2,0),1) 2825 & -29.0d0*AVEDISP(ITRI(10,1,0),1) ) / 32.0d0 2826 EPPP = (AVEDISP(ITRI(10,2,0),1) 2827 & -15.0d0 * AVEDISP(ITRI(10,0,0),1) ) / 36.0d0 2828 EPPPP=(AVEDISP(ITRI(10,4,1),1)-9.0d0*AVEDISP(ITRI(10,2,1),1) 2829 & -14.0d0*AVEDISP(ITRI(10,2,0),1) 2830 & +61.0d0*AVEDISP(ITRI(10,1,0),1) ) / 36.0d0 2831 2832 ERR1 = 32.d0*E - AVEDISP(ITRI(10,0,0),1) 2833 ERR2 = 160.d0*E - AVEDISP(ITRI(10,1,0),1) 2834 ERR3 = 480.d0*E+ 36.d0*EPPP - AVEDISP(ITRI(10,2,0),1) 2835 ERR4 = 800.d0*E+ 32.d0*EP+ 72.d0*EPPP 2836 & - AVEDISP(ITRI(10,2,1),1) 2837 ERR5 = 960.d0*E+ 144.d0*EPPP - AVEDISP(ITRI(10,3,0),1) 2838 ERR6 = 2240.d0*E+128.d0*EP+ 432.d0*EPPP 2839 & - AVEDISP(ITRI(10,3,1),1) 2840 ERR7 = 1440.d0*E+288.d0*EPPP - AVEDISP(ITRI(10,4,0),1) 2841 ERR8 = 4160.d0*E+288.d0*EP+1152.d0*EPPP+ 36.d0*EPPPP 2842 & - AVEDISP(ITRI(10,4,1),1) 2843 ERR9 = 5760.d0*E+480.d0*EP+32.d0*EPP+1728.d0*EPPP+ 2844 & 72.d0*EPPPP - AVEDISP(ITRI(10,4,2),1) 2845 ERR10 = 1632.d0*E+ 360.d0*EPPP - AVEDISP(ITRI(10,5,0),1) 2846 ERR11= 5600.d0*E+416.d0*EP+1800.d0*EPPP+ 108.d0*EPPPP 2847 & - AVEDISP(ITRI(10,5,1),1) 2848 ERR12= 9600.d0*E+960.d0*EP+96.d0*EPP+3600.d0*EPPP+ 2849 & 324.d0*EPPPP - AVEDISP(ITRI(10,5,2),1) 2850 ERR13=11360.d0*E+1184.d0*EP+128.d0*EPP+4536.d0*EPPP+ 2851 & 504.d0*EPPPP - AVEDISP(ITRI(10,6,2),1) 2852 ERR14=14080.d0*E+1632.d0*EP+224.d0*EPP+6048.d0*EPPP+ 2853 & 792.d0*EPPPP - AVEDISP(ITRI(10,6,3),1) 2854 2855 E = E / GAMMA0 2856 EP = EP / GAMMA0 2857 EPP = EPP / GAMMA0 2858 EPPP = EPPP / GAMMA0 2859 EPPPP = EPPPP / GAMMA0 2860 2861 ERROR = ( DABS(ERR1) + DABS(ERR2) + DABS(ERR3) + DABS(ERR4) + 2862 & DABS(ERR5) + DABS(ERR6) + DABS(ERR7) + DABS(ERR8) + 2863 & DABS(ERR9) + DABS(ERR10)+ DABS(ERR11)+ DABS(ERR12)+ 2864 & DABS(ERR13)+ DABS(ERR14) ) / DABS(14.0d0*E) 2865 IF ( ERROR .GT. THRSDISP ) THEN 2866 WRITE (LUPRI,*) 2867 & 'WARNING: LARGE NUMERICAL (?) ERRORS DURING', 2868 & ' CALCULATION OF E COEFFICIENTS:' 2869 WRITE (LUPRI,*) 'TEST #1 GAVE:',ERR1 2870 WRITE (LUPRI,*) 'TEST #2 GAVE:',ERR2 2871 WRITE (LUPRI,*) 'TEST #3 GAVE:',ERR3 2872 WRITE (LUPRI,*) 'TEST #4 GAVE:',ERR4 2873 WRITE (LUPRI,*) 'TEST #5 GAVE:',ERR5 2874 WRITE (LUPRI,*) 'TEST #6 GAVE:',ERR6 2875 WRITE (LUPRI,*) 'TEST #7 GAVE:',ERR7 2876 WRITE (LUPRI,*) 'TEST #8 GAVE:',ERR8 2877 WRITE (LUPRI,*) 'TEST #9 GAVE:',ERR9 2878 WRITE (LUPRI,*) 'TEST #10 GAVE:',ERR10 2879 WRITE (LUPRI,*) 'TEST #11 GAVE:',ERR11 2880 WRITE (LUPRI,*) 'TEST #12 GAVE:',ERR12 2881 WRITE (LUPRI,*) 'TEST #13 GAVE:',ERR13 2882 WRITE (LUPRI,*) 'TEST #14 GAVE:',ERR14 2883 LERROR = .TRUE. 2884 END IF 2885 ABCDE(KE,1) = E 2886 ABCDE(KEP,1) = EP 2887 ABCDE(KEPP,1) = EPP 2888 ABCDE(KEPPP,1) = EPPP 2889 ABCDE(KEPPPP,1) = EPPPP 2890 IF (LOCDBG) THEN 2891 WRITE (LUPRI,*) " E, E', E'', E''', E'''' coefficients:", 2892 & E, EP, EPP, EPPP, EPPPP 2893 END IF 2894 END IF 2895 END IF 2896 2897 RETURN 2898 END 2899*---------------------------------------------------------------------* 2900* END OF SUBROUTINE CCCRDISP * 2901*---------------------------------------------------------------------* 2902c /* deck CCCDISPRT */ 2903*=====================================================================* 2904 SUBROUTINE CCCDISPRT (DISPCF, AVEDISP,ABCDE, 2905 & THGCF, ESHGCF, KERRCF, DFWMCF, 2906 & AVETHG, AVESHG, AVEKERR, AVEDFWM, LERROR) 2907*---------------------------------------------------------------------* 2908* 2909* Purpose: print dispersion coefficients second hyperpolarizabilities 2910* 2911* 2912* Written by Christof Haettig in March 1998. 2913* 2914*=====================================================================* 2915#if defined (IMPLICIT_NONE) 2916 IMPLICIT NONE 2917#else 2918# include "implicit.h" 2919#endif 2920#include "priunit.h" 2921#include "ccorb.h" 2922#include "ccsdinp.h" 2923#include "cccrinf.h" 2924#include "ccroper.h" 2925 2926 INTEGER KA, 2927 & KB, KBP, 2928 & KC, KCP, KCPP, 2929 & KD, KDP, KDPP, KDPPP, 2930 & KE, KEP, KEPP, KEPPP, KEPPPP 2931 PARAMETER (KA=1, 2932 & KB=2, KBP=3, 2933 & KC=4, KCP=5, KCPP=6, 2934 & KD=7, KDP=8, KDPP=9, KDPPP=10, 2935 & KE=11,KEP=12,KEPP=13,KEPPP=14,KEPPPP=15) 2936 2937 2938#if defined (SYS_CRAY) 2939 REAL DISPCF( (NCRDSPE+3)*(NCRDSPE+2)*(NCRDSPE+1)/6,*) 2940 REAL AVEDISP((NCRDSPE+3)*(NCRDSPE+2)*(NCRDSPE+1)/6,4) 2941 REAL ABCDE(15,2) 2942 REAL THGCF( NCRDSPE+1,NCROPER) 2943 REAL ESHGCF(NCRDSPE+1,NCROPER) 2944 REAL KERRCF(NCRDSPE+1,NCROPER) 2945 REAL DFWMCF(NCRDSPE+1,NCROPER) 2946 REAL AVETHG( NCRDSPE+1,4) 2947 REAL AVESHG( NCRDSPE+1,4) 2948 REAL AVEKERR(NCRDSPE+1,4) 2949 REAL AVEDFWM(NCRDSPE+1,4) 2950#else 2951 DOUBLE PRECISION DISPCF( (NCRDSPE+3)*(NCRDSPE+2)*(NCRDSPE+1)/6,*) 2952 DOUBLE PRECISION AVEDISP((NCRDSPE+3)*(NCRDSPE+2)*(NCRDSPE+1)/6,4) 2953 DOUBLE PRECISION ABCDE(15,2) 2954 DOUBLE PRECISION THGCF( NCRDSPE+1,NCROPER) 2955 DOUBLE PRECISION ESHGCF(NCRDSPE+1,NCROPER) 2956 DOUBLE PRECISION KERRCF(NCRDSPE+1,NCROPER) 2957 DOUBLE PRECISION DFWMCF(NCRDSPE+1,NCROPER) 2958 DOUBLE PRECISION AVETHG( NCRDSPE+1,4) 2959 DOUBLE PRECISION AVESHG( NCRDSPE+1,4) 2960 DOUBLE PRECISION AVEKERR(NCRDSPE+1,4) 2961 DOUBLE PRECISION AVEDFWM(NCRDSPE+1,4) 2962#endif 2963 2964 CHARACTER*5 BLANKS 2965 CHARACTER*80 STRING 2966 LOGICAL LERROR 2967 INTEGER ISYMA, ISYMB, ISYMC, ISYMD 2968 INTEGER IOPER, MN, LMN 2969 2970C 2971 INTEGER ITRI, L, M, N 2972 ITRI(L,M,N) = L*(L+1)*(L+2)/6 + M*(M+1)/2 + N + 1 2973C 2974 2975*---------------------------------------------------------------------* 2976* print header for hyperpolarizability section 2977*---------------------------------------------------------------------* 2978 BLANKS = ' ' 2979 STRING = ' RESULTS FOR DISPERSION COEFFICIENTS' 2980 2981 IF (CCS) THEN 2982 CALL AROUND( BLANKS//'FINAL CCS'//STRING(1:45)//BLANKS ) 2983 ELSE IF (CC2) THEN 2984 CALL AROUND( BLANKS//'FINAL CC2'//STRING(1:45)//BLANKS ) 2985 ELSE IF (CCSD) THEN 2986 CALL AROUND( BLANKS//'FINAL CCSD'//STRING(1:45)//BLANKS ) 2987 ELSE IF (CC3) THEN 2988 CALL AROUND( BLANKS//'FINAL CC3'//STRING(1:45)//BLANKS ) 2989 ELSE 2990 CALL QUIT('CCCDISPRT called for an unknown Coupled '// 2991 & 'Cluster model.') 2992 END IF 2993 2994 WRITE(LUPRI,'(/1X,4(A," op.",6X),4(4X,A),/,74("-"))') 2995 & 'A','B','C','D','L','M','N','D_ABCD' 2996 2997 2998 DO IOPER = 1, NCROPER 2999 ISYMA = ISYOPR(IACROP(IOPER)) 3000 ISYMB = ISYOPR(IBCROP(IOPER)) 3001 ISYMC = ISYOPR(ICCROP(IOPER)) 3002 ISYMD = ISYOPR(IDCROP(IOPER)) 3003 3004 IF (MULD2H(ISYMA,ISYMB).EQ.MULD2H(ISYMC,ISYMD)) THEN 3005 WRITE(LUPRI,'(/1X,4(A8,3X),3I5,1X,G16.8)') 3006 & LBLOPR(IACROP(IOPER)), LBLOPR(IBCROP(IOPER)), 3007 & LBLOPR(ICCROP(IOPER)), LBLOPR(IDCROP(IOPER)),0,0,0, 3008 & -DISPCF(ITRI(0,0,0),IOPER) 3009 ELSE 3010 IF (IPRCHYP.GT.0) THEN 3011 WRITE(LUPRI,'(/1X,4(A8,3X),4A)') LBLOPR(IACROP(IOPER)), 3012 & LBLOPR(IBCROP(IOPER)), LBLOPR(ICCROP(IOPER)), 3013 & LBLOPR(IDCROP(IOPER)),' - ',' - ', ' - ', '---' 3014 END IF 3015 END IF 3016 3017 IF (MULD2H(ISYMA,ISYMB).EQ.MULD2H(ISYMC,ISYMD)) THEN 3018 DO LMN = 2, NCRDSPE, 2 3019 DO MN = 0, LMN 3020 DO N = 0, MN 3021 L = LMN - MN 3022 M = MN - N 3023 WRITE(LUPRI,'(1X,4(8X,3X),3I5,1X,G16.8)') 3024 & L,M,N, -DISPCF(ITRI(LMN,MN,N),IOPER) 3025 END DO 3026 END DO 3027 END DO 3028 END IF 3029 3030 END DO 3031 3032 IF (GAMMA_PAR) THEN 3033 WRITE(LUPRI,'(/1X,A10,34X,3I5,1X,G16.8)') 3034 & 'gamma_{||} ',0,0,0, -AVEDISP(ITRI(0,0,0),1) 3035 DO LMN = 2, NCRDSPE, 2 3036 DO MN = 0, LMN 3037 DO N = 0, MN 3038 L = LMN - MN 3039 M = MN - N 3040 WRITE(LUPRI,'(6X,3(8X,5X),3I5,1X,G16.8)') 3041 & L,M,N, -AVEDISP(ITRI(LMN,MN,N),1) 3042 END DO 3043 END DO 3044 END DO 3045 END IF 3046 3047 IF (GAMMA_PAR .AND. GAMMA_ORT) THEN 3048 WRITE(LUPRI,'(/1X,A10,34X,3I5,1X,G16.8)') 'gamma_K ',0,0,0, 3049 & -1.5d0*(AVEDISP(ITRI(0,0,0),1)-AVEDISP(ITRI(0,0,0),2)) 3050 DO LMN = 2, NCRDSPE, 2 3051 DO MN = 0, LMN 3052 DO N = 0, MN 3053 L = LMN - MN 3054 M = MN - N 3055 WRITE(LUPRI,'(1X,10X,34X,3I5,1X,G16.8)') L,M,N, 3056 & -1.5d0*(AVEDISP(ITRI(LMN,MN,N),1)-AVEDISP(ITRI(LMN,MN,N),2)) 3057 END DO 3058 END DO 3059 END DO 3060 END IF 3061 3062 IF (GAMMA_ORT) THEN 3063 WRITE(LUPRI,'(/1X,A11,33X,3I5,1X,G16.8)') 3064 & 'gamma_{_|_}',0,0,0, -AVEDISP(ITRI(0,0,0),2) 3065 DO LMN = 2, NCRDSPE, 2 3066 DO MN = 0, LMN 3067 DO N = 0, MN 3068 L = LMN - MN 3069 M = MN - N 3070 WRITE(LUPRI,'(6X,3(8X,5X),3I5,1X,G16.8)') 3071 & L,N,M, -AVEDISP(ITRI(LMN,MN,N),2) 3072 END DO 3073 END DO 3074 END DO 3075 3076 WRITE(LUPRI,'(/1X,A10,34X,3I5,1X,G16.8)') 3077 & 'gamma_{ms} ',0,0,0, -AVEDISP(ITRI(0,0,0),3) 3078 DO LMN = 2, NCRDSPE, 2 3079 DO MN = 0, LMN 3080 DO N = 0, MN 3081 L = LMN - MN 3082 M = MN - N 3083 WRITE(LUPRI,'(6X,3(8X,5X),3I5,1X,G16.8)') 3084 & L,N,M, -AVEDISP(ITRI(LMN,MN,N),3) 3085 END DO 3086 END DO 3087 END DO 3088 END IF 3089 3090 WRITE(LUPRI,'(74("-"),//)') 3091 3092 WRITE(LUPRI,'(////1X,4(A," op.",6X),2(4X,A),/,112("-"))') 3093 & 'A','B','C','D','N', 3094 & 'D^THG D^ESHG D^DFWM D^KERR' 3095 3096 DO IOPER = 1, NCROPER 3097 ISYMA = ISYOPR(IACROP(IOPER)) 3098 ISYMB = ISYOPR(IBCROP(IOPER)) 3099 ISYMC = ISYOPR(ICCROP(IOPER)) 3100 ISYMD = ISYOPR(IDCROP(IOPER)) 3101 3102 IF (MULD2H(ISYMA,ISYMB).EQ.MULD2H(ISYMC,ISYMD)) THEN 3103 WRITE(LUPRI,'(/1X,4(A8,3X),I5,1X,4G16.8)') 3104 & LBLOPR(IACROP(IOPER)), 3105 & LBLOPR(IBCROP(IOPER)), 3106 & LBLOPR(ICCROP(IOPER)), 3107 & LBLOPR(IDCROP(IOPER)),0, 3108 & -THGCF(1,IOPER), -ESHGCF(1,IOPER), 3109 & -DFWMCF(1,IOPER), -KERRCF(1,IOPER) 3110 ELSE 3111 IF (IPRCHYP.GT.0) THEN 3112 WRITE(LUPRI,'(/1X,4(A8,3X),A,4(A,8X))') 3113 & LBLOPR(IACROP(IOPER)), 3114 & LBLOPR(IBCROP(IOPER)), 3115 & LBLOPR(ICCROP(IOPER)), 3116 & LBLOPR(IDCROP(IOPER)),' - ', 3117 & '---','---','---','---' 3118 END IF 3119 END IF 3120 3121 IF (MULD2H(ISYMA,ISYMB).EQ.MULD2H(ISYMC,ISYMD)) THEN 3122 DO N = 2, NCRDSPE, 2 3123 WRITE(LUPRI,'(1X,4(8X,3X),I5,1X,4G16.8)') N, 3124 & -THGCF(N+1,IOPER), -ESHGCF(N+1,IOPER), 3125 & -DFWMCF(N+1,IOPER),-KERRCF(N+1,IOPER) 3126 END DO 3127 END IF 3128 3129 END DO 3130 3131 IF (GAMMA_PAR) THEN 3132 WRITE(LUPRI,'(/1X,A10,34X,I5,1X,4G16.8)') 'gamma_{||} ',0, 3133 & -AVETHG(1,1),-AVESHG(1,1),-AVEDFWM(1,1),-AVEKERR(1,1) 3134 DO N = 2, NCRDSPE, 2 3135 WRITE(LUPRI,'(1X,44X,I5,1X,4G16.8)') N, 3136 & -AVETHG(N+1,1),-AVESHG(N+1,1),-AVEDFWM(N+1,1),-AVEKERR(N+1,1) 3137 END DO 3138 END IF 3139 IF (GAMMA_PAR .AND. GAMMA_ORT) THEN 3140 WRITE(LUPRI,'(/1X,A10,34X,I5,1X,4G16.8)') 'gamma^K ',0, 3141 & -1.5d0*(AVETHG(1,1) -AVETHG(1,2)), 3142 & -1.5d0*(AVESHG(1,1) -AVESHG(1,2)), 3143 & -1.5d0*(AVEDFWM(1,1)-AVEDFWM(1,2)), 3144 & -1.5d0*(AVEKERR(1,1)-AVEKERR(1,2)) 3145 DO N = 2, NCRDSPE, 2 3146 WRITE(LUPRI,'(1X,44X,I5,1X,4G16.8)') N, 3147 & -1.5d0*(AVETHG(N+1,1) -AVETHG(N+1,2)), 3148 & -1.5d0*(AVESHG(N+1,1) -AVESHG(N+1,2)), 3149 & -1.5d0*(AVEDFWM(N+1,1)-AVEDFWM(N+1,2)), 3150 & -1.5d0*(AVEKERR(N+1,1)-AVEKERR(N+1,2)) 3151 END DO 3152 END IF 3153 IF (GAMMA_ORT) THEN 3154 WRITE(LUPRI,'(/1X,A11,33X,I5,1X,4G16.8)') 'gamma_{_|_}',0, 3155 & -AVETHG(1,2),-AVESHG(1,2),-AVEDFWM(1,2),-AVEKERR(1,2) 3156 DO N = 2, NCRDSPE, 2 3157 WRITE(LUPRI,'(1X,44X,I5,1X,4G16.8)') N, 3158 & -AVETHG(N+1,2),-AVESHG(N+1,2),-AVEDFWM(N+1,2),-AVEKERR(N+1,2) 3159 END DO 3160 WRITE(LUPRI,'(/1X,A10,34X,I5,1X,4G16.7)') 'gamma_{ms} ',0, 3161 & -AVETHG(1,3),-AVESHG(1,3),-AVEDFWM(1,3),-AVEKERR(1,3) 3162 DO N = 2, NCRDSPE, 2 3163 WRITE(LUPRI,'(1X,44X,I5,1X,4G16.8)') N, 3164 & -AVETHG(N+1,3),-AVESHG(N+1,3),-AVEDFWM(N+1,3),-AVEKERR(N+1,3) 3165 END DO 3166 END IF 3167 3168 WRITE(LUPRI,'(112("-"),//)') 3169 3170 IF (GAMMA_PAR) THEN 3171 WRITE(LUPRI,'(////1X,A,/,78("-"))') 3172 & 'Process independent dispersion coefficients for gamma_{||}:' 3173 WRITE(LUPRI,'(1X,A,G16.8)') 'gamma_0',-AVEDISP(1,1) 3174 IF (NCRDSPE.GE.2) THEN 3175 WRITE(LUPRI,'(1X,A,2X,G16.8,5X)') "A ",ABCDE(KA,1) 3176 END IF 3177 IF (NCRDSPE.GE.4) THEN 3178 WRITE(LUPRI,'(1X,2(A,2X,G16.8,5X))') 3179 & "B ",ABCDE(KB,1), "B' ",ABCDE(KBP,1) 3180 END IF 3181 IF (NCRDSPE.GE.6) THEN 3182 WRITE(LUPRI,'(1X,3(A,2X,G16.8,5X))') 3183 & "C ",ABCDE(KC,1), "C' ",ABCDE(KCP,1), 3184 & "C'' ",ABCDE(KCPP,1) 3185 END IF 3186 IF (NCRDSPE.GE.8) THEN 3187 WRITE(LUPRI,'(1X,4(A,2X,G16.8,5X))') 3188 & "D ",ABCDE(KD,1), "D' ",ABCDE(KDP,1), 3189 & "D'' ",ABCDE(KDPP,1), "D''' ",ABCDE(KDPPP,1) 3190 END IF 3191 IF (NCRDSPE.GE.10) THEN 3192 WRITE(LUPRI,'(1X,5(A,2X,G16.8,5X))') 3193 & "E ",ABCDE(KE,1), "E' ",ABCDE(KEP,1), 3194 & "E'' ",ABCDE(KEPP,1), "E''' ",ABCDE(KEPPP,1), 3195 & "E''''",ABCDE(KEPPPP,1) 3196 END IF 3197 IF (NCRDSPE.GE.2 .AND. LERROR) THEN 3198 WRITE(LUPRI,'(/1X,A,/)') 'large numerical (?) errors ', 3199 & 'encountered ... check output for warnings.' 3200 END IF 3201 WRITE(LUPRI,'(78("-"))') 3202 END IF 3203 3204 IF (GAMMA_ORT) THEN 3205 WRITE(LUPRI,'(///1X,A,/,78("-"))') 3206 & 'Process independent dispersion coefficients for gamma_{ms}:' 3207 IF (NCRDSPE.GE.2) 3208 & WRITE(LUPRI,'(1X,A,2X,G16.8)') "A ",ABCDE(KA,2) 3209 IF (NCRDSPE.GE.4) THEN 3210 WRITE(LUPRI,'(1X,A,2X,G16.8,5X,A,2X,G16.8)') 3211 & "B ",ABCDE(KB,2), "B' ",ABCDE(KBP,2) 3212 END IF 3213 IF (NCRDSPE.GE.6) THEN 3214 WRITE(LUPRI,'(1X,3(A,2X,G16.8,5X))') "C ",ABCDE(KC,2), 3215 & "C' ",ABCDE(KCP,2),"C'' ",ABCDE(KCPP,2) 3216 END IF 3217 IF (NCRDSPE.GE.2 .AND. LERROR) THEN 3218 WRITE(LUPRI,'(/1X,A,/)') 'large numerical (?) errors ', 3219 & 'encountered ... check output for warnings.' 3220 END IF 3221 WRITE(LUPRI,'(78("-"))') 3222 END IF 3223 3224 WRITE(LUPRI,'(///)') 3225 3226 RETURN 3227 END 3228*---------------------------------------------------------------------* 3229* END OF SUBROUTINE CCCDISPRT * 3230*---------------------------------------------------------------------* 3231