PROGRAM SDR IMPLICIT REAL*8 (A-H,O-Z) C C Calculate the spectral densities ratio J1(w)/J2(2w) using the C Full Biaxial diffusional model, as a function of the parameter C epsilon=(Dx-Dy)/(Dx+Dy) C C........please quote when using: C R. Tarroni, C. Zannoni,J. Chem. Phys. 95, 4550-4564 (1991) C input parameters: C P2 order parameter of the phase C CLAMBDA biaxiality of the molecule (e.g. solute) C OMEGA frequency of the spectral densities (in GHz) C BETA tilt angle of the quadrupole principal frame with C respect the Z axis of the molecule C RHO parameter of the Full Biaxial model corresponding to C (Dx+Dy)/2 C ETA parameter of the Full Biaxial model corresponding to C 2Dz/(Dx+Dy) C EPSSTEP step for the asymmetry parameter of the diffusional C tensor epsilon=(Dx-Dy)/(Dx+Dy) C EPSSTART starting value for epsilon C EPSEND last value for epsilon C READ (5,200) P2 READ (5,200) CLAMBDA READ (5,200) OMEGA READ (5,200) BETA READ (5,200) RHO READ (5,200) ETA READ (5,200) EPSSTEP READ (5,200) EPSSTART READ (5,200) EPSEND C CALL D2BAR (P2,D202,CLAMBDA,A20,A22,2) CALL D4BAR (D400,D402,D404,A20,A22) C WRITE (6,100) A20,A22,CLAMBDA WRITE (6,101) P2,D202,D400,D402,D404 WRITE (6,102) BETA,OMEGA WRITE (6,103) RHO,ETA WRITE (6,*) C OMEGA2= 2.*OMEGA A40 = 0. NSTEP = (EPSEND-EPSSTART)/EPSSTEP + 1 WRITE (6,105) WRITE (6,104) WRITE (6,105) C DO 50 I = 1,NSTEP EPS = EPSSTART + (I-1) * EPSSTEP DX = RHO * (1.D0 + EPS) DY = RHO * (1.D0 - EPS) DZ = RHO * ETA C CALL JM(1,OMEGA,BETA,A20,A22,A40,DX,DY,DZ,SDENS1) CALL JM(2,OMEGA2,BETA,A20,A22,A40,DX,DY,DZ,SDENS2) C R = SDENS1/SDENS2 WRITE (6,110) EPS,SDENS1,SDENS2,R 50 CONTINUE WRITE (6,105) C 100 FORMAT (' a20 =',F10.5,' a22 =',F10.5,' lambda =',F10.5) 101 FORMAT (' =',F10.5,' =',F10.5,' =',F10.5, > ' =',F10.5,' =',F10.5) 102 FORMAT (' BETA =',F10.5,' OMEGA =',F10.5) 103 FORMAT (' RHO =',F10.5,' ETA =',F10.5) 104 FORMAT (' EPSILON J1 J2 J1/J2') 105 FORMAT (' ------------------------------------------------') 110 FORMAT (8F12.6) 200 FORMAT (F10.0) STOP END C ----------------------------------------------------------------------------- subroutine JM (M,OMEGA,BETA,R20,R22,R40,DX,DY,DZ,SDENS) implicit REAL*8 (A-H,O-Z) DATA PI /3.141592653589793/ RAD6 = DSQRT(6.D0) LMAX = 10 NMAX = 4 C C M rank of the spectral density C OMEGA frequency (in GHz) C BETA tilt angle of the quadrupole principal frame C R20,R22,R40 coefficients of the anisotropic potential C DX,DY,DZ diffusional coefficients (in ns-1) C SDENS spectral density C BETARAD = BETA * PI / 180.D0 C 2 C calculate the small Wigner matrices d (BETA) C 0n BETAHALF = BETARAD / 2.D0 C = dcos(BETAHALF) S = dsin(BETAHALF) D00 = 1.D0 + 6.D0 * C * C * (C * C - 1.D0) D0M1 = RAD6 * C * S * (1.D0 - 2.D0 * C * C) D01 = - D0M1 D0M2 = RAD6 * C * C * (1.D0 - C * C) D02 = D0M2 C C calculate the JM00 C call JDENS(LMAX,NMAX,M,0,0,R20,R40,R22,DX,DY,DZ,OMEGA,1,SJ00) C C calculate the JM2-2 C call JDENS(LMAX,NMAX,M,2,-2,R20,R40,R22,DX,DY,DZ,OMEGA,1,SJ2M2) C C calculate the JM20 C call JDENS(LMAX,NMAX,M,2,0,R20,R40,R22,DX,DY,DZ,OMEGA,1,SJ20) C C calculate the JM22 C call JDENS(LMAX,NMAX,M,2,2,R20,R40,R22,DX,DY,DZ,OMEGA,1,SJ22) C C calculate the JM1-1 C call JDENS(LMAX,NMAX,M,1,-1,R20,R40,R22,DX,DY,DZ,OMEGA,1,SJ1M1) C C calculate the JM11 C call JDENS(LMAX,NMAX,M,1,1,R20,R40,R22,DX,DY,DZ,OMEGA,1,SJ11) C SDENS = D00 * D00 * SJ00 + 2.D0 * D01 * D0M1 * SJ1M1 + > 2.D0 * D01 * D01 * SJ11 + 2.D0 * D02 * D0M2 * SJ2M2 + > 4.D0 * D00 * D02 * SJ20 + 2.D0 * D02 * D02 * SJ22 return end C ----------------------------------------------------------------------------- subroutine JDENS(LMAX,NMAX,M,N0,N1,R20,R40,R22,DX,DY,DZ,OMEGA, > NFREQ,S) C C Written by: R.TARRONI C C Last revised: 11.10.89 C C Calculate the spectral densities Jmnn'(w) using the Full Biaxial C Diffusional model C R. Tarroni, C. Zannoni,J. Chem. Phys. 95, 4550-4564 (1991) C Input: C LMAX,NMAX Maximum value for L and N in the diffusional C matrix C M,N0,N1 Labels of the spectral density C R20,R40,R22 Parameters of the anisotropic potential C DX,DY,DZ Diffusional coefficients C OMEGA Array containing the frequencies (in GHz) C NFREQ Number of frequencies in OMEGA C C Output: C S Array containing the spectral densities C evaluated for the various frequencies C implicit REAL*8 (A-H,O-Z) dimension COEF(400),EIGVAL(400),OMEGA(1),S(1) C W22 = R22 / 2.D0 RHO = (DX + DY) / 2.D0 EPS = (DX - DY) / (DX + DY) ETA = 2.D0 * DZ / (DX + DY) C call BXCOR (LMAX,NMAX,2,M,N0,N1,R20,R40,W22,ETA,EPS,COEF, > EIGVAL,MXD) K = 0 IPRT = iabs(mod(N0,2)) if (M .eq. 0 .and. IPRT .eq. 0) goto 6 goto 5 C Pick up the vanishing eigenvalue! (if m=0 and parity=even) 6 SMALL = 99999. K = 99999 do 2 I = 1,MXD if (dabs(EIGVAL(I)) .lt. dabs(SMALL)) then K = I SMALL = EIGVAL(I) endif 2 continue if (dabs(SMALL) .gt. 1.D-3) then write (6,100) SMALL stop 100 format (' ZERO EIGENVALUE NOT FOUND! SMALL=',F12.8,/, > ' INCREASE THE DIMENSION OF THE DIFFUSIONAL MATRIX') endif C 5 do 20 I = 1,NFREQ S(I) = 0.D0 do 10 J = 1,MXD C skip zero eigenvalue if (J .eq. K) goto 10 C skip vanishing terms if (dabs(COEF(J)) .lt. 1.D-12) goto 10 ALPHA = RHO * EIGVAL(J) S(I) = S(I) + (COEF(J) * ALPHA)/(ALPHA*ALPHA + OMEGA(I)*OMEGA(I)) 10 continue 20 continue return end C ------------------------------------------------------------------------------ SUBROUTINE D2BAR (D200,D202,LAMBDA,A20,A22,IFLAG) C C R.T. 11.10.89 C C Calculate the autoconsistent values of D200,D202,LAMBDA,A20,A22 C according to IFLAG C IFLAG = 1 : calculate D200,D202 and Lambda given A20,A22 C IFLAG = 2 : calculate A20,A22 and D202 given D200,LAMBDA C IFLAG = 3 : calculate A20,A22 and LAMBDA given D200,D202 C IMPLICIT REAL*8 (A-H,O-Z) REAL*8 LAMBDA DIMENSION CF(2),W(2),DS(2) EXTERNAL FX1,FX2,FUNZ,FUN200,FUN202 COMMON /VAR/ R38,A,B DATA NN,NPSP,W(1),W(2) /2,2,1.D0,1.D0/ DATA V1,V2,ERR /0.D0,1.D0,1.D-6/ R38 = DSQRT(3.D0/8.D0) C GOTO (100,200,300) IFLAG 100 A = -A20 B = -A22 LAMBDA = A22 /( 2. * A20) ZAB = GAUSS(FUNZ,V1,V2,ERR) R200 = GAUSS(FUN200,V1,V2,ERR) R202 = GAUSS(FUN202,V1,V2,ERR) D200 = R200/ZAB D202 = R202/ZAB RETURN 200 CF(1) = 1. CF(2) = 1. DS(1) = D200 DS(2) = LAMBDA C CALL SIMPLEX (FX1,CF,CHI,NN,1.D-12,ITC,1000,DS,W,NPSP) A20 = CF(1) A22 = CF(2) A = -CF(1) B = -CF(2) ZAB = GAUSS(FUNZ,V1,V2,ERR) R202 = GAUSS(FUN202,V1,V2,ERR) D202 = R202/ZAB RETURN 300 CF(1) = 1. CF(2) = 1. DS(1) = D200 DS(2) = D202 c CALL SIMPLEX (FX2,CF,CHI,NN,1.D-12,ITC,1000,DS,W,NPSP) A20 = CF(1) A22 = CF(2) LAMBDA = A22 / (2. * A20) RETURN END C ------------------------------------------------------------------------------ SUBROUTINE D4BAR (D400,D402,D404,A20,A22) C C R.T. 11.10.89 C C Calculate D400,D402,D404 given the potential coefficients A20,A22 C IMPLICIT REAL*8 (A-H,O-Z) COMMON /VAR/ R38,A,B DATA V1,V2,ERR /0.D0,1.D0,1.D-6/ EXTERNAL FUNZ,FUN400,FUN402,FUN404 A = - A20 B = - A22 ZAB = GAUSS(FUNZ,V1,V2,ERR) R40 = GAUSS(FUN400,V1,V2,ERR) R42 = GAUSS(FUN402,V1,V2,ERR) R44 = GAUSS(FUN404,V1,V2,ERR) D400 = R40 / ZAB D402 = R42 / ZAB D404 = R44 / ZAB RETURN END C ----------------------------------------------------------------------------- DOUBLE PRECISION FUNCTION FX1(I,CF,NN) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION CF(2) COMMON /VAR/ R38,A,B EXTERNAL FUNZ,FUN200,FUN202 DATA V1,V2,ERR /0.D0,1.D0,1.D-6/ GOTO (100,200),I 100 A = -CF(1) B = -CF(2) ZAB = GAUSS(FUNZ,V1,V2,ERR) R200 = GAUSS(FUN200,V1,V2,ERR) FX1 = R200 / ZAB RETURN 200 FX1 = CF(2)/(2.*CF(1)) RETURN END C ------------------------------------------------------------------------------ DOUBLE PRECISION FUNCTION FX2(I,CF,NN) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION CF(2) COMMON /VAR/ R38,A,B EXTERNAL FUNZ,FUN200,FUN202 DATA V1,V2,ERR /0.D0,1.D0,1.D-6/ GOTO (100,200),I 100 A = -CF(1) B = -CF(2) ZAB = GAUSS(FUNZ,V1,V2,ERR) R200 = GAUSS(FUN200,V1,V2,ERR) FX2 = R200 / ZAB RETURN 200 A = -CF(1) B = -CF(2) ZAB = GAUSS(FUNZ,V1,V2,ERR) R202 = GAUSS(FUN202,V1,V2,ERR) FX2 = R202 / ZAB RETURN END C ------------------------------------------------------------------------------ FUNCTION FUNZ(X) IMPLICIT REAL*8 (A-H,O-Z) REAL*8 MMBSI0 COMMON /VAR/ R38,A,B IOPT=1 XX=X*X ARG=B*R38*(1.-XX) FI0 = MMBSI0(IOPT, ARG,IER) FUNZ=FI0*DEXP(A*(1.5D0*XX-.5D0)) RETURN END C ------------------------------------------------------------------------------ FUNCTION FUN200(X) IMPLICIT REAL*8 (A-H,O-Z) REAL*8 MMBSI0 COMMON /VAR/ R38,A,B IOPT=1 XX=X*X ARG=R38*B*(1.D0-XX) FI0=MMBSI0(IOPT,ARG,IER) P2=1.5D0*XX-.5D0 FUN200=P2*DEXP(A*P2)*FI0 RETURN END C ------------------------------------------------------------------------------ FUNCTION FUN202(X) IMPLICIT REAL*8 (A-H,O-Z) REAL*8 MMBSI1 COMMON /VAR/ R38,A,B IOPT=1 SS=1.D0-X*X ARG=B*R38*SS FI1= MMBSI1(IOPT,ARG,IER) FUN202=SS*R38*DEXP(A*(1.5D0*X*X-.5D0))*FI1 RETURN END C ----------------------------------------------------------------------------- FUNCTION FUN400(X) IMPLICIT REAL*8 (A-H,O-Z) REAL*8 MMBSI0 COMMON /VAR/ R38,A,B IOPT=1 XX=X*X ARG=R38*B*(1.-XX) FI0= MMBSI0(IOPT,ARG,IER) P2=1.5D0*XX-.5D0 P4=(35.D0*X**4-30.D0*X**2+3.D0)/8.D0 FUN400=P4*DEXP(A*P2)*FI0 RETURN END C ----------------------------------------------------------------------------- FUNCTION FUN402(X) IMPLICIT REAL*8 (A-H,O-Z) REAL*8 MMBSI1 COMMON /VAR/ R38,A,B PARAMETER (R12=0.5D0) RAD10 = DSQRT(10.D0) IOPT=1 SS=1.D0-X*X ARG=B*R38*SS CBM2 = R12 + R12*X SBM2 = R12 - R12*X W402 = RAD10*CBM2*SBM2*(3.D0 - 14.D0*CBM2*SBM2) FI1= MMBSI1(IOPT,ARG,IER) FUN402= W402*DEXP(A*(1.5D0*X*X-.5D0))*FI1 RETURN END C ------------------------------------------------------------------------------ FUNCTION FUN404(X) IMPLICIT REAL*8 (A-H,O-Z) REAL*8 MMBSI0,MMBSI1 COMMON /VAR/ R38,A,B PARAMETER (R12=0.5D0) RAD70 = DSQRT(70.D0) IOTP = 1 XX = X*X ARG=B*R38*(1.D0-XX) IF (ARG .EQ. 0.D0) THEN FUN404 = 0.D0 RETURN ENDIF CBM2 = R12 + R12*X SBM2 = R12 - R12*X W404 = RAD70*CBM2*CBM2*SBM2*SBM2 FI0 = MMBSI0(IOTP,ARG,IER) FI1 = MMBSI1(IOTP,ARG,IER) FI2 = FI0 - 2.D0 * FI1 / ARG FUN404= W404*DEXP(A*(1.5D0*XX-.5D0))*FI2 RETURN END C ------------------------------------------------------------------------------ DOUBLE PRECISION FUNCTION GAUSS(F,A,B,EPS) C......DOUBLE PRECIS. ADAPT OF CERN D103 FOR VAX11-780 C......C. Z. 8.8.84 IMPLICIT REAL*8 (A-H,O-Z) DIMENSION W(12),X(12) PARAMETER (CONST=1.0D-12, R12=0.5D0) DATA W / 1 0.10122 85362 90376D0 , 0.22238 10344 53374D0 2 , 0.31370 66458 77887D0 , 0.36268 37833 78362D0 3 , 0.02715 24594 11754D0 , 0.06225 35239 38648D0 4 , 0.09515 85116 82493D0 , 0.12462 89712 55534D0 5 , 0.14959 59888 16577D0 , 0.16915 65193 95003D0 6 , 0.18260 34150 44924D0 , 0.18945 06104 55069D0/ DATA X / 1 0.96028 98564 97536D0 , 0.79666 64774 13627D0 2 , 0.52553 24099 16329D0 , 0.18343 46424 95650D0 3 , 0.98940 09349 91650D0 , 0.94457 50230 73233D0 4 , 0.86563 12023 87832D0 , 0.75540 44083 55003D0 5 , 0.61787 62444 02644D0, 0.45801 67776 57227D0 6 , 0.28160 35507 79259D0 , 0.09501 25098 37637D0/ DELTA= CONST* DABS(A-B) GAUSS= 0.D0 AA=A 5 Y= B-AA IF ( DABS(Y) .LE. DELTA) RETURN 2 BB=AA+Y C1= R12 *(AA+BB) C2=C1-AA S8=0.D0 S16=0.D0 DO 1 I=1,4 U= X(I)*C2 1 S8= S8+ W(I)* ( F(C1+U) + F(C1-U) ) DO 3 I= 5,12 U= X(I)*C2 3 S16=S16 + W(I) *( F(C1+U) + F(C1-U) ) S8= S8*C2 S16= S16*C2 IF ( DABS(S16-S8) .GT. EPS*(1.0D0 +DABS(S16)) ) GOTO 4 GAUSS = GAUSS+ S16 AA=BB GOTO 5 4 Y= R12 * Y IF (DABS(Y) .GT. DELTA ) GOTO 2 WRITE(6,7000) GAUSS=0.D0 RETURN 7000 FORMAT( 1X, '...ERR IN GAUSS: TOO HIGH ACCURACY REQUIRED') END C ------------------------------------------------------------------------------ DOUBLE PRECISION FUNCTION MMBSI0 (IOPT,ARG,IER) INTEGER IOPT,IER DOUBLE PRECISION ARG DOUBLE PRECISION A,CONST,DEXP40,P(15),PP(8),XSMALL,XX,CON, 1 Q(5),QQ(7),SUMP,SUMQ,X,XINF,XMAX C DATA XMAX /91.200D0/ DATA DEXP40/ .235385266837019985D+18/ DATA XSMALL/8.67361738D-19/ DATA XINF /1.7D+38/ DATA CONST / .666666666666666667D-01/ C COEFFICIENTS FOR XSMALL .LE. C ABS(ARG) .LT. 15.0 C DATA P /-.524878303872486911D-17,-.159822176963152107D-13, * -.268434353985728351D-10,-.305172129516732270D-07, * -.251726344203303272D-04,-.154539719077661595D-01, * -.709353218849083214D+01,-.241251875635161846D+04, * -.595456062371357883D+06,-.103130633806846395D+09, * -.119127423465868193D+11,-.849250748990540193D+12, * -.329400775241734956D+14,-.550503528993272799D+15, * -.223355758483236763D+16/ DATA Q /-.372775563866849905D+04, .651584945870435991D+07, * -.656265446265324359D+10, .376041775067259003D+13, * -.970879130649400689D+15/ C COEFFICIENTS FOR 15.0 .LE. ABS(ARG) C DATA PP/-.398437500000000000D+00, .292053801142487550D+01, * -.247084602559363172D+01, .479148550995084342D+00, * -.373847878179342018D-02,-.268014935222659065D-02, * .991686607024382073D-04,-.218771393882892240D-05/ DATA QQ/-.314466876245542596D+02, .855395457249381542D+02, * -.602279794663117256D+02, .139825866959438771D+02, * -.111517488948088167D+01, .325476719383372820D-01, * -.551943584858776969D-03/ C FIRST EXECUTABLE STATEMENT IER = 0 XCON = 85.0D0 IF (IOPT.NE.1.AND.IOPT.NE.2) GO TO 30 X = DABS(ARG) IF (X.LT.XSMALL) GO TO 20 IF ((IOPT.EQ.1).AND.(X.GT.XMAX)) GO TO 25 IF (X.GE.15.0D0) GO TO 10 C XSMALL .LE. ABS(ARG) .LT. 15.0 XX = X*X SUMP = P(1) C DO 5 I=2,15 SUMP = SUMP*XX+P(I) 5 CONTINUE C XX = XX-225.D0 SUMQ = ((((XX+Q(1))*XX+Q(2))*XX+Q(3))*XX+Q(4))*XX+Q(5) MMBSI0 = SUMP/SUMQ IF (IOPT.EQ.2) MMBSI0 = MMBSI0*DEXP(-X) GO TO 9005 C 15.0 .LE. ABS(ARG) 10 XX = 1.0D0/X-CONST SUMP = ((((((PP(1)*XX+PP(2))*XX+PP(3))*XX+PP(4))*XX+PP(5))*XX+ *PP(6))*XX+PP(7))*XX+PP(8) SUMQ = ((((((XX+QQ(1))*XX+QQ(2))*XX+QQ(3))*XX+QQ(4))*XX+QQ(5))*XX+ *QQ(6))*XX+QQ(7) MMBSI0 = (SUMP/SUMQ-PP(1))/DSQRT(X) IF (IOPT.EQ.2) GO TO 9005 IF (X.GT.XCON) GO TO 15 MMBSI0 = MMBSI0*DEXP(X) GO TO 9005 15 A = MMBSI0*DEXP(X-40.0D0) MMBSI0 = A*DEXP40 GO TO 9005 20 MMBSI0 = 1.0D0 GO TO 9005 25 MMBSI0 = XINF IER = 130 GO TO 9000 30 MMBSI0 = XINF IER = 129 9000 CONTINUE STOP 99999 9005 RETURN END C ------------------------------------------------------------------------------ DOUBLE PRECISION FUNCTION MMBSI1 (IOPT,ARG,IER) INTEGER IOPT,IER DOUBLE PRECISION ARG DOUBLE PRECISION A,CONST,DEXP40,P(15),PBAR,PP(8),XCON, 1 Q(5),QQ(6),SUMP,SUMQ,X,XINF,XMAX,XMAX1,XMINS, 2 XSMALL, XX INTEGER J C DATA XCON/85.0D0/ DATA XMAX1 /.18014398509481984D17/ DATA XMAX /91.200D0/ DATA DEXP40/ .235385266837019985D+18/ DATA XSMALL/8.67361738D-19/ DATA XINF /1.7D+38/ DATA CONST / .666666666666666667D-01/ DATA XMINS/1.058791184D-22/ C COEFFICIENTS FOR XSMALL .LE. C ABS(ARG) .LT. 15.0 C DATA P /-.197052918025351399D-18,-.652455155831519027D-15, * -.119287889036032387D-11,-.148319049359946477D-08, * -.134668298276351529D-05,-.917464432878175013D-03, * -.472070908273101624D+00,-.182259466316573159D+03, * -.518940919823080174D+05,-.105885507247693471D+08, * -.148282676066123661D+10,-.133574376822754930D+12, * -.698767796480100898D+13,-.177320378407915915D+15, * -.145771802781434636D+16/ DATA Q /-.400768646799041899D+04, .748105803566550690D+07, * -.800595189986197647D+10, .485447142582736225D+13, * -.132181683073214422D+16/ C COEFFICIENTS FOR 15.0 .LE. ABS(ARG) C DATA PP/-.604371590561376000D-01, .457481229019334590D+00, * -.428437669033048064D+00, .973560001508866121D-01, * -.324577239744655683D-02,-.363952647121217953D-03, * .162586618674408364D-04,-.363475784046082235D-06/ DATA QQ/-.388065867215565934D+01, .325937148890369963D+01, * -.850174764632179244D+00, .742120108131865300D-01, * -.228356244894925126D-02, .375104331119228246D-04/ DATA PBAR / .398437500000000000D+00/ IER = 0 IF (IOPT.NE.1.AND.IOPT.NE.2) GO TO 30 X = DABS(ARG) IF (X.LT.XSMALL) GO TO 20 IF ((IOPT.EQ.1).AND.(X.GT.XMAX)) GO TO 25 IF (X.GE.15.0D0) GO TO 10 C XSMALL .LE. ABS(ARG) .LT. 15.0 XX = X*X SUMP = P(1) C DO 5 J=2,15 SUMP = SUMP*XX+P(J) 5 CONTINUE C XX = XX-225.0D0 SUMQ = ((((XX+Q(1))*XX+Q(2))*XX+Q(3))*XX+Q(4))*XX+Q(5) MMBSI1 = (SUMP/SUMQ)*X IF (IOPT.EQ.2) MMBSI1 = MMBSI1*DEXP(-X) GO TO 35 C 15.0 .LE. ABS(ARG) 10 XX = 1.0D0/X-CONST SUMP = ((((((PP(1)*XX+PP(2))*XX+PP(3))*XX+PP(4))*XX+PP(5))*XX+ *PP(6))*XX+PP(7))*XX+PP(8) SUMQ = (((((XX+QQ(1))*XX+QQ(2))*XX+QQ(3))*XX+QQ(4))*XX+QQ(5))*XX+ *QQ(6) MMBSI1 = (SUMP/SUMQ+PBAR)/DSQRT(X) IF (IOPT.EQ.2) GO TO 35 IF (X.GT.XCON) GO TO 15 MMBSI1 = MMBSI1*DEXP(X) GO TO 35 C CALCULATION REFORMULATED TO AVOID C PREMATURE OVERFLOW AND TO FOIL THE C COMPILER OPTIMIZATION 15 A = MMBSI1*DEXP(X-40.0D0) MMBSI1 = A*DEXP40 GO TO 35 C RETURN FOR ABS(ARG) .LT. XSMALL 20 IF (X*XMAX1.LT.XMINS) X = 0.0D0 MMBSI1 = 0.5D0*X GO TO 35 C ERROR RETURN FOR ABS(ARG) .GT. XMAX 25 MMBSI1 = XINF IER = 130 GO TO 35 30 IER = 129 MMBSI1 = XINF 35 IF (ARG.LT.0.0D0) MMBSI1 = -MMBSI1 9000 CONTINUE 9005 RETURN END C ----------------------------------------------------------------------------- SUBROUTINE BXCOR(LMAX,NMAX,LR,MR,NR1,NR2,A20,A40,A22,ETA, > EPS,COEF,EIGVAL,MXD) C C Written by: R.TARRONI, Dipartimento di Chimica Fisica ed Inorganica, C Universita' di Bologna C C Last updated: 26.7.89 C C LR,LR C Calculation of the orient. correlations G (t) ; LR = 1,2 C MR,NR1,MR,NR2 C using the C C 'F U L L B I A X I A L' Diffusional Model C C Anisotropic potential of the type: C C 2 2 2 4 C U(r)/kT = a20 D (r) + a22 (D (r) + D (r)) + a40 D (r) C 00 02 0-2 00 C C Versione in doppia precisione per VAX 11/780 C N. massimo di funzioni di base: 400 (max dimensione matrici) C C C LR,MR,NR1,NR2 : Indici della funzione di correlazione da calcolare C LMAX,NMAX : Indici L,n massimi della matrice diffusionale. Il C programma e' dimensionato per accettare LMAX=40, C NMAX=10 come valori massimi, ai quali corrisponde una C matrice diffusionale di dimensioni 400x400 C A20,A40,A22, : Parametri del potenziale anisotropo C EPS,ETA Parametri del tensore di diffusione. C MXD : Dimensione della matrice diffusionale corrispondente a C LMAX,NMAX C EIGVAL : Vettore contenente gli autovalori della matrice C diffusionale, cambiati di segno e ordinati in base ai C corrispondenti fattori di intensita' C COEF : Fattori di intensita' relativi agli autovalori C C N.B. Il programma assume rho=1, per cui gli autovalori in uscita devono C essere scalati per un rho diverso C IMPLICIT REAL*8(A-H,O-Z) DIMENSION COEF(1),EIGVAL(1) COMMON /DIFF/ QMAT(400,100),V(400,3),B0(400) DIMENSION DG(400),SBDG(400),AUX1(400),AUX2(400),AUX3(400) EQUIVALENCE (QMAT(1,1),DG(1)),(QMAT(1,2),SBDG(1)) COMMON /DIM/ MXDIM,LBDIM,MXPAR,LXDIM,NXDIM COMMON /AUX2/ MXD2,NMAX2,KA2(400,2),KB2(0:40,-10:10) COMMON /AUX3/ MXD3,NMAX3,KA3(400,2),KB3(0:40,-10:10) COMMON /SAVE/ ZA20,ZA40,ZA22,ZETA,ZEPS,NPRT,NMR,LV(-2:2), > LZ(0:1,3),LQ(0:1) DATA MXDIM,LBDIM,MXPAR,NXDIM,LXDIM /400,100,400,10,40/ DATA LV(-2),LV(-1),LV(0),LV(1),LV(2) /1,1,2,2,3/ DATA LZ(0,1),LZ(0,2),LZ(0,3) /-2,0,2/ DATA LZ(1,1),LZ(1,2) /-1,1/ DATA LQ(0),LQ(1) /3,2/ DATA ZA20,ZA40,ZA22,ZETA,ZEPS/-9999.,-9999.,-9999.,-9999.,-9999./ DATA NMR,NPRT /9999,9999/ C IERR = 0 IF (LR .LT. 1 .OR. LR .GT. 2) IERR = 5 IF (IABS(MR) .GT. LR) IERR = 6 IF (IABS(NR1) .GT. LR) IERR = 6 IF (IABS(NR2) .GT. LR) IERR = 6 C ...controlla la parita' IF (MOD(NR1+NR2,2) .EQ. 1) IERR = 7 IF (LMAX .GT. LXDIM .OR. NMAX .GT. NXDIM) IERR = 8 IF (IERR .NE. 0) GOTO 2000 C C controlla se e' richiesto nel potenziale il contributo di tipo P4 C IP4 = 0 IF (DABS(A40) .GT. 1.D-12) IP4 = 1 IPRT = IABS(MOD(NR1,2)) C 1/2 C calcolo dei coefficienti della distribuzione all'equilibrio P (beta,gamma) C per il potenziale richiesto C IF (DABS(A20-ZA20) .GT. 1.D-12) GOTO 4 IF (DABS(A40-ZA40) .GT. 1.D-12) GOTO 4 IF (DABS(A22-ZA22) .GT. 1.D-12) GOTO 4 GOTO 12 4 CALL AUX(IP4,LMAX,NMAX,0,0,KA3,KB3,MXD3,LBANDA3,NMAX3,IERR) IF (IERR .NE. 0) GOTO 2000 CALL STOC(A20,A40,A22,0.,0.,0,KA3,KB3,MXD3,LMAX,NMAX3) NSTEP = MXD3 / 2 CALL PEQUIL(QMAT,MXDIM,LBDIM,LBANDA3,MXD3,NSTEP,JSTEP,B0,AUX1, > AUX2,AUX3) ZA20 = A20 ZA40 = A40 ZA22 = A22 GOTO 14 C C ...parita' della matrice diffusionale C 12 IF (DABS(ZETA-ETA) .GT. 1.D-12) GOTO 14 IF (DABS(ZEPS-EPS) .GT. 1.D-12) GOTO 14 IF (IPRT .NE. NPRT) GOTO 14 IF (MR .NE. NMR) GOTO 14 GOTO 16 14 CALL AUX(IP4,LMAX,NMAX,MR,IPRT,KA2,KB2,MXD2,LBANDA2,NMAX2,IERR) C C riempi la matrice diffusionale C CALL STOC(A20,A40,A22,EPS,ETA,MR,KA2,KB2,MXD2,LMAX,NMAX2) ZETA = ETA ZEPS = EPS NPRT = IPRT NMR = MR C C calcola gli elementi dei vettori V C NVECT = LQ(IPRT) DO 22 K = 1,NVECT NX = LZ(IPRT,K) DO 21 I = 1,MXD2 I0 = KA2(I,1) N0 = KA2(I,2) V(I,K) = 0. N = NX - N0 IF (IABS(N) .GT. NMAX3) GOTO 21 C1 = I0 * 2. + 1. L1 = IABS(I0-LR) L2 = MIN0(I0+LR,LMAX) DO 20 L0 = L1,L2 IF (IABS(N) .GT. L0) GOTO 20 CCG = DSQRT(C1/(L0*2.+1.)) * CG(LR,I0,L0,MR,-MR) * > CG(LR,I0,L0,NX,-N0) * (-1.)**NX V(I,K) = V(I,K) + CCG * B0(KB3(L0,N)) 20 CONTINUE 21 CONTINUE 22 CONTINUE C C Per la diagonalizzazione vengono utilizzate le routine riportate sulla tesi C di dottorato di C.F.Polnaszek, specifiche per l'applicazione a matrici a C banda reali e simmetriche. Per risparmiare tempo calcolo e memoria i vettori C di V vengono ruotati all'interno delle routine stesse C RSQZ applica l'algoritmo di Rutishauser e trasforma la matrice a banda in C una matrice tridiagonale C RQRT calcola gli autovalori della matrice tridiagonale applicando C l'algoritmo QR C CALL RSQZ(QMAT,V,MXDIM,LBDIM,MXDIM,NVECT,MXD2,LBANDA2,NVECT,SQTOL) C C controlla se la matrice e' gia' diagonale C SBDGMAX = 0.D0 DO 18 I = 1,MXD2 IF (DABS(SBDG(I)) .GT. DABS(SBDGMAX)) SBDGMAX=SBDG(I) 18 CONTINUE IF (DABS(SBDGMAX) .LT. 1.D-12) GOTO 16 EIGVAL (1) = 0.D0 MM = MXDIM * NVECT CALL RQRT(DG,SBDG,V,MXDIM,MXD2,MXDIM,NVECT,MM,TOL,EIGVAL,IERR) IF (IERR .NE. 0) GOTO 2000 C C calcola la funzioni di correlazione C 16 IK = LV(NR1) JK = LV(NR2) DO 24 I = 1,MXD2 COEF(I) = V(I,IK) * V(I,JK) EIGVAL(I) = -DG(I) 24 CONTINUE C C sistema i fattori di peso in ordine crescente C DO 26 I = 1,MXD2 DO 26 J = I,MXD2 IF (DABS(COEF(I)) .GT. DABS(COEF(J))) GOTO 26 SAVE = COEF(I) COEF(I) = COEF(J) COEF(J) = SAVE SAVE = EIGVAL(I) EIGVAL(I) = EIGVAL(J) EIGVAL(J) = SAVE 26 CONTINUE MXD = MXD2 GOTO 62 C 2000 GOTO (171,172,173,174,175,176,177,178),IERR 171 WRITE (6,71) GOTO 62 172 WRITE (6,72) GOTO 62 173 WRITE (6,73) GOTO 62 174 WRITE (6,74) GOTO 62 175 WRITE (6,75) GOTO 62 176 WRITE (6,76) GOTO 62 177 WRITE (6,77) GOTO 62 178 WRITE (6,78) 62 RETURN C 71 FORMAT (' ***** ERRORE! RQRT non ha raggiunto convergenza nel C numero di iterazioni consentito') 72 FORMAT (' ***** ERRORE! Numero dei termini dell''espansione C superiore alla dimensione della matrice diffusionale') 73 FORMAT (' ***** ERRORE! Larghezza della banda della matrice C diffusionale superiore ai limiti fissati per QMAT') 74 FORMAT (' ***** ERRORE! Dimensione della matrice C diffusionale superiore ai limiti fissati per QMAT') 75 FORMAT (' ***** ERRORE! solo funzioni di rango 1 e 2') 76 FORMAT (' ***** ERRORE! LR1,LR2,MR,NR1,NR2 incompatibili') 77 FORMAT (' ***** ERRORE! NR1 + NR2 deve essere pari') 78 FORMAT (' ***** ERRORE! LMAX oppure NMAX superiori i limiti C fissati') END C C ----------------------------------------------------------------------------- C SUBROUTINE AUX (IP4,LMAX,NMAX,MR,IPRT,KA,KB,MXD,LB,NP,IERR) C C Riempimento delle matrici KA e KB, necessarie per effettuare le conversioni C tra indici della matrice diffusionale e numeri quantici L,N. C C --( KA )--> C indice I numeri quantici L,N C <--( KB )-- C C Vengono calcolate inoltre la dimensione della matrice diffusionale C e la larghezza della banda; qualora tali dimensioni siano superiori C a quelle consentite viene segnalato un errore. C C MXD : dimensione della matrice corrspondente a LMAX,NMAX C LB : larghezza della banda C NP : N massimo consentito da NMAX e parita' C IERR : codice errore C COMMON /DIM/ MXDIM,LBDIM,MXPAR,LXDIM,NXDIM DIMENSION KA(MXDIM,2),KB(0:LXDIM,-NXDIM:NXDIM) DO 2 I = 0,LXDIM DO 2 J = -NXDIM,NXDIM KB(I,J) = 0 2 CONTINUE DO 4 I = 1,MXDIM KA(I,1) = 0 KA(I,2) = 0 4 CONTINUE C C calcola la larghezza della banda C IF (MOD(NMAX,2) .EQ. 0) THEN NP = NMAX - IPRT ELSE NP = NMAX + IPRT - 1 ENDIF IF (IP4 .EQ. 0) LB = 5 * (NP + 1) IF (IP4 .EQ. 1) LB = 9 * (NP + 1) IF (LB .GT. LBDIM) THEN IERR = 3 RETURN ENDIF LI = MAX0(MR,IPRT) K = 0 DO 10 L = LI,LMAX IF (IPRT .EQ. 0) THEN NNX = L - MOD(L,2) ELSE NNX = L + MOD(L,2) - IPRT ENDIF NX = MIN0(NNX,NP) DO 10 N = -NX,NX,2 K = K + 1 IF (K .GT. MXDIM) THEN IERR = 4 RETURN ENDIF KA(K,1) = L KA(K,2) = N KB(L,N) = K 10 CONTINUE MXD = K RETURN END C C ----------------------------------------------------------------------------- C SUBROUTINE STOC(A20,A40,A22,EPS,ETA,MR,KA,KB,MXD,LMAX,NMAX) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION KA(MXDIM,2),KB(0:LXDIM,-NXDIM:NXDIM) COMMON /DIM/ MXDIM,LBDIM,MXPAR,LXDIM,NXDIM COMMON /DIFF/ QMAT(400,100),V(400,3),B0(400) DIMENSION S(37),CST(10) C C Riempie la matrice diffusionale C N.B. La matrice diffusionale e' simmetrica e a banda, per cui vengono C calcolati solo gli elementi relativi alla semibanda superiore C C ... CG(2,2,0,0,0) S(1) = DSQRT(1.D0/5.D0) C ... CG(2,2,0,1,-1) S(2) = - DSQRT(1.D0/5.D0) C ... CG(2,2,0,2,-2) S(3) = DSQRT(1.D0/5.D0) C ... CG(2,2,2,0,0) S(4) = - DSQRT(2.D0/7.D0) C ... CG(2,2,2,1,-1) S(5) = DSQRT(1.D0/14.D0) C ... CG(2,2,2,2,-2) S(6) = DSQRT(2.D0/7.D0) C ... CG(2,2,2,1,1) S(7) = - DSQRT(3.D0/7.D0) C ... CG(2,2,4,0,0) S(8) = DSQRT(18.D0/35.D0) C ... CG(2,2,4,1,-1) S(9) = DSQRT(8.D0/35.D0) C ... CG(2,2,4,2,-2) S(10) = DSQRT(1.D0/70.D0) C ... CG(2,2,4,1,1) S(11) = DSQRT(4.D0/7.D0) C ... CG(2,2,4,2,2) S(12) = 1.D0 C ... CG(2,4,0,0,0) S(13) = 0. C ... CG(2,4,0,1,-1) S(14) = 0. C ... CG(2,4,2,0,0) S(15) = DSQRT(2.D0/7.D0) C ... CG(2,4,2,1,-1) S(16) = - DSQRT(5.D0/21.D0) C ... CG(2,4,2,1,1) S(17) = - DSQRT(5.D0/126.D0) C ... CG(2,4,4,0,0) S(18) = - DSQRT(20.D0/77.D0) C ... CG(2,4,4,1,-1) S(19) = DSQRT(3.D0/154.D0) C ... CG(2,4,4,1,1) S(20) = - DSQRT(243.D0/1540.D0) C ... CG(2,4,6,0,0) S(21) = DSQRT(5.D0/11.D0) C ... CG(2,4,6,1,-1) S(22) = DSQRT(8.D0/33.D0) C ... CG(2,4,6,1,1) S(23) = DSQRT(224.D0/495.D0) C ... CG(4,4,0,0,0) S(24) = DSQRT(1.D0/9.D0) C ... CG(4,4,0,1,-1) S(25) = - DSQRT(1.D0/9.D0) C ... CG(4,4,2,0,0) S(26) = - DSQRT(100.D0/693.D0) C ... CG(4,4,2,1,-1) S(27) = DSQRT(289.D0/2772.D0) C ... CG(4,4,2,1,1) S(28) = - DSQRT(50.D0/231.D0) C ... CG(4,4,4,0,0) S(29) = DSQRT(162.D0/1001.D0) C ... CG(4,4,4,1,-1) S(30) = - DSQRT(81.D0/2002.D0) C ... CG(4,4,4,1,1) S(31) = DSQRT(180.D0/1001.D0) C ... CG(4,4,6,0,0) S(32) = - DSQRT(20.D0/99.D0) C ... CG(4,4,6,1,-1) S(33) = - DSQRT(1.D0/1980.D0) C ... CG(4,4,6,1,1) S(34) = - DSQRT(7.D0/33.D0) C ... CG(4,4,8,0,0) S(35) = DSQRT(490.D0/1287.D0) C ... CG(4,4,8,-1,1) S(36) = DSQRT(1568.D0/6435.D0) C ... CG(4,4,8,1,1) S(37) = DSQRT(56.D0/143.D0) C CST(1) = 5.D0 * A40 * A40 * S(24) * S(25) > + 1.5D0 * A20 * A20 * S(1) * S(2) > + A22 * A22 * S(1) * S(2) > - 2.D0 * ETA * A22 * A22 * S(1) * S(3) > + EPS * DSQRT(6.D0) * A20 * A22 * S(1) * S(2) CST(2) = - 3.D0 * A20 > + 1.5D0 * A20 * A20 * S(4) * S(5) > + A22 * A22 * S(4) * S(5) > - 2.D0 * ETA * A22 * A22 * S(4) * S(6) > + DSQRT(30.D0) * A20 * A40 * S(15) * S(16) > + 5.D0 * A40 * A40 * S(26) * S(27) > - EPS * DSQRT(6.D0) * A22 > + EPS * DSQRT(6.D0) * A20 * A22 * S(4) * S(5) > + EPS * DSQRT(20.D0) * A22 * A40 * S(15) * S(16) CST(3) = - (1.D0 + 2.D0 * ETA) * A22 > + DSQRT(1.5D0) * A20 * A22 * S(4) * S(7) > + DSQRT(5.D0) * A40 * A22 * S(15) * S(17) > - 0.5D0 * DSQRT(6.D0) * EPS * A20 > + 0.25D0*EPS*(3.D0*A20*A20+2.D0*A22*A22)*S(4)*S(7) > + 2.5D0 * EPS * A40 * A40 * S(26) * S(28) > + 0.25D0 * DSQRT(120.D0) * EPS * A20 * A40 * S(15)*S(17) CST(4) = - 10.D0 * A40 > + 1.5D0 * A20 * A20 * S(8) * S(9) > + A22 * A22 * S(8) * S(9) > - 2.D0 * ETA * A22 * A22 * S(8) * S(10) > + DSQRT(30.D0) * A20 * A40 * S(18) * S(19) > + 5.D0 * A40 * A40 * S(29) * S(30) > + EPS * DSQRT(6.D0) * A20 * A22 * S(8) * S(9) > + EPS * DSQRT(20.D0) * A22 * A40 * S(18) * S(19) CST(5) = DSQRT(1.5D0) * A20 * A22 * S(8) * S(11) > + DSQRT(5.D0) * A40 * A22 * S(18) * S(20) > - 1.5D0 * EPS * DSQRT(10.D0) * A40 > + 0.25D0*DSQRT(120.D0)*EPS*A20*A40*S(18)*S(20) > + 0.25D0*EPS*(3.D0*A20*A20+2.D0*A22*A22)*S(8)*S(11) > + 2.5D0 * EPS * A40 * A40 * S(29) * S(31) CST(6) = ETA * A22 * A22 * S(8) * S(12) CST(7) = 5.D0 * A40 * A40 * S(32) * S(33) > + DSQRT(30.D0) * A20 * A40 * S(21) * S(22) > + EPS * DSQRT(20.D0) * A22 * A40 * S(21) * S(22) CST(8) = 2.5D0 * EPS * A40 * A40 * S(32) * S(34) > + 0.25D0*DSQRT(120.D0)*EPS*A20*A40*S(21)*S(23) > + DSQRT(5.D0) * A40 * A22 * S(21) * S(23) CST(9) = 5.D0 * A40 * A40 * S(35) * S(36) CST(10) = 2.5D0 * EPS * A40 * A40 * S(35) * S(37) C DO 2 I = 1,MXDIM DO 2 J = 1,LBDIM QMAT(I,J) = 0.D0 2 CONTINUE DO 4 I = 1,MXD L0 = KA(I,1) N0 = KA(I,2) QMAT(I,1) = - (L0 * (L0+1.D0)+ N0 * N0 * (ETA-1.D0)) + CST(1) N = N0 + 2 IN = IABS(N) IF (IN .GT. L0 .OR. IN .GT. NMAX) GOTO 4 K = KB(L0,N) - I + 1 IF (K .LT. 1) GOTO 4 QMAT(I,K) = QMAT(I,K) - 0.5D0 * EPS * DSQRT((L0*(L0+1.D0) - > (N0+2.D0)*(N0+1.D0)) * (L0*(L0+1.D0) - (N0+1.D0)*N0)) 4 CONTINUE C CALL CALC(KA,KB,MXD,LMAX,NMAX,MR,2,0,CST(2)) CALL CALC(KA,KB,MXD,LMAX,NMAX,MR,2,2,CST(3)) CALL CALC(KA,KB,MXD,LMAX,NMAX,MR,4,0,CST(4)) CALL CALC(KA,KB,MXD,LMAX,NMAX,MR,4,2,CST(5)) CALL CALC(KA,KB,MXD,LMAX,NMAX,MR,4,4,CST(6)) IF (DABS(A40) .LT. 1E-12) RETURN CALL CALC(KA,KB,MXD,LMAX,NMAX,MR,6,0,CST(7)) CALL CALC(KA,KB,MXD,LMAX,NMAX,MR,6,2,CST(8)) CALL CALC(KA,KB,MXD,LMAX,NMAX,MR,8,0,CST(9)) CALL CALC(KA,KB,MXD,LMAX,NMAX,MR,8,2,CST(10)) RETURN END C C ----------------------------------------------------------------------------- C SUBROUTINE CALC(KA,KB,MXD,LMAX,NMAX,MR,LK,NK,CST) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION KA(MXDIM,2),KB(0:LXDIM,-NXDIM:NXDIM) COMMON /DIM/ MXDIM,LBDIM,MXPAR,LXDIM,NXDIM COMMON /DIFF/ QMAT(400,100),V(400,3),B0(400) C DO 40 JJ = 1,MXD L0 = KA(JJ,1) N0 = KA(JJ,2) L2 = MAX0(IABS(L0-LK),MR,L0) L3 = MIN0(L0+LK,LMAX) C2 = 2.*L0 + 1 DO 40 L1 = L2,L3 C3 = DSQRT((2.*L1+1.)/C2)*CG(L1,LK,L0,MR,0) N = N0 - NK IN = IABS(N) IF (IN.GT.L1.OR.IN.GT.NMAX) GOTO 45 JCK = KB(L1,N) - JJ + 1 IF (JCK.LT.1) GOTO 45 ST = C3 * CG(L1,LK,L0,N,NK) QMAT(JJ,JCK) = QMAT(JJ,JCK) + ST * CST 45 IF (NK.EQ.0) GOTO 40 N = N0 + NK IN = IABS(N) IF (IN.GT.L1.OR.IN.GT.NMAX) GOTO 40 JCK = KB(L1,N) - JJ + 1 IF (JCK.LT.1) GOTO 40 ST = C3 * CG(L1,LK,L0,N,-NK) QMAT(JJ,JCK) = QMAT(JJ,JCK) + ST * CST 40 CONTINUE RETURN END C ------------------------------------------------------------------------------ SUBROUTINE PEQUIL (A,MXDIM,LBDIM,LB,N,NSTEP,JSTEP,B,R,W,Q) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(MXDIM,LBDIM),B(1),R(1),W(1),Q(1) C C Calculate the equilibrium distribution by solving a linear equation system C on the diffusion matrix A labelled with MR=0,IPRT=even C The linear equation system is solved using a coniugate gradient procedure C (see G.H.Golub, C.F.Van Loan, "Matrix Computation", pag.372, C North Oxford Academic press, 1983) C JSTEP = 1 M = N - 1 C C define the starting vector R and the first BETA C DO 2 I = 1,N R(I) = 0.D0 W(I) = 0.D0 Q(I) = 0.D0 B(I) = 0.D0 2 CONTINUE B(1) = 1.D0 RHO1 = 0.D0 DO 4 I = 2,LB R(I-1) = A(1,I) RHO1 = RHO1 + R(I-1) * R(I-1) 4 CONTINUE C DO 6 I = 1,M 6 Q(I) = R(I) GOTO 28 18 BETA = RHO1/RHO0 DO 26 I = 1,M 26 Q(I) = R(I) + BETA * Q(I) C T C calculate the vector W = A Q and the new GAMMA = Q A Q C 28 GAMMA = 0.D0 S = 0.D0 DO 8 K = 1,LB S = S + A(2,K) * Q(K) 8 CONTINUE W(1) = S GAMMA = GAMMA + Q(1) * W(1) DO 10 I = 3,N S = 0.D0 JJ = MAX0(I-LB+1,2) DO 12 K = JJ,I-1 S = S + A(K,I-K+1) * Q(K-1) 12 CONTINUE JJ = MIN0(I+LB-1,N) DO 14 K = I,JJ S = S + A(I,K-I+1) * Q(K-1) 14 CONTINUE W(I-1) = S GAMMA = GAMMA + Q(I-1) * W(I-1) 10 CONTINUE C IF (GAMMA .LT. 1.D-12) GOTO 20 ALPHA = RHO1 / GAMMA RHO0 = RHO1 RHO1 = 0.D0 DO 16 I = 1,M B(I+1) = B(I+1) - ALPHA * Q(I) R(I) = R(I) - ALPHA * W(I) RHO1 = RHO1 + R(I) * R(I) 16 CONTINUE JSTEP = JSTEP + 1 IF (RHO1 .LT. 1.D-12) GOTO 20 IF (JSTEP .GT. M) GOTO 20 IF (JSTEP .GT. NSTEP) GOTO 20 GOTO 18 C C calculate the equilibrium distribution by normalizing the vector B C 20 FNORM = 0.D0 DO 22 I = 1,N FNORM = FNORM + B(I) * B(I) 22 CONTINUE FNORM = DSQRT(FNORM) DO 24 I = 1,N B(I) = B(I) / FNORM 24 CONTINUE RETURN END C C ----------------------------------------------------------------------------- C DOUBLE PRECISION FUNCTION CG(L1,L2,L3,M1,M2) C C Restituisce il coefficiente di Clebsch-Gordan C(L1,L2,L3;M1,M2,M1+M2) C corrispondente C IMPLICIT REAL*8 (A-H,O-Z) A = DFLOAT(L1) B = DFLOAT(L2) C = DFLOAT(L3) X = DFLOAT(M1) Y = DFLOAT(M2) Z = DFLOAT(M1+M2) IF (IDINT(C-A-B)) 301,301,302 302 CG = 0.D0 GOTO 399 301 IF (IDINT(A-B-C)) 303,303,302 303 IF (IDINT(B-A-C)) 304,304,302 304 IF (IDINT(Z-X-Y)) 302,305,302 305 IF (DABS(Z) .GT. C) GOTO 302 IF (DABS(Y) .GT. B) GOTO 302 IF (DABS(X) .GT. A) GOTO 302 IF (IDINT(B) .NE. 0) GOTO 5 IF (IDINT(A) .EQ. IDINT(C)) THEN CG = 1.D0 ELSE CG = 0.D0 ENDIF GOTO 399 5 IF (IDINT(C) .NE. 0) GOTO 6 IF (IDINT(A) .EQ. IDINT(B)) THEN FMULT = (-1.D0)**(IDINT(A-X)) CG = FMULT/DSQRT(2.D0*A + 1.D0) ELSE CG = 0.D0 ENDIF GOTO 399 6 NMAX = IDINT(DMIN1(C-A+B,C+Z)) + 1 NMIN = IDINT(DMAX1(0.D0,B-A+Z)) + 1 SUM = 0.D0 CG = 1.D0 X1 = A + X Y1 = B + Y Z1 = C + Z X2 = A - X Y2 = B - Y Z2 = C - Z Z3 = 2.D0 * C + 1.D0 R = A + B + C + 1.D0 S = C + B + X T = B + C - A J = IDINT (Y1) DO 306 I = NMIN,NMAX N = I - 1 U = DFLOAT(N) 306 SUM = SUM + (-1.D0)**N * BICO(Z1,U) * BICO(S-U,X1) * C BICO(X2+U,Y1) SUM = SUM * (-1.D0)**J IF (IDINT(X)) 11,12,13 11 G = 1.D0 / BICO(X2,X1) L = - 1 GOTO 14 12 G = 1.D0 L = 1 GOTO 14 13 G = BICO(X1,X2) L = 1 14 IF (IDINT(Y)) 15,16,17 15 G = G / BICO(Y2,Y1) L = L - 2 GOTO 18 16 L = L + 2 GOTO 18 17 G = G * BICO(Y1,Y2) L = L + 2 18 IF (IDINT(Z)) 19,20,21 19 G = G * BICO(Z2,Z1) L = L - 3 GOTO 22 20 L = L + 1 GOTO 22 21 G = G / BICO(Z1,Z2) L = L + 1 22 L = (8 + L) / 2 GOTO (23,24,25,26,27,28),L 23 G = G * BICO(-2.D0*Z,-2.D0*X) GOTO 307 24 G = G / BICO(-2.D0*Y,-2.D0*Z) GOTO 307 25 G = G / BICO(-2.D0*X,-2.D0*Z) GOTO 307 26 G = G * BICO(2.D0*X,2.D0*Z) GOTO 307 27 G = G * BICO(2.D0*Y,2.D0*Z) GOTO 307 28 G = G / BICO(2.D0*Z,2.D0*X) 307 CG = CG * SUM * DSQRT(G*Z3/(BICO(R-1.D0,2.D0*C) * C BICO(2.D0*C,T) * R )) 399 RETURN END C ----------------------------------------------------------------------------- DOUBLE PRECISION FUNCTION BICO(A,B) IMPLICIT REAL*8 (A-H,O-Z) X = A - B IF (IDINT(X)) 201,202,203 201 BICO = 0.D0 GOTO 207 202 BICO = 1.D0 GOTO 207 203 IF (IDINT(B)) 201,202,204 204 Y = DMIN1(B,X) BICO = A / Y J = IDINT(Y) - 1 IF (J .LT. 1) GOTO 207 DO 205 I = 1,J U = DFLOAT(I) 205 BICO = BICO * (A-U) / (Y - U) 207 RETURN END C C ----------------------------------------------------------------------------- C SUBROUTINE RSQZ (A,AMP,IDIM,JDIM,NROW,NCOL,N,M,NORT,SQTOL) IMPLICIT REAL*8 (A-H,O-Z) C C Subroutine CSQZ modified for real matrices; see Polnaszeck PhD thesis C DIMENSION A(IDIM,JDIM),AMP(NROW,NCOL) SQTOL = 0.D0 GRASS = 0.D0 DO 101 LJM = 1,N GRASS = GRASS + 1.D0 101 SQTOL = DABS (A(LJM,1)) + SQTOL DO 102 KID = 2,M LOW = N + 1 - KID DO 102 LJM = 1,LOW GRASS = GRASS + 2.D0 102 SQTOL = 2.D0 * DABS (A(LJM,KID)) + SQTOL SQTOL = 1.D-15 * SQTOL / DSQRT (GRASS) MX = M-2 DO 20 L = 1,MX NR = M - L NQ = NR + 1 NX = N - NR IF (NX) 20,20,5 5 DO 15 K=1,NX DO 25 J = K,NX,NR IF (J-K) 40,50,40 50 IF (DABS(A(J,NQ)) - SQTOL) 15,15,60 60 B = - A(J,NR) / A(J,NQ) GOTO 80 40 IF (G**2 - SQTOL**2) 15,15,70 70 B = - A(J-1,NQ) / G 80 SX = 1.D0 / (1.D0 + B**2) S = DSQRT(SX) C = B * S CX = C**2 CS = C * S I = J + NR - 1 DO 201 IORT = 1,NORT TEMP = AMP(I,IORT) AMP(I,IORT) = C * TEMP - S * AMP(I+1,IORT) 201 AMP(I+1,IORT) = C * AMP(I+1,IORT) + S * TEMP U = A(I,1) V = A(I+1,1) A(I,1) = CX * U - 2.D0 * CS * A(I,2) + SX * V A(I+1,1) = SX * U + 2.D0 * CS * A(I,2) + CX * V A(I,2) = CS * (U-V) + (CX-SX) * A(I,2) IZ = I - 1 DO 90 IX = J,IZ IK = I - IX + 1 U = A(IX,IK) V = A(IX,IK+1) A(IX,IK) = C * U - S * V 90 A(IX,IK+1) = S * U + C * V IF(J-K) 95,85,95 95 A(J-1,NQ) = C * A(J-1,NQ) - S * G 85 IR = MIN0(NR,N-1) IF (IR-2) 55,45,45 45 DO 10 IX = 2,IR U = A(I,IX+1) A(I,IX+1) = C * U - S * A(I+1,IX) 10 A(I+1,IX) = S * U + C * A(I+1,IX) 55 IF (I+NR-N) 30,25,25 30 G = - S * A(I+1,NQ) A(I+1,NQ) = C * A(I+1,NQ) 25 CONTINUE 15 CONTINUE 20 CONTINUE RETURN END C C ----------------------------------------------------------------------------- C SUBROUTINE RQRT(A,B,AMP,IDIM,N,NROW,NORT,M,TOL,EIGVAL,IERR) IMPLICIT REAL*8 (A-H,O-Z) C C Subroutine CQRT modified for real matrices; see Polnaszeck PhD thesis C INTEGER HELL DIMENSION A(IDIM),B(IDIM),AMP(M),EIGVAL(IDIM) NX = N NI = 1 IERR = 0 TOL = 0.D0 GRASS = DFLOAT(N) DO 200 KID = 1,NX,2 200 TOL = DABS(A(KID)) + TOL + 2.D0 * DABS(B(KID)) TOL = 1.0D-15 * TOL / DSQRT(3.D0 * GRASS) C SHR will contain the total shift SHR = 0.D0 K = 0 ASSIGN 100 TO HELL IF (DABS(EIGVAL(1)) .NE. 0.) ASSIGN 101 TO HELL GOTO HELL, (100,101) 100 K = K + 1 ATR = A(NX-1) + A(NX) SC = ATR**2 - 4.D0 * (A(NX-1) * A(NX) - B(NX-1)**2) SC = DSQRT(DABS(SC)) STR = (ATR - SC) * 0.5D0 STTR = (ATR + SC) * 0.5D0 TEMP = DABS(A(NX) - STR) TEMP0 = DABS(A(NX) - STTR) IF (TEMP0 .GT. TEMP) GOTO 102 STR = STTR GOTO 102 101 K = K + 1 IF (K .GE. M) GOTO 999 STR = EIGVAL(NX) - SHR 102 IF (K .LE. M) GOTO 104 999 WRITE (6,103) M 103 FORMAT(1H0,' ***** ERROR! RQRT HAS NOT CONVERGED IN',I3, C ' ITERATIONS') IERR = 1 RETURN 104 SHR = SHR + STR DO 20 I = 1,NX 20 A(I) = A(I) - STR CR = 1.D0 CSR = 0.D0 UR = 0.D0 DO 10 I = NI,NX GR = A(I) - UR IF (I .LT. 2) GOTO 17 QR = CR * A(I)-CSR*B(I-1) GO TO 16 17 QR = CR * A(I) 16 IF (I .EQ. NX) GOTO 10 SC = QR**2 + B(I)**2 SC = DSQRT(SC) RR = SC IF (I .EQ. NI) GOTO 18 B(I-1) = SR * RR 18 RS = 1.D0 / (RR**2) SR = RR * B(I) * RS CSR = CR * SR CR = RR * QR * RS TR = GR + A(I + 1) SXR = SR**2 UR = SXR * TR A(I) = GR + UR DO 301 IORT = 1,NORT L = I + (IORT-1) * NROW TAR = AMP (L) TBR = AMP (L+1) AMP(L) = TAR * CR + TBR * SR 301 AMP(L+1) = TBR * CR - TAR * SR 10 CONTINUE B(NX-1) = SR * QR A(NX) = GR 85 IT = NX 30 IT = IT - 1 IF (DABS(B(IT)) .LE. TOL) GOTO 40 IF (IT-NI) 100,100,30 40 IF (NX-IT-2) 50,60,70 50 A(NX) = A(NX) + SHR NX = NX - 1 K = 0 GO TO 80 60 ALR = - B(NX-1) AMR = 0.5D0 * (A(NX-1) - A(NX)) SC = ALR**2 + AMR**2 SC = DSQRT(SC) ANR = SC TAR = AMR + ANR TBR = ANR - AMR SIG = 1.0D0 IF (TAR**2 .GT. TBR**2) GOTO 65 TAR = TBR SIG = - 1.0D0 65 TBR = 0.5D0 / (ANR**2) SC = ANR * TAR * TBR SC = DSQRT(SC) CR = SC TAR = ANR * CR TBR = 0.5D0 / (TAR**2) SR = SIG * TAR * ALR * TBR TAR = A(NX-1) TBR = A(NX) TCR = B(NX-1) CXR = CR**2 SXR = SR**2 CSR = CR * SR TR = CXR - SXR UR = 2.D0 * TCR * CSR A (NX-1) = TAR * CXR + TBR * SXR - UR + SHR A (NX) = TAR * SXR + TBR * CXR + UR + SHR B(NX-1) = 2.0D0 * AMR * CSR + TCR * TR I = NX - 1 DO 302 IORT = 1,NORT L = I + (IORT-1) * NROW TAR = AMP(L) TBR = AMP(L+1) AMP(L) = TAR * CR - TBR * SR 302 AMP(L+1) = TBR * CR + TAR * SR K = 0 NX = NX - 2 GOTO 80 70 NI = IT + 1 GOTO HELL,(100,101) 80 IF (NX .LT. NI) GOTO 90 95 IF (NX-NI-1) 50,60,85 90 IF (NI .EQ. 1) RETURN NI = 1 GO TO 95 END