1! 2! Dalton, a molecular electronic structure program 3! Copyright (C) by the authors of Dalton. 4! 5! This program is free software; you can redistribute it and/or 6! modify it under the terms of the GNU Lesser General Public 7! License version 2.1 as published by the Free Software Foundation. 8! 9! This program is distributed in the hope that it will be useful, 10! but WITHOUT ANY WARRANTY; without even the implied warranty of 11! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 12! Lesser General Public License for more details. 13! 14! If a copy of the GNU LGPL v2.1 was not distributed with this 15! code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html. 16! 17! 18C 19#ifdef revlog 20!=========================================================================== 21!921011-pj+hjaaj: EFPOLL added C6IFC test 22!920722-Hinne Hettema + HJAaJ 23!inserted Hinnes changes in EFPOLL,POLIM (DIFMAX test of convergence, 24!changed print levels, new format for LURSP7, fixed CFREQ bug) 25!=========================================================================== 26#endif 27C /* Deck rspc6 */ 28 SUBROUTINE RSPC6(CMO,UDV,PV,FC,FV,FCAC,H2AC,XINDX,WRK,LWRK) 29C 30#include "implicit.h" 31C 32 LOGICAL LDONE 33#include "iratdef.h" 34C 35 CHARACTER*8 BLANK 36 DIMENSION CMO(*),UDV(*),PV(*),FC(*),FV(*),FCAC(*),H2AC(*) 37 DIMENSION XINDX(*),WRK(LWRK) 38C 39#include "priunit.h" 40#include "infpri.h" 41#include "infrsp.h" 42#include "wrkrsp.h" 43#include "rspprp.h" 44#include "infc6.h" 45#include "inflr.h" 46#include "infdim.h" 47#include "inforb.h" 48C 49 PARAMETER ( ZERO = 0.0D0 , DM1 = -1.0D0 , BIG = 1.0D10 ) 50 PARAMETER ( DLRTST = 1.0D-6, BLANK = ' ' ) 51C 52C DETERMINE POLARIZABILITIES AT IMAGINARY FREQUENCIES. 53C 54C EXPAND V.(E[2]-OMEGA*S[2])**(-1).Vt = alpha AS CAUCHY-TYPE SERIES 55C AND FIND A BASIS OF PSEUDO-STATES BY RECURSIVE SOLUTION OF 56C LINEAR EQUATIONS ARISING AT EACH ORDER OF OMEGA. 57C FINALLY DIAGONALISE WTHIN REDUCED BASIS AND OBTAIN 16 58C POLARISABILITIES AS SUMS OF PDOSD CONTRIBUTIONS. 59C SIZE OF PROBLEM IS GOVERNED BY A CUTOFF ON THE NORM OF THE 60C PSEUDOSTATE VECTORS S[2]*LAMBDA[N] OR BY A PRESET MAXIMUM 61C NUMBER OF MOMENTS IN THE EXPANSION. 62C 63C NUMBER OF VECTORS OF SYMMETRY KSYMOP 64C 65 NGPNUM = NGPC6(KSYMOP) 66 NGRDMX = MAX(16,NGRID) 67C 68 KREDE = 1 69 KREDS = KREDE + MAXRM*MAXRM 70 KIBTYP = KREDS + MAXRM*MAXRM 71 KEIVAL = KIBTYP + MAXRM 72 KRESID = KEIVAL + MAXRM 73 KEIVEC = KRESID + MAXRM 74 KREDGP = KEIVEC + MAXRM*MAXRM 75 KGP = KREDGP + MAXRM 76 KSOL = KGP + KZYVAR*NGPNUM 77 KXMOM = KSOL + KZYVAR*NGPNUM 78 KOLDAL = KXMOM + MAXRM*NGPNUM 79 KWRK1 = KOLDAL + NGPNUM*(NGRDMX+1) 80 LWRK1 = LWRK - KWRK1 81C 82 LBVMAX = MAX(KZCONF,KZYWOP) 83 IF (LWRK1.LT.2*KZYVAR+LBVMAX) THEN 84 WRITE (LUPRI,9100) LWRK1,2*KZYVAR+LBVMAX 85 CALL QTRACE(LUPRI) 86 CALL QUIT('RSPC6 ERROR, INSUFFICIENT SPACE TO SOLVE LIN. EQ.S') 87 ENDIF 88 9100 FORMAT(/' RSPC6, work space too small for 2.5 (z,y)-vectors', 89 * /' had:',I10,', need more than:',I10) 90C 91C WORK SPACE FOR RSPEVE 92C 93 KBVECS = KWRK1 94 KWRKE = KBVECS + KZYVAR 95 LWRKE = LWRK - KWRKE 96 IF (LWRKE.LT.0) CALL ERRWRK('RSPC6',KWRKE-1,LWRK) 97C 98 KZRED = 0 99 KZYRED = 0 100 IPRRSP = IPRC6 101 MAXIT = MAXITC 102C 103 CALL DZERO(WRK(KOLDAL),NGPNUM*(NGRDMX+1)) 104C 105 IF ( IPRRSP.GT.83) THEN 106 WRITE(LUPRI,*)' THCC6,IPRC6,MAXITC',THCC6,IPRC6,MAXITC 107 WRITE(LUPRI,*)' NGPNUM',NGPNUM 108 WRITE(LUPRI,*)' MAXMOM',MAXMOM 109 END IF 110C 111C Call RSPCTL to solve linear set of response equations 112C 113 KGPOFF = KGP 114 KSOLOF = KSOL 115 DO 600 IOP = 1,NGPNUM 116 WRITE(LUPRI,'(//A)') ' Solving RSPC6 for' 117 WRITE(LUPRI,'(A,A/A,I4)') 118 * ' OPERATOR TYPE: ',LBLC6(KSYMOP,IOP), 119 * ' OPERATOR SYMMETRY:',KSYMOP 120 CALL GETGPV(LBLC6(KSYMOP,IOP),FC,FV,CMO,UDV,PV,XINDX,ANTSYM, 121 * WRK(KGPOFF),LWRK1) 122 KGPOFF = KGPOFF + KZYVAR 123 600 CONTINUE 124 NTOTAL = KZYVAR*NGPNUM 125 DO 610 I = 1,NTOTAL 126 WRK(KSOL-1+I) = WRK(KGP-1+I) 127 610 CONTINUE 128C CALL DCOPY(KZYVAR*NGPNUM,WRK(KGP),1,WRK(KSOL),1) 129 KEXSIM = 1 130 KEXCNV = 1 131 ITC6 = 0 132 C6NEW = BIG 133 C6DIF = BIG 134 THCRSP = THCC6 135 LDONE = .FALSE. 136 100 CONTINUE 137 WRITE(LUPRI,'(//A,I3/)') 138 & ' **** RSPC6 moment iteration no.',ITC6 139 KSOLOF = KSOL 140 KGPOFF = KGP 141 DO 200 IOP = 1,NGPNUM 142 WRITE(LUPRI,'(A,A/A,I4)') 143 * ' OPERATOR TYPE :',LBLC6(KSYMOP,IOP), 144 * ' OPERATOR SYMMETRY:',KSYMOP 145 WRK(KEIVAL) = ZERO 146 IF (IPRRSP.GT.110) THEN 147 WRITE(LUPRI,'(/A,I5)') 148 * ' GRADIENT VECTOR FOR CAUCHY SERIES NUMBER',ITC6 149 CALL OUTPUT(WRK(KSOLOF),1,2,1,KZVAR,2,KZVAR,-1,LUPRI) 150 END IF 151 XLRTST = DDOT(KZYVAR,WRK(KSOLOF),1,WRK(KSOLOF),1) 152 IF( SQRT(XLRTST).LE.DLRTST) THEN 153 WRITE(LUPRI,'(/A/A)') 154 * ' *** RIGHT HAND SIDE OF LINEAR EQUATION IS ZERO', 155 * ' *** THE LINEAR EQUATION IS NOT SOLVED' 156 GO TO 200 157 END IF 158 CALL RSPCTL(CMO,UDV,PV,FC,FV,FCAC,H2AC, 159 * .TRUE.,LBLC6(KSYMOP,IOP),BLANK,WRK(KSOLOF), 160 * WRK(KREDGP),WRK(KREDE),WRK(KREDS), 161 * WRK(KIBTYP),WRK(KEIVAL),WRK(KRESID),WRK(KEIVEC), 162 * XINDX,WRK(KWRK1),LWRK1) 163C CALL RSPCTL(CMO,UDV,PV,FC,FV,FCAC,H2AC, 164C LINEQ,GP,REDGP,REDE,REDS, 165C * IBTYP,EIVAL,EIVEC,XINDX,WRK,LWRK) 166C 167 RESTLR = .TRUE. 168 IF ( IPRRSP.GT.5) THEN 169 WRITE(LUPRI,'(/A,1P,G15.2)') 170 * ' Threshold for linear equations THCRSP:',THCRSP 171 END IF 172 CALL RSPEVE(WRK(KIBTYP),WRK(KEIVAL),WRK(KEIVEC), 173 * WRK(KBVECS),WRK(KWRKE),1,0) 174C CALL RSPEVE(IBTYP,EIVAL,EIVEC,BVECS,WRK,NBX,IBOFF) 175 IF (IPRRSP.GT.100) THEN 176 WRITE(LUPRI,'(/A,I5/2A)') 177 * ' SOLUTION VECTOR FOR lamda(k) k =',ITC6, 178 * ' OPERATOR ',LBLC6(KSYMOP,IOP) 179 CALL OUTPUT(WRK(KBVECS),1,2,1,KZVAR,2,KZVAR,-1,LUPRI) 180 END IF 181 GPNORM = DDOT(KZYVAR,WRK(KSOLOF),1,WRK(KBVECS),1) 182 IF (IPRRSP.GT.12) WRITE(LUPRI,'(/A,1P,G15.8)') 183 * ' GPNORM',GPNORM 184 IF (IPRRSP.GT.10) THEN 185 CAUCMO = DDOT(KZYVAR,WRK(KGPOFF),1,WRK(KBVECS),1) 186 WRITE(LUPRI,'(/A,I3,A,1P,G15.8)') 187 * ' CAUCHY MOMENT itc6 =',ITC6,' CAUCMO = ',CAUCMO 188 END IF 189C 190C CONSTRUCT S[2]*X WHERE X IS THE SOLUTION VECTOR 191C 192 CALL DZERO(WRK(KSOLOF),KZYVAR) 193 IF (KZCONF.GT.0) THEN 194 CALL RSPSLI(1,0,WRK(KBVECS+KZVAR),DUMMY, 195 * UDV,WRK(KSOLOF),XINDX,WRK(KWRKE),LWRKE) 196 CALL DSWAP(KZVAR,WRK(KSOLOF),1,WRK(KSOLOF+KZVAR),1) 197 CALL DSCAL(KZYVAR,DM1,WRK(KSOLOF),1) 198 CALL RSPSLI(1,0,WRK(KBVECS),DUMMY, 199 * UDV,WRK(KSOLOF),XINDX,WRK(KWRKE),LWRKE) 200 END IF 201 IF (KZWOPT.GT.0) THEN 202 DO 191 II = 1,KZWOPT 203 WRK(KBVECS+KZVAR-1+II) = 204 * WRK(KBVECS+KZVAR+KZCONF-1+II) 205 191 CONTINUE 206C CALL DCOPY(KZWOPT,WRK(KBVECS+KZVAR+KZCONF),1, 207C * WRK(KBVECS+KZVAR),1) 208 CALL RSPSLI(0,1,WRK(KBVECS+KZCONF),WRK(KBVECS+KZCONF), 209 * UDV,WRK(KSOLOF),XINDX,WRK(KWRKE),LWRKE) 210 END IF 211 KGPOFF = KGPOFF + KZYVAR 212 KSOLOF = KSOLOF + KZYVAR 213 200 CONTINUE 214 ITC6 = ITC6 + 1 215 IF ( ITC6.GT.MAXMOM ) THEN 216 WRITE(LUPRI,'(2A/A,I4,A,I4)') 217 * ' *** WARNING: Stopping ITC6 iterations, but', 218 * ' polarizabilities not converged.', 219 * ' ITC6 = ', ITC6, ' MAXMOM = ', MAXMOM 220 NWARN = NWARN + 1 221 GO TO 400 222 ENDIF 223 CALL POLIM(WRK(KIBTYP),WRK(KREDS),WRK(KREDE), 224 * WRK(KXMOM),WRK(KEIVAL),WRK(KEIVEC),WRK(KGP), 225 * WRK(KOLDAL),UDV,XINDX,WRK(KWRK1),LWRK1,LDONE) 226 IF (LDONE) THEN 227 WRITE(LUPRI,'(A,/A,I4,/A)') ' ', 228 * ' Done. Number of ITC6 iterations used :', ITC6, 229 * ' =============================================' 230 GOTO 400 231 ENDIF 232 GO TO 100 233C loop 100 finish 234 400 CONTINUE 235C 236 LDONE = .TRUE. 237 IF ( MOD(ITC6,2).EQ.0 ) ITC6 = ITC6 - 1 238 CALL POLIM(WRK(KIBTYP),WRK(KREDS),WRK(KREDE), 239 * WRK(KXMOM),WRK(KEIVAL),WRK(KEIVEC),WRK(KGP), 240 * WRK(KOLDAL),UDV,XINDX,WRK(KWRK1),LWRK1,LDONE) 241 242 RESTLR = .FALSE. 243C 244C *** end of RSPC6 -- 245C 246 RETURN 247 END 248C 249C /* Deck polim */ 250 SUBROUTINE POLIM(IBTYP,REDS,REDE,XMOM,EIVAL,EIVEC, 251 * GP,OLDAL,UDV,XINDX,WRK,LWRK,LDONE) 252C 253C 254C LDONE is a flag. LDONE=.false. THE QUADRATURE IS CARRIED OUT. 255C .true. THE CAUCHY PROCEDURE HAS 256C CONVERGED AND THE PSEUDOSTATES ARE USED 257C TO EVALUATE A SET OF REAL-FREQUENCY 258C POLARISABILITIES, OSCILLATOR SUMS 259C AND THE W4 COEFFICIENT. 260C 261C OUTPUT EIVAL(KZRED) AND XMOM(KZRED,NOP) 262C 263#include "implicit.h" 264C 265#include "priunit.h" 266#include "inforb.h" 267#include "infdim.h" 268#include "rspprp.h" 269#include "infc6.h" 270#include "infpri.h" 271#include "infrsp.h" 272#include "wrkrsp.h" 273C 274 LOGICAL LDONE 275 DIMENSION XMOM(MAXRM,*) 276 DIMENSION IBTYP(*),REDS(*),REDE(*),UDV(*),XINDX(*),WRK(LWRK) 277 DIMENSION EIVAL(*),EIVEC(*),GP(KZYVAR,*) 278 DIMENSION OLDAL(*) 279C 280 PARAMETER ( ZERO = 0.0D0 , D1 = 1.0D0 , 281 * D2 = 2.0D0 , D3 = 3.0D0, D10= 10.0D0, DUMMY = 1.0D20 ) 282 PARAMETER ( COMPLX = 1.0D7 ) 283C 284 IF (IPRRSP.GE.10) WRITE (LUPRI,'(//A,L2/)') 285 & ' *** Output from POLIM *** done =',LDONE 286 CALL PPRST(IBTYP,REDS,REDE,UDV,XINDX,WRK,LWRK) 287 CALL RSPRED(2,.FALSE.,0,IBTYP,DUMMY,DUMMY,REDE,REDS, 288 * EIVAL,EIVEC,WRK,WRK,UDV,WRK, 289 * XINDX,WRK,LWRK) 290C 291C CALCULATE EIGENVECTORS AND TRANSITION MOMENTS 292C 293C ALLOCATE WORK SPACE 294C 295C MAXIMUM NUMBER OF TRIAL VECTORS 296C 297 NSIMMA = (LWRK-KZYVAR)/KZYVAR 298 IF (NSIMMA.GT.KZRED) THEN 299 NSIM = KZRED 300 ELSE 301 NSIM = NSIMMA 302 END IF 303 IF (NSIM .LE. 0) THEN 304 WRITE (LUPRI,*) 'ERROR in POLIM, NSIM = 0' 305 WRITE (LUPRI,*) 'More memory required.' 306 CALL QUIT('Insufficient memory in POLIM') 307 END IF 308 KBVECS = 1 309 KWRK2 = KBVECS + NSIM*KZYVAR 310 LWRK2 = LWRK - KWRK2 311 IF (LWRK2.LT.0) CALL ERRWRK('POLIM',KWRK2-1,LWRK) 312C 313 DO 500 ISIM = 1,KZRED,NSIM 314 NBX = MIN( NSIM,(KZRED+1-ISIM) ) 315 CALL RSPEVE(IBTYP,EIVAL,EIVEC,WRK(KBVECS), 316 * WRK(KWRK2),NBX,(ISIM-1)) 317C CALL RSPEVE(IBTYP,EIVAL,EIVEC,BVECS,WRK,NBX,IBOFF) 318 DO 550 INUM = 1,NBX 319 DO 560 IOP = 1,NGPC6(KSYMOP) 320 XMOM(ISIM-1+INUM,IOP) = DDOT(KZYVAR,GP(1,IOP),1, 321 * WRK(KBVECS+(INUM-1)*KZYVAR),1) 322 560 CONTINUE 323 550 CONTINUE 324 500 CONTINUE 325 IF (IPRRSP.GE.30) THEN 326 DO 600 IOP = 1,NGPC6(KSYMOP) 327 WRITE(LUPRI,'(/A/2A/A,I4)') 328 * ' Test output of pseudo spectrum for', 329 * ' OPERATOR TYPE: ',LBLC6(KSYMOP,IOP), 330 * ' OPERATOR SYMMETRY:',KSYMOP 331 WRITE(LUPRI,'(/2A)') 332 * ' STATE NO: *TRANSITION MOMENT: ENERGY(au)' 333 DO 700 IST = 1,KZRED 334 WRITE(LUPRI,'(1X,I10,2F15.5)')IST,XMOM(IST,IOP),EIVAL(IST) 335 700 CONTINUE 336 600 CONTINUE 337 END IF 338C 339C Print of pseudo-spectrum 340C 341 IF (LDONE) THEN 342 LIMPRI = 5 343 ELSE 344 LIMPRI = 20 345 END IF 346 IF (IPRRSP .GE. LIMPRI) THEN 347 DO 800 IOP = 1,NGPC6(KSYMOP) 348 IF (LBLC6(KSYMOP,IOP)(2:4) .EQ. 'DIP') THEN 349 CALL C6PRPS(LDONE,LBLC6(KSYMOP,IOP),KSYMOP, 350 & KZRED,XMOM,EIVAL) 351 END IF 352 800 CONTINUE 353 END IF 354C 355 NVEL = 0 356 NCAR = 0 357 NSPH = 0 358 DO 900 IOP = 1,NGPC6(KSYMOP) 359 IF(LBLC6(KSYMOP,IOP)(2:5).EQ.'DIPL' ) NCAR = NCAR + 1 360 IF(LBLC6(KSYMOP,IOP)(3:6).EQ.'QUAD' ) NCAR = NCAR + 1 361 IF(LBLC6(KSYMOP,IOP)(4:6).EQ.'MOM' ) NCAR = NCAR + 1 362 IF(LBLC6(KSYMOP,IOP)(2:5).EQ.'DIPV' ) NVEL = NVEL + 1 363 IF(LBLC6(KSYMOP,IOP)(1:2).EQ.'SM' ) NSPH = NSPH + 1 364 900 CONTINUE 365 IF (IPRRSP.GT.15) THEN 366 WRITE(LUPRI,'(3(/A,I3))') 367 * ' NUMBER OF CARTESIAN Cn LENGTH OPERATORS : NCAR =',NCAR, 368 * ' NUMBER OF SPHERICAL Cn LENGTH OPERATORS : NSPH =',NSPH, 369 * ' NUMBER OF Cn VELOCITY OPERATORS : NVEL =',NVEL 370 END IF 371 NGRDMX = MAX(16,NGRID) 372 NMOMC = 20 373C -- NMOMC is the number of Cauchy moments computed 374C 375 KGRID = 1 376 KALPHA = KGRID + NGRDMX+1 377 KCAUCH = KALPHA + (NGRDMX+1) * NGPC6(KSYMOP) * NGPC6(KSYMOP) 378 KWRK1 = KCAUCH + (NMOMC+1) * NGPC6(KSYMOP) * NGPC6(KSYMOP) 379 LWRK1 = LWRK - KWRK1 380 IF (NCAR.EQ.NGPC6(KSYMOP).OR.(NSPH.EQ.NGPC6(KSYMOP))) THEN 381 CALL EFPOLL(LDONE,WRK(KGRID),WRK(KALPHA),WRK(KCAUCH),NGRDMX, 382 * NMOMC,NGPC6(KSYMOP),EIVAL,XMOM,OLDAL) 383 ELSE IF (NVEL.EQ.NGPC6(KSYMOP)) THEN 384 CALL EFPOLV(LDONE,WRK(KGRID),WRK(KALPHA),WRK(KCAUCH),NGRDMX, 385 * NMOMC,NGPC6(KSYMOP),EIVAL,XMOM) 386 ELSE 387 WRITE (LUPRI,'(//A,I4,A,I4,A,I4)') 388 & 'ERROR: BOTH DIPLEN AND DIPVEL OPERATORS SPECIFIED. '// 389 & ' NCAR=',NCAR,' NVEL=',NVEL, 'NSPH= ',NSPH 390 CALL QUIT(' POLIM: Illegal option, operator types are mixed') 391 END IF 392 RETURN 393 END 394C /* Deck efpoll */ 395 SUBROUTINE EFPOLL(LDONE,GRIDSQ,ALPHA,CAUCHY,NGRDMX,NMOMC,NSYMOP, 396 * EIVAL,XMOM,OLDAL) 397c------------------------------------------------------------------------ 398c This subroutine calculates the frequency dependent polarisabilities on 399c a precomputed grid and optionally on real frequencies. 400c The default is the Gauss Chebychev grid, the alternative is the Gauss 401c Legendre grid. Details can be found in: 402c default: W.Rijks and P.E.S.Wormer JCP Vol.88, 5704 (1988). 403c G-Legendre: Amos et al, J.Phys.Chem. Vol.89, 2186 (1985) 404c 405c We predefine 10 Cauchy moments. On change also change the memory layout 406c in the C8 module. Now NMOMC input parameter (900718/hjaaj). 407c 408c Note: loops to NGRID go from 0(!) to NGRID, because in ALPHA(I,J,0) 409c we keep the static polarisability. We do not need this in the computa- 410c tion of the Cn coefficients. 411c 412c When this subroutine is called, we assume that we only have operators 413c of one type present, so that the interface file to the coupling program 414c is consistent. This is ensured by the IF THEN ELSE statement around the 415c call. 416c 417c We require too much space for the polarisability, but it facilitates 418c addressing and allows at the same time for properties which are not 419c symmetric in the indices. 420c------------------------------------------------------------------------ 421c 422#include "implicit.h" 423c 424#include "pi.h" 425#include "dummy.h" 426c 427#include "priunit.h" 428#include "infrsp.h" 429#include "rspprp.h" 430#include "infc6.h" 431#include "wrkrsp.h" 432C 433#include "inforb.h" 434#include "infpri.h" 435c 436 DIMENSION EIVAL(*),XMOM(MAXRM,*) 437 DIMENSION GRIDSQ(0:NGRDMX),ALPHA(1:NSYMOP,1:NSYMOP,0:NGRDMX) 438 DIMENSION CAUCHY(1:NSYMOP,1:NSYMOP,1:NMOMC) 439 DIMENSION OLDAL(NSYMOP,0:NGRDMX) 440 DIMENSION OMEGA(0:16) 441 CHARACTER*8 IC,JC 442 LOGICAL LDONE,END 443 LOGICAL FIRST, LC6IFC 444 DATA FIRST /.TRUE./ 445 SAVE FIRST 446 DATA ITYPE /3/ 447 DATA THRESHW /1.D-08/ 448c 449 DATA OMEGA/0.0D0, 450 * 0.0000011354D0,0.0000324954D0,0.0002074939D0, 451 * 0.0007766098D0,0.0022314002D0,0.0055272185D0, 452 * 0.0125684274D0,0.0273216537D0,0.0585616091D0, 453 * 0.1273031181D0,0.2894765239D0,0.7170385627D0, 454 * 2.0602366551D0,7.7110713698D0,49.2377677794D0, 455 * 1409.1898779477D0/ 456c -- parameters 457 PARAMETER ( ZERO = 0.0D0 , D1 = 1.0D0, D2 = 2.0D0 ) 458c 459c compute squared gridpoints 460c 461 IF(GSLEGN) THEN 462c -- use Gauss Legendre grid from old C6 module 463 NGRID=16 464 DO 5 I=0,NGRID 465 GRIDSQ(I) = OMEGA(I) 466 5 CONTINUE 467 ELSE 468c-- compute gridpoints according to Gauss Chebychev scheme 469 GRIDSQ(0) = ZERO 470 N= NGRID 471 NN=2*N 472 FAC=PI/(2*NN) 473 DO 15 I=1,N 474 X=FAC*(2*I-1) 475 GRIDSQ(N-I+1)=(D1/TAN(X))**2 476 15 CONTINUE 477 ENDIF 478c 479c write header output file when first call 480c 481 LURSP7 = -1 482 CALL GPOPEN(LURSP7,'RESPONSE.C8','UNKNOWN',' ','FORMATTED',IDUMMY, 483 & .FALSE.) 484 IF ( FIRST .AND. C6IFC ) THEN 485 WRITE(LURSP7,'(A)') '* SIRIUS OUTPUT FILE' 486 WRITE(LURSP7,'(A,I8)') '* NOCC',NOCCT 487 WRITE(LURSP7,'(A,I8)') '* NACT',NASHT 488 WRITE(LURSP7,'(A,I8)') '* NVIR',NSSHT 489 WRITE(LURSP7,'(A)') '* GROUP NONE' 490 WRITE(LURSP7,'(A,I8)') '* NOMEGA',NGRID 491 DO I=0,NGRID+1,5 492 JMAX = MIN(I+4,NGRID) 493 WRITE(LURSP7,'(A,5F12.6)') 494 & '*',(SQRT(GRIDSQ(J)),J=I,JMAX) 495 ENDDO 496 FIRST = .FALSE. 497 ENDIF 498c 499c initialise alpha and cauchy moments to zero 500c 501 CALL DZERO(ALPHA,NSYMOP*NSYMOP*(NGRID+1)) 502 CALL DZERO(CAUCHY,NSYMOP*NSYMOP*NMOMC) 503c 504 END = LDONE 505 LDONE = .TRUE. 506c 507c Compute the polarisability and cauchy moments from the effective spectrum. 508c We only need the lower triangle, hence we have j<= i in the second loop 509c 510c We use that the LBLC6 has the following format: "SM02-02" for 511c the l=2,m=-2 multipole matrix to extract the l,m 512c 513 IF (IPRRSP.GE.1) WRITE(LUPRI,'(/1X,A,/1X,A)') 514 & 'Static polarisabilities (length)', 515 & '------------------------------------------------' 516 DO 200 I=1,NSYMOP 517 IC=LBLC6(KSYMOP,I) 518 IF (C6IFC .AND. IC(1:2) .EQ. 'SM') THEN 519 READ(IC,'(2X,I2,I3)') L,M 520 LC6IFC = .TRUE. 521 ELSE 522 LC6IFC = .FALSE. 523 END IF 524 DO 210 J=1,I 525 JC=LBLC6(KSYMOP,J) 526 DO 100 IEFF=1,KZRED 527 OSCSTR = D2 * (XMOM(IEFF,I)*XMOM(IEFF,J)*EIVAL(IEFF)) 528 DO 110 NG=0,NGRID 529 ALPHA(I,J,NG) = ALPHA(I,J,NG) 530 & + OSCSTR/(EIVAL(IEFF)**2 + GRIDSQ(NG)) 531 110 CONTINUE 532 DO 111 KK=1,NMOMC 533 CAUCHY(I,J,KK) = CAUCHY(I,J,KK) 534 & + OSCSTR/(EIVAL(IEFF)**(2*KK) ) 535 111 CONTINUE 536 100 CONTINUE 537 IF (END .AND. LC6IFC) THEN 538 IF (JC(1:2) .EQ. 'SM' .AND. 539 & ALPHA(I,J,0) .GT. THRESHW ) THEN 540 READ(JC,'(2X,I2,I3)') LP,MP 541 DO, K=0, NGRID 542 WRITE(LURSP7,'(4I3,I5,1P,D17.8,I6)') 543 & L,LP,M,MP,K,ALPHA(I,J,K),ITYPE 544 END DO 545 ENDIF 546 ENDIF 547 IF (IPRRSP.GE.1) THEN 548 IF (I.EQ.J) WRITE(LUPRI,'(A)') ' ' 549 WRITE(LUPRI,'(1X,A8,1X,A8,F20.10)') IC,JC,ALPHA(I,J,0) 550 END IF 551 IF (END.AND.IPRRSP.GE.5) THEN 552 WRITE(LUPRI,*) ' GRIDSQ ALPHA' 553 DO, II=0,NGRID 554 WRITE(LUPRI,'(F15.7,F20.7)')GRIDSQ(II),ALPHA(I,J,II) 555 END DO 556 ENDIF 557c 558c compare old and new polarizabilities and decide on convergence 559c 560 IF (I.EQ.J .AND. .NOT.END) THEN 561 DIFFMX = ZERO 562 DO, II=0,NGRID 563 DIFF = ABS(OLDAL(I,II)-ALPHA(I,I,II)) 564 DIFFMX = MAX(DIFF,DIFFMX) 565 LDONE=(LDONE.AND.(DIFF.LT.THCC6)) 566 OLDAL(I,II) = ALPHA(I,I,II) 567 END DO 568 IF (C6IFC) WRITE(LURSP7,'(A,F14.7,3X,A)') 569 * '* DIFFERENCE FOUND: ', DIFFMX, LBLC6(KSYMOP,I) 570 IF (IPRRSP .GT. 1) THEN 571 WRITE(LUPRI,'(A,F14.7,3X,A/)') 572 * ' Max difference in grid polarizabilities: ', 573 * DIFFMX, LBLC6(KSYMOP,I) 574 ENDIF 575 ENDIF 576c 577 210 CONTINUE 578 200 CONTINUE 579 CALL GPCLOSE(LURSP7,'KEEP') 580c 581c Now do something very much the same to print the Cauchy moments 582c 583 IF ((END.AND.IPRRSP.GT.1) .OR. IPRRSP .GT. 10) THEN 584 WRITE(LUPRI,'(//A/A)') ' Cauchy moments (length)', 585 & ' ------------------------------------------------' 586 DO, I=1,NSYMOP 587 IC=LBLC6(KSYMOP,I) 588 DO, J=1,I 589 JC=LBLC6(KSYMOP,J) 590 WRITE(LUPRI,'(//1X,A8,1X,A8)') IC,JC 591 WRITE(LUPRI,'(/1X,A4,5X,A25)') ' k ','Moment = S(-k-2)' 592 DO, KK=1,NMOMC 593 WRITE(LUPRI,'(I5,1P,E33.20)') 2*KK-2,CAUCHY(I,J,KK) 594 END DO 595 END DO 596 END DO 597 END IF 598c 599c Now compute the polarisabilities at real frequencies 600c 601 IF((NCFREQ.GT.0).AND. (END.OR.IPRRSP.GT.10)) THEN 602 WRITE(LUPRI,'(/A/A/A,F10.4,A)') 603 & ' Polarisabilities at real frequencies (length)', 604 & ' ---------------------------------------------------', 605 & ' (lowest pole in effective spectrum:',EIVAL(1),'au )' 606 DO 400 I=1,NSYMOP 607 IC=LBLC6(KSYMOP,I) 608 DO 410 J=1,I 609 JC=LBLC6(KSYMOP,J) 610 WRITE(LUPRI,'(/1X,A8,2X,A8)') IC,JC 611 DO 420 NG=1,NCFREQ 612 POLRL = ZERO 613 DO, IEFF=1,KZRED 614 POLRL = POLRL 615 & + D2 * (XMOM(IEFF,I)*XMOM(IEFF,J)*EIVAL(IEFF)) 616 & /(EIVAL(IEFF)**2 - CFREQ(NG)**2) 617 END DO 618 WRITE(LUPRI,'(1X,2F20.10)') CFREQ(NG),POLRL 619 420 CONTINUE 620 410 CONTINUE 621 400 CONTINUE 622 END IF 623c 624 RETURN 625 END 626C /* Deck efpolv */ 627 SUBROUTINE EFPOLV(LDONE,GRIDSQ,ALPHA,CAUCHY,NGRDMX,NMOMC,NSYMOP, 628 * EIVAL,XMOM) 629c----------------------------------------------------------------------- 630c This subroutine calculates the frequency dependent polarisabilities in 631c the velocity representation. For details on grid etc. see subroutine 632c EFPOLL 633c The Cauchy moments begin with S(0) here, instead of S(-2) as in EFPOLL 634c----------------------------------------------------------------------- 635c 636#include "implicit.h" 637c 638#include "pi.h" 639c 640#include "priunit.h" 641#include "infrsp.h" 642#include "rspprp.h" 643#include "infc6.h" 644#include "wrkrsp.h" 645#include "infpri.h" 646c 647 DIMENSION EIVAL(*),XMOM(MAXRM,*) 648 DIMENSION GRIDSQ(0:NGRDMX),ALPHA(1:NSYMOP,1:NSYMOP,0:NGRDMX) 649 DIMENSION CAUCHY(1:NSYMOP,1:NSYMOP,0:NMOMC) 650 DIMENSION OMEGA(0:16) 651 CHARACTER*4 IC,JC 652 LOGICAL LDONE 653c 654 DATA OMEGA/0.0D0, 655 * 0.0000011354D0,0.0000324954D0,0.0002074939D0, 656 * 0.0007766098D0,0.0022314002D0,0.0055272185D0, 657 * 0.0125684274D0,0.0273216537D0,0.0585616091D0, 658 * 0.1273031181D0,0.2894765239D0,0.7170385627D0, 659 * 2.0602366551D0,7.7110713698D0,49.2377677794D0, 660 * 1409.1898779477D0/ 661c -- parameters 662 PARAMETER ( ZERO = 0.0D0 , D1 = 1.0D0, D2 = 2.0D0 ) 663c 664c compute squared gridpoints 665c 666 IF(GSLEGN) THEN 667c -- use Gauss Legendre grid from old C6 module 668 NGRID=16 669 DO, I=0,NGRID 670 GRIDSQ(I) = OMEGA(I) 671 END DO 672 ELSE 673c-- compute gridpoints according to Gauss Chebychev scheme 674 GRIDSQ(0) = ZERO 675 N= NGRID 676 NN=2*N 677 FAC=PI/(2*NN) 678 DO, I=1,N 679 X=FAC*(2*I-1) 680 GRIDSQ(N-I+1)=(D1/TAN(X))**2 681 END DO 682 ENDIF 683c 684c initialise alpha and cauchy moments to zero 685c 686 CALL DZERO(ALPHA ,NSYMOP*NSYMOP*(NGRID+1)) 687 CALL DZERO(CAUCHY,NSYMOP*NSYMOP*(NMOMC+1)) 688c 689c Compute the polarisability and cauchy moments from the effective spectrum. 690c We only need the lower triangle, hence we have j<= i in the second loop 691c 692 IF (IPRRSP .GE. 1) WRITE(LUPRI,'(//1X,A/,1X,A)') 693 & 'S(0) and static polarisabilities (velocity)', 694 & '------------------------------------------------' 695 IC = ' ' 696 JC = ' ' 697 DO 200 I=1,NSYMOP 698 IF(LBLC6(KSYMOP,I)(2:5).EQ.'DIPV' ) IC=LBLC6(KSYMOP,I)(1:1) 699 DO 210 J=1,I 700 IF(LBLC6(KSYMOP,J)(2:5).EQ.'DIPV' ) JC=LBLC6(KSYMOP,J)(1:1) 701 DO 100 IEFF=1,KZRED 702 OSCSTR = D2 * (XMOM(IEFF,I)*XMOM(IEFF,J)/EIVAL(IEFF)) 703 DO, NG=0,NGRID 704 ALPHA(I,J,NG) = ALPHA(I,J,NG) 705 & + OSCSTR/(EIVAL(IEFF)**2 + GRIDSQ(NG)) 706 END DO 707 CAUCHY(I,J,0) = CAUCHY(I,J,0) + OSCSTR 708 DO, KK=1,NMOMC 709 CAUCHY(I,J,KK) = CAUCHY(I,J,KK) 710 & + OSCSTR/(EIVAL(IEFF)**(2*KK) ) 711 END DO 712 100 CONTINUE 713 IF (IPRRSP.GE.1) THEN 714 WRITE(LUPRI,'(1X,A4,1X,A4,F20.10)') IC,JC,CAUCHY(I,J,0) 715 WRITE(LUPRI,'(1X,A4,1X,A4,F20.10)') IC,JC,ALPHA(I,J,0) 716 END IF 717 IF(LDONE .AND. IPRRSP.GE.10) THEN 718 WRITE(LUPRI,*) ' GRIDSQ ALPHA' 719 DO, II=0,NGRID 720 WRITE(LUPRI,'(1X,F14.7,F20.7)')GRIDSQ(II),ALPHA(I,J,II) 721 END DO 722 ENDIF 723 210 CONTINUE 724 200 CONTINUE 725c 726c Now do something very much the same to print the Cauchy moments 727c 728 IF ((LDONE.AND.IPRRSP.GT.1) .OR. IPRRSP .GT. 10) THEN 729 WRITE(LUPRI,'(//1X,A/,1X,A)') 'Cauchy moments (velocity)', 730 & '------------------------------------------------' 731 IC = ' ' 732 JC = ' ' 733 DO 300 I=1,NSYMOP 734 IF(LBLC6(KSYMOP,I)(2:5).EQ.'DIPV' ) IC=LBLC6(KSYMOP,I)(1:1) 735 DO 310 J=1,I 736 IF(LBLC6(KSYMOP,J)(2:5).EQ.'DIPV' ) JC=LBLC6(KSYMOP,J)(1:1) 737 WRITE(LUPRI,'(//1X,A4,1X,A4)') IC,JC 738 WRITE(LUPRI,'(/1X,A4,5X,A25)') ' k ','Moment = S(-k)' 739 DO, KK=0,NMOMC 740 WRITE(LUPRI,'(1X,I4,5X,1P,E30.20)') 2*KK,CAUCHY(I,J,KK) 741 END DO 742 310 CONTINUE 743 300 CONTINUE 744 END IF 745c 746c Now compute the polarisabilities at real frequencies 747c 748 IF((NCFREQ.GT.0).AND. (LDONE.OR.IPRRSP.GT.10)) THEN 749 WRITE(LUPRI,'(/1X,A/,1X,A/1X,A,F10.4,A)') 750 & 'Polarisabilities at real frequencies (velocity)', 751 & '---------------------------------------------------', 752 & '(lowest pole in effective spectrum:',EIVAL(1),'au )' 753 DO 400 I=1,NSYMOP 754 IF(LBLC6(KSYMOP,I)(2:5).EQ.'DIPV' ) IC=LBLC6(KSYMOP,I)(1:1) 755 DO 410 J=1,I 756 IF(LBLC6(KSYMOP,J)(2:5).EQ.'DIPV')JC=LBLC6(KSYMOP,J)(1:1) 757 WRITE(LUPRI,'(/1X,A4,3X,A4)') IC,JC 758 DO 420 NG=1,NCFREQ 759 POLRL = ZERO 760 DO, IEFF=1,KZRED 761 POLRL = POLRL 762 & + D2 * (XMOM(IEFF,I)*XMOM(IEFF,J)) 763 & /(EIVAL(IEFF)*(EIVAL(IEFF)**2 - CFREQ(NG)**2)) 764 END DO 765 WRITE(LUPRI,'(1X,2F20.10)') CFREQ(NG),POLRL 766 420 CONTINUE 767 410 CONTINUE 768 400 CONTINUE 769 END IF 770c 771 RETURN 772 END 773C /* Deck c6prps */ 774 SUBROUTINE C6PRPS(LDONE,OPTYPE,KSYMOP,KZRED,XMOM,EIVAL) 775C 776C 25-Jun-1990 PJ+HJAAJ 777C 778C Print pseudo-spectrum 779C 780#include "implicit.h" 781#include "codata.h" 782 LOGICAL LDONE 783 CHARACTER*8 OPTYPE 784 DIMENSION XMOM(KZRED),EIVAL(KZRED) 785 PARAMETER ( NSUM = 15 ) 786 DIMENSION SUMOSC(NSUM) 787 PARAMETER ( D0 = 0.0D0, D2 = 2.0D0 ) 788C 789#include "priunit.h" 790#include "infpri.h" 791C 792 IF (LDONE) THEN 793 WRITE(LUPRI,'(//A)') ' Final pseudo-spectrum for' 794 ELSE 795 WRITE(LUPRI,'(//A)') ' Intermediate pseudo-spectrum for' 796 END IF 797 WRITE(LUPRI,'(/A,A/A,I4)') 798 * ' operator type: ',OPTYPE, 799 * ' operator symmetry:',KSYMOP 800 IF (OPTYPE(2:4) .NE. 'DIP') THEN 801 NWARN = NWARN + 1 802 WRITE (LUPRI,'(//3A/A)') 803 & ' C6PRPS WARNING: oscillator strength is not', 804 & ' defined for ',OPTYPE, 805 & ' WARNING: C6PRPS cannot continue.' 806 RETURN 807 END IF 808 CALL DZERO(SUMOSC,NSUM) 809 WRITE(LUPRI,'(/A)') 810 * ' State no: Transition moment:'// 811 * ' 3*Oscillator strength: Energy (eV):' 812C:State no: Transition moment: 3*Oscillator strength: Energy (eV): 813C: 5 0.123456789012 0.123456789012 0.123456789012 814C(I5,1P,G26.12,G24.12,G20.12) 815 DO 500 INUM = 1,KZRED 816 XTEMP=XMOM(INUM) 817 ETEMP=EIVAL(INUM) 818 IF (OPTYPE(2:5) .EQ. 'DIPL') THEN 819 OSCSTR=D2*XTEMP*XTEMP*ETEMP 820 ELSE IF (OPTYPE(2:5) .EQ. 'DIPV') THEN 821 OSCSTR=D2*XTEMP*XTEMP/ETEMP 822 ELSE 823 OSCSTR=D0 824 END IF 825 WRITE(LUPRI,'(I5,1P,G26.12,G24.12,G20.12)') 826 * INUM,XTEMP,OSCSTR, ETEMP*XTEV 827 DO 400 ISUM = 1,NSUM 828 SUMOSC(ISUM) = SUMOSC(ISUM) 829 * + D2*XTEMP*XTEMP*ETEMP**(-ISUM) 830 400 CONTINUE 831 500 CONTINUE 832 IF (OPTYPE(2:5) .EQ. 'DIPL') THEN 833 WRITE (LUPRI,'(/A)') ' Sum rules (length)' 834 WRITE (LUPRI,'(/A)') ' k S(k)' 835 WRITE (LUPRI,'(I5,1P,G20.12)') 836 * ( (-ISUM-1), SUMOSC(ISUM), ISUM = 1,NSUM) 837 ELSE IF (OPTYPE(2:5) .EQ. 'DIPV') THEN 838 WRITE (LUPRI,'(/A)') ' Sum rules (velocity)' 839 WRITE (LUPRI,'(/A)') ' k S(k)' 840 WRITE (LUPRI,'(I5,1P,G20.12)') 841 * ( (-ISUM+1), SUMOSC(ISUM), ISUM = 1,NSUM) 842 END IF 843 RETURN 844 END 845C /* Deck paddy */ 846#if defined (VAR_PADDY) 847C920722: Old C6 code by Patrick Fowler 848 SUBROUTINE PADDY(NCAUCH,NI,CAUCHY,COMI,WRK ) 849#include "implicit.h" 850#include "priunit.h" 851#include "infpri.h" 852#include "infrsp.h" 853 PARAMETER ( D0 = 0.0D0 , D1 = 1.0D0 ) 854C 855 DIMENSION PLVAL(16) 856 DIMENSION CAUCHY(*),COMI(*),WRK(*) 857C 858C DETERMINE DIMENSION OF LINEAR SET OF EQUATIONS 859C 860 IF (MOD(NCAUCH,2).GT.0) THEN 861 K = 0 862 N = (NCAUCH-1)/2 863 M = N 864 ELSE 865 K = -1 866 N = NCAUCH/2 867 M = N - 1 868 END IF 869C 870 KCAUNE = 1 871 KCVEC = KCAUNE + NCAUCH 872 KCMAT = KCVEC + N 873 KBVEC = KCMAT + N*N 874 KAVEC = KBVEC + N +1 875 KWRK1 = KAVEC + M + 1 876C 877 DO 50 I = 1,NCAUCH 878 WRK(KCAUNE-1+I) = CAUCHY(I) 879 50 CONTINUE 880 DO 100 I = 1,N 881 WRK(KCVEC-1+I) = -WRK(KCAUNE+N+K+I) 882 DO 200 J = 1,N 883 WRK(KCMAT-1+(J-1)*N+I) = WRK(KCAUNE+N+K+I-J) 884 200 CONTINUE 885 100 CONTINUE 886C 887 CALL INVMAT(WRK(KCMAT),WRK(KWRK1),N,N) 888C 889 CALL DZERO(WRK(KBVEC),N+1) 890 CALL DGEMM('N','N',N,1,N,1.D0, 891 & WRK(KCMAT),N, 892 & WRK(KCVEC),N,1.D0, 893 & WRK(KBVEC+1),N) 894 WRK(KBVEC) = D1 895 DO 300 I = 1,N+K+1 896 WRK(KAVEC-1+I) = D0 897 DO 400 J = 1,I 898 WRK(KAVEC-1+I) = WRK(KAVEC-1+I) + WRK(KBVEC-1+J) 899 * *WRK(KCAUNE+I-J) 900 400 CONTINUE 901 300 CONTINUE 902 WRITE(LUPRI,'(/A,I2,A,I2,A,I3,A/)') 903 *' [',N,',',M,'] Pade approximant for ',NCAUCH,' moments' 904 IF (IPRRSP.GT.5) THEN 905 WRITE(LUPRI,'(/A)') 906 * ' M+1 coefficients in P numerator' 907 CALL OUTPUT(WRK(KAVEC),1,M+1,1,1,M+1,1,-1,LUPRI) 908 WRITE(LUPRI,'(/A)') 909 * ' N+1 coefficients in Q denominator' 910 CALL OUTPUT(WRK(KBVEC),1,N+1,1,1,N+1,1,-1,LUPRI) 911 END IF 912 DO 500 I = 1,NI 913 CALL POLVAL(M+1,WRK(KAVEC),COMI(I),VALP) 914 CALL POLVAL(N+1,WRK(KBVEC),COMI(I),VALQ) 915 PLVAL(I) = VALP/VALQ 916 WRITE(LUPRI,'(I5,A,1P,G16.8,A,G16.8)') 917 * I, ' alpha at imaginary frequency',SQRT(COMI(I)),' = ' 918 * ,PLVAL(I) 919 500 CONTINUE 920 CALL C6INT(PLVAL,PLVAL,C6TEMP) 921 WRITE(LUPRI,'(/A,I2,A,I2,A,1P,G16.8) ') 922 *' C6 coefficient for [',N,',',M,'] Pade approximant = ', 923 *C6TEMP 924 RETURN 925 END 926 SUBROUTINE POLVAL(NUM,POLCOF,FREQ,VALUE) 927#include "implicit.h" 928 DIMENSION POLCOF(NUM) 929 VALUE = 0.0D0 930 DO 100 I = 1,NUM 931 VALUE = VALUE + POLCOF(I)*((-FREQ)**(I-1)) 932 100 CONTINUE 933 RETURN 934 END 935 SUBROUTINE C6INT(A1,A2,SUM) 936#include "implicit.h" 937C 938#include "pi.h" 939C 940 DIMENSION RLOW(8),WLOW(8) 941 DIMENSION W(16),A1(16),A2(16) 942 DATA RLOW/ 943 &0.97891421016235D+00,0.89222197421380D+00, 944 &0.74931737854740D+00,0.57063582016217D+00,0.38177105339712D+00, 945 &0.20977936861551D+00,0.79300559811486D-01,0.90273770256471D-02/ 946 DATA WLOW/ 947 &0.27152459411755D-01,0.62253523938648D-01, 948 &0.95158511682493D-01,0.12462897125553D+00,0.14959598881658D+00, 949 &0.16915651939500D+00,0.18260341504492D+00,0.18945061045507D+00/ 950C 951C THIS PROGRAM CALCULATES THE C6 COEFFICIENT FOR THE INTERACTION 952C OF SYSTEM 1 WITH SYSTEM 2, GIVEN THE POLARISABILITY OF EACH 953C SYSTEM AT A SET OF 16 PRE-DEFINED IMAGINARY FREQUENCIES: 954C 955C C6 = (3/PI) INTEGRAL(0,INFINITY) ALPHA1(IW) ALPHA2(IW) DW 956C 957 OMEGA0=0.2D0 958C 959C FOR EACH PARTNER READ THE POLARISABILITY AT THE 16 FREQUENCIES. 960C THE FREQUENCIES REQUIRED ARE LISTED AS THE NEGATIVES OF THEIR 961C SQUARES BELOW. 962C 963C 0.0000011354 964C 0.0000324954 965C 0.0002074939 966C 0.0007766098 967C 0.0022314002 968C 0.0055272185 969C 0.0125684274 970C 0.0273216537 971C 0.0585616091 972C 0.1273031181 973C 0.2894765239 974C 0.7170385627 975C 2.0602366551 976C 7.7110713698 977C 49.2377677794 978C 1409.1898779475 979C 980C FOR A DESCRIPTION OF THE QUADRATURE SCHEME, SEE AMOS ET AL. 981C J. PHYS. CHEM. 1985 VOL. 89 PAGE 2186. 982C 983 DO 1 I=1,8 984 WW=WLOW(I)*OMEGA0/PI 985 RR=SQRT(RLOW(I)) 986 I2=16+1-I 987 W(I)=WW/(1.0D0+RR)**2 988 W(I2)=WW/(1.0D0-RR)**2 9891 CONTINUE 990 SUM=0.0D0 991 DO 2 I=1,16 9922 SUM=SUM+A1(I)*A2(I)*W(I) 993 SUM=6.0D0*SUM 994 RETURN 995 END 996#endif 997C /* Deck invmat */ 998#if defined (VAR_PADDY) 999C920722: Old C6 code by Patrick Fowler 1000C These two routines are used to invert a matrix in PADDY 1001 SUBROUTINE INVMAT(A,B,MATDIM,NDIM) 1002C FIND INVERSE OF MATRIX A 1003C INPUT 1004C A : MATRIX TO BE INVERTED 1005C B : SCRATCH ARRAY 1006C MATDIM : PHYSICAL DIMENSION OF MATRICES 1007C NDIM : DIMENSION OF SUBMATRIX TO BE INVERTED 1008C 1009C OUTPUT : A : INVERSE MATRIX ( ORIGINAL MATRIX THUS DESTROYED ) 1010C WARNINGS ARE ISSUED IN CASE OF CONVERGENCE PROBLEMS ) 1011C 1012#include "implicit.h" 1013#include "priunit.h" 1014 DIMENSION A(MATDIM,MATDIM),B(MATDIM,MATDIM) 1015C 1016 DETERM=0.0D0 1017 EPSIL=0.0D0 1018 ITEST=0 1019 CALL BNDINV(A,B,NDIM,DETERM,EPSIL,ITEST,MATDIM) 1020C 1021 IF( ITEST .NE. 0 ) THEN 1022 WRITE (LUPRI,'(A,I3)') ' INVERSION PROBLEM NUMBER..',ITEST 1023 END IF 1024 NTEST = 0 1025 IF ( NTEST .NE. 0 ) THEN 1026 WRITE(LUPRI,*) ' INVERTED MATRIX ' 1027 CALL WRTMAT(A,NDIM,NDIM,MATDIM,MATDIM) 1028 END IF 1029C 1030 RETURN 1031 END 1032 SUBROUTINE BNDINV(A,EL,N,DETERM,EPSIL,ITEST,NSIZE) 1033C 1034C DOUBLE PRECISION MATRIX INVERSION SUBROUTINE 1035C FROM "DLYTAP". 1036C 1037C* DOUBLE PRECISION E,F 1038C* DOUBLE PRECISION A,EL,D,DSQRT,C,S,DETERP 1039#include "implicit.h" 1040 DIMENSION A(NSIZE,1),EL(NSIZE,1) 1041 IF(N.LT.2)GO TO 140 1042 ISL2=0 1043 K000FX=2 1044 IF(ISL2.EQ.0)INDSNL=2 1045 IF(ISL2.EQ.1)INDSNL=1 1046C CALL SLITET(2,INDSNL) 1047C CALL OVERFL(K000FX) 1048C CALL DVCHK(K000FX) 1049C 1050C SET EL = IDENTITY MATRIX 1051 DO 30 I=1,N 1052 DO 10 J=1,N 1053 10 EL(I,J)=0.0D0 1054 30 EL(I,I)=1.0D0 1055C 1056C TRIANGULARIZE A, FORM EL 1057C 1058 N1=N-1 1059 M=2 1060 DO 50 J=1,N1 1061 DO 45 I=M,N 1062 IF(A(I,J).EQ.0.0D0)GO TO 45 1063 D=SQRT(A(J,J)*A(J,J)+A(I,J)*A(I,J)) 1064 C=A(J,J)/D 1065 S=A(I,J)/D 1066 38 DO 39 K=J,N 1067 D=C*A(J,K)+S*A(I,K) 1068 A(I,K)=C*A(I,K)-S*A(J,K) 1069 A(J,K)=D 1070 39 CONTINUE 1071 DO 40 K=1,N 1072 D=C*EL(J,K)+S*EL(I,K) 1073 EL(I,K)=C*EL(I,K)-S*EL(J,K) 1074 EL(J,K)=D 1075 40 CONTINUE 1076 45 CONTINUE 1077 50 M=M+1 1078C CALL OVERFL(K000FX) 1079C GO TO (140,51),K000FX 1080C 1081C CALCULATE THE DETERMINANT 1082 51 DETERP=A(1,1) 1083 DO 52 I=2,N 1084 52 DETERP=DETERP*A(I,I) 1085 DETERM=DETERP 1086C CALL OVERFL(K000FX) 1087C GO TO (140,520,520),K000FX 1088C 1089C IS MATRIX SINGULAR 1090 520 F=A(1,1) 1091 E=A(1,1) 1092 DO 58 I=2,N 1093 IF(ABS(F).LT.ABS(A(I,I)))F=A(I,I) 1094 IF(ABS(E).GT.ABS(A(I,I)))E=A(I,I) 1095 58 CONTINUE 1096 EPSILP=EPSIL 1097 IF(EPSILP.LE.0)EPSILP=1.0E-8 1098 RAT=E/F 1099 IF(ABS(RAT).LT.EPSILP)GO TO 130 1100C 1101C INVERT TRIANGULAR MATRIX 1102 J=N 1103 DO 100 J1=1,N 1104C CALL SLITE(2) 1105 I=J 1106 ISL2=1 1107 DO 90 I1=1,J 1108C CALL SLITET(2,K000FX) 1109 IF(ISL2.EQ.0)K000FX=2 1110 IF(ISL2.EQ.1)K000FX=1 1111 IF(ISL2.EQ.1)ISL2=0 1112 GO TO (70,75),K000FX 1113 70 A(I,J)=1.0D0/A(I,I) 1114 GO TO 90 1115 75 KS=I+1 1116 D=0.0D0 1117 DO 80 K=KS,J 1118 80 D=D+A(I,K)*A(K,J) 1119 A(I,J)=-D/A(I,I) 1120 90 I=I-1 1121 100 J=J-1 1122C CALL OVERFL(K000FX) 1123C GO TO (140,103,103),K000FX 1124 1125C103 CALL DVCHK(K000FX) 1126C GO TO (140,105),K000FX 1127C 1128C PREMULTIPLY EL BY INVERTED TRIANGULAR MATRIX 1129 105 M=1 1130 DO 120 I=1,N 1131 DO 118 J=1,N 1132 D=0.0D0 1133 DO 107 K=M,N 1134 107 D=D+A(I,K)*EL(K,J) 1135 EL(I,J)=D 1136 118 CONTINUE 1137 120 M=M+1 1138C CALL OVERFL(K000FX) 1139C GO TO (140,123,123),K000FX 1140C 1141C RECOPY EL TO A 1142 123 DO 124 I=1,N 1143 DO 124 J=1,N 1144 124 A(I,J)=EL(I,J) 1145 ITEST=0 1146C126 IF(INDSNL.EQ.1)CALL SLITE(2) 1147 126 IF(INDSNL.EQ.1)ISL2=1 1148 RETURN 1149C 1150 130 ITEST=1 1151 GO TO 126 1152 140 ITEST=-1 1153 GO TO 126 1154 END 1155#endif 1156! -- end of rsp/rspc8.F -- 1157