molscat
|
README,
README.jkl,
README.v12,
dblas.f,
dblas.f.Z,
dcs_save.f,
diag_eispack.f,
ghm_save.f,
ghm_subs.f,
ghm_vib.f,
ident2disting.f,
lapack.f,
lapack.f.Z,
prbr_save.f,
prbr_vib.f,
read_isigu.f,
restrt.v12.f,
sbe.doc,
sbe_save.f,
sig_save.f,
spline.f,
syminv.f,
test1.input,
test1.v12.out,
test1.v14.out,
test2.input,
test2.v12.out,
test2.v14.out,
test3.f,
test3.input,
test3.v12.input,
test3.v12.out,
test3.v14.out,
test5.f,
test5.input,
test5.v12.out,
test5.v14.out,
test6.input,
test6.v12.out,
test6.v14.out,
test8.input,
test8.v12.out,
test8.v14.out,
timers.f,
timers_c.c,
v12.f,
v14.doc.tar,
v14.f,
v14.f.Z,
version_12.doc,
version_14.doc,
version_14.tutorial,
|
|
|
SUBROUTINE ADD0(JT1,N1,SR1,SI1,L11,L21,JT2,N2,SR2,SI2,L12,L22,
1 JIN,JFN,IQ,AR,AI)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION SR1(N1),SI1(N1),L11(N1),L21(N1),
1 SR2(N2),SI2(N2),L12(N2),L22(N2)
DATA LAMBDA/0/
FN(J1,J2)=DBLE((J1+J1+1)*(J2+J2+1))
C
AR=0.D0
AI=0.D0
DO 1000 I1=1,N1
DO 1000 I2=1,N2
IF (L11(I1).NE.L12(I2).OR.L21(I1).NE.L22(I2)) GO TO 1000
IF (L11(I1).EQ.L21(I1)) THEN
DELTA=1.D0
ELSE
DELTA=0.D0
ENDIF
CR=DELTA-(SR1(I1)*SR2(I2)+SI1(I1)*SI2(I2))
CI=SR1(I1)*SI2(I2)-SR2(I2)*SI1(I1)
XL=SQRT(FN(L11(I1),L21(I1)))
LL=L11(I1)+L21(I1)
IF (LL-2*(LL/2).NE.0) XL=-XL
XJ=X12J(JIN,JIN,L11(I1),L21(I1),JFN,L11(I1),JFN,L21(I1),
1 IQ,JT2,JT1,LAMBDA)
AR=AR+XL*XJ*CR
AI=AI+XL*XJ*CI
1000 CONTINUE
FACT=FN(JT1,JT2)
AR=AR*FACT
AI=AI*FACT
RETURN
END
SUBROUTINE ADD1(JT1,N1,SR1,SI1,L11,L21,JT2,N2,SR2,SI2,L12,L22,
1 JIN,JFN,IQ,AR,AI)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION SR1(N1),SI1(N1),L11(N1),L21(N1),
1 SR2(N2),SI2(N2),L12(N2),L22(N2)
DATA LAMBDA/1/
FN(J1,J2)=DBLE((J1+J1+1)*(J2+J2+1))
C
AR=0.D0
AI=0.D0
DO 1000 I1=1,N1
DO 1000 I2=1,N2
C SKIP ON 3J TRIANGLE RELATIONSHIPS, ELSE EVALUATE 3J'S
IF (ABS(L11(I1)-L12(I2)).GT.LAMBDA.OR.
1 ABS(L21(I1)-L22(I2)).GT.LAMBDA) GO TO 1000
TJ=THREEJ(L11(I1),L12(I2),LAMBDA)*THREEJ(L21(I1),L22(I2),LAMBDA)
IF (L11(I1).EQ.L21(I1).AND.L12(I2).EQ.L22(I2)) THEN
DELTA=1.D0
ELSE
DELTA=0.D0
ENDIF
CR=DELTA-(SR1(I1)*SR2(I2)+SI1(I1)*SI2(I2))
CI=SR1(I1)*SI2(I2)-SR2(I2)*SI1(I1)
C ADJUST FOR (-I)**(L-L'-L-BAR+L-BAR')
LL=L11(I1)-L21(I1)-L12(I2)+L22(I2)
2001 IF (LL.GT.0) GO TO 2002
LL=LL+400
GO TO 2001
2002 LL=LL-4*(LL/4)
IF (LL.EQ.1) THEN
TT=CR
CR=CI
CI=-TT
ELSEIF (LL.EQ.2) THEN
CR=-CR
CI=-CI
ELSEIF (LL.EQ.3) THEN
TT=CR
CR=-CI
CI=TT
ENDIF
XL=SQRT(FN(L11(I1),L21(I1))*FN(L12(I2),L22(I2)))
XJ=X12J(JIN,JIN,L12(I2),L22(I2),JFN,L11(I1),JFN,L21(I1),
1 IQ,JT2,JT1,LAMBDA)
AR=AR+(TJ*XL*XJ)*CR
AI=AI+(TJ*XL*XJ)*CI
1000 CONTINUE
FACT=FN(JT1,JT2)
AR=AR*FACT
AI=AI*FACT
RETURN
END
SUBROUTINE SWRITE(IU,N,S)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION S(N,N)
C
WRITE(IU) ((S(I,J),J=1,I),I=1,N)
RETURN
C
ENTRY SREAD(IU,N,S,IEND)
IEND=0
READ(IU,END=9999) ((S(I,J),J=1,I),I=1,N)
DO 1000 I=1,N
DO 1000 J=1,I-1
1000 S(J,I)=S(I,J)
RETURN
9999 IEND=1
RETURN
END
FUNCTION X12J(J1,J2,J3,J4,L1,L2,L3,L4,K1,K2,K3,K4)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C CALCULATION OF 12J SYMBOL OF 2ND KIND AS DEFINED BY
C YUTSIS ET AL., EQ. (19.3), P. 63
C CALLS J6J (MOLSCAT SHULTEN-GORDON ALGORITHM) FOR 6J'S
C VERSION 2. CHECKS FOR REPEATED 6J INDICES (SG, 22 NOV 94)
C
PARAMETER (MXDIM=100)
DIMENSION X6J1(MXDIM),X6J2(MXDIM),X6J3(MXDIM),X6J4(MXDIM)
DATA TENTH/0.1D0/
C
C CHECK FOR K4=0. THIS WILL ARISE FROM LOU MONCHICK'S GHM EQS.
C ONE MIGHT WANT TO CHECK FOR OTHER ZEROS AND REARRANGE
IF (K4.EQ.0) THEN
C USE YUTSIS EQ. (19.10)
IF (J4.NE.L4.OR.J3.NE.L2) THEN
X12J=0.D0
RETURN
ENDIF
X12J=SIXJ(L1,L2,K1,K3,K2,J1)*SIXJ(K2,L3,K3,J2,L4,K1)/
1 SQRT(DBLE((L2+L2+1)*(L4+L4+1)))
LP=J1+L1-J2-L3
IF (LP-2*(LP/2).NE.0) X12J=-X12J
RETURN
ENDIF
C
C GENERAL CASE: GET MIN, MAX VALUES OF SUMMATION INDEX
IXMIN=MAX(ABS(K1-K2),ABS(K3-K4),ABS(J1-J3),ABS(J2-J4))
IXMAX=MIN(K1+K2,K3+K4,J1+J3,J2+J4)
IF (IXMIN.GT.IXMAX) THEN
C 12J MUST VANISH IN THIS CASE
X12J=0.D0
RETURN
ENDIF
C
C CONVERT EVERYONE TO DOUBLE PRECISION FOR INPUT TO J6J
XJ1=J1
XJ2=J2
XJ3=J3
XJ4=J4
XL1=L1
XL2=L2
XL3=L3
XL4=L4
XK1=K1
XK2=K2
XK3=K3
XK4=K4
C
C GET 4 SETS OF 6J'S -- CHECK TO SEE IF SOME ARE SAME
IVAL1=MXDIM
CALL J6J(XK1,XK2,XL1,XJ3,XJ1,IVAL1,XMIN1,X6J1)
IF (K3.EQ.K1.AND.K4.EQ.K2.AND.L2.EQ.L1) THEN
IVAL2=IVAL1
XMIN2=XMIN1
DO 2001 II=1,IVAL2
2001 X6J2(II)=X6J1(II)
ELSE
IVAL2=MXDIM
CALL J6J(XK3,XK4,XL2,XJ3,XJ1,IVAL2,XMIN2,X6J2)
ENDIF
IF (L3.EQ.L1.AND.J4.EQ.J3.AND.J2.EQ.J1) THEN
IVAL3=IVAL1
XMIN3=XMIN1
DO 2002 II=1,IVAL3
2002 X6J3(II)=X6J1(II)
ELSEIF (J4.EQ.J3.AND.J2.EQ.J1.AND.L3.EQ.L2.AND.
1 K1.EQ.K3.AND.K2.EQ.K4) THEN
IVAL3=IVAL2
XMIN3=XMIN2
DO 2003 II=1,IVAL3
2003 X6J3(II)=X6J2(II)
ELSE
IVAL3=MXDIM
CALL J6J(XK1,XK2,XL3,XJ4,XJ2,IVAL3,XMIN3,X6J3)
ENDIF
IF (L4.EQ.L2.AND.J4.EQ.J3.AND.J2.EQ.J1) THEN
IVAL4=IVAL2
XMIN4=XMIN2
DO 2004 II=1,IVAL4
2004 X6J4(II)=X6J2(II)
ELSEIF (K3.EQ.K1.AND.K4.EQ.K2.AND.L4.EQ.L3) THEN
IVAL4=IVAL3
XMIN4=XMIN3
DO 2005 II=1,IVAL4
2005 X6J4(II)=X6J3(II)
ELSEIF (K3.EQ.K1.AND.K4.EQ.K2.AND.L4.EQ.L1.AND.
1 J4.EQ.J3.AND.J2.EQ.J1) THEN
IVAL4=IVAL1
XMIN4=XMIN1
DO 2006 II=1,IVAL4
2006 X6J4(II)=X6J1(II)
ELSE
IVAL4=MXDIM
CALL J6J(XK3,XK4,XL4,XJ4,XJ2,IVAL4,XMIN4,X6J4)
ENDIF
C
C SUM OVER IX=IXMIN,IXMAX
SUM=0.D0
DO 1000 IX=IXMIN,IXMAX
IND1=1+IX-INT(XMIN1+TENTH)
IF (IND1.LT.1.OR.IND1.GT.IVAL1) GO TO 1000
IND2=1+IX-INT(XMIN2+TENTH)
IF (IND2.LT.1.OR.IND2.GT.IVAL2) GO TO 1000
IND3=1+IX-INT(XMIN3+TENTH)
IF (IND3.LT.1.OR.IND3.GT.IVAL3) GO TO 1000
IND4=1+IX-INT(XMIN4+TENTH)
IF (IND4.LT.1.OR.IND4.GT.IVAL4) GO TO 1000
SUM=SUM+X6J1(IND1)*X6J2(IND2)*X6J3(IND3)*X6J4(IND4)*DBLE(IX+IX+1)
1000 CONTINUE
C GET PARITY(L1-L2-L3+L4)
LP=L1-L2-L3+L4
IF (LP-2*(LP/2).NE.0) SUM=-SUM
X12J=SUM
RETURN
END
FUNCTION SIXJ(J1,J2,J5,J4,J3,J6)
C
C CALCULATES 6-J SYMBOL: _(J1 J2 J3 )_
C (J4 J5 J6 )
C INTERFACE TO J6J ROUTINE.
C MODIFIED BY S. GREEN 20 AUG 93; PASS DIMENSION OF XJ6J FOR CHECKING
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
PARAMETER(MXDIM=200)
DIMENSION XJ6J(MXDIM)
IVAL=MXDIM
CALL J6J( DBLE(J2),DBLE(J3),
1 DBLE(J4),DBLE(J5),DBLE(J6),
3 IVAL,XJ1MIN,XJ6J)
IND=1+J1-INT(XJ1MIN+0.1D0)
SIXJ=0.D0
IF(IND.GE.1 .AND. IND.LE.IVAL) SIXJ=XJ6J(IND)
RETURN
END
SUBROUTINE J6J(J2,J3,L1,L2,L3,IVAL,J1MIN,D6J)
IMPLICIT DOUBLE PRECISION (A-H,J-Z)
DIMENSION D6J(2)
DATA ZERO/0.D0/,TENTH/0.1D0/,HALF/0.5D0/,ONE/1.D0/,TWO/2.D0/,
$ CONST/1.0D-12/
E(J1S)=SQRT((J1S-MJ23S)*(J23P1S-J1S)*(J1S-ML23S)*(L23P1S-J1S))
F(J1,JJP1)=(TWO*J1+ONE)*(JJP1*(FACT-JJP1-TWO*LLP1)+FACT2)
C
C THIS ROUTINE CALCULATES THE 6-J COEFFICIENTS FOR ALL PERMISSIBLE
C VALUES OF J1 FOR FIXED VALUES OF J2, J3, L1, L2, AND L3 USING THE
C RECURSIVE ALGORITHM OF K. SCHULTEN AND R. G. GORDON, J. MATH. PHYS.
C VOL. 16, P. 1961, (1975).
C PROGRAMMED BY D. E. FITZ, 10/22/79
C MODIFIED BY S. GREEN 20 AUG 93 TO TEST DIMENSION ON D6J
C
MXDIM=IVAL
JJP2=J2*(J2+ONE)
JJP3=J3*(J3+ONE)
LLP1=L1*(L1+ONE)
LLP2=L2*(L2+ONE)
LLP3=L3*(L3+ONE)
MJ23S=(J2-J3)**2
ML23S=(L2-L3)**2
J23P1S=(J2+J3+ONE)**2
L23P1S=(L2+L3+ONE)**2
FACT2=(LLP2-LLP3)*(JJP2-JJP3)
FACT=JJP2+JJP3+LLP2+LLP3
J1MIN=MAX(ABS(J2-J3),ABS(L2-L3))
J1MAX=MIN(J2+J3,L2+L3)
IVAL=INT(J1MAX-J1MIN+ONE+TENTH)
IF (IVAL.GT.MXDIM) THEN
WRITE(6,*) 'J6J: ARRAY D6J TOO SMALL. NEEDS ',IVAL,' BUT ONLY ',
1 MXDIM,' SUPPLIED'
STOP
ENDIF
C
C TEST FOR OTHER TRIANGULAR INEQUALITES.
C
IL1=INT(TWO*L1+TENTH)
IL2=INT(TWO*L2+TENTH)
IL3=INT(TWO*L3+TENTH)
IJ2=INT(TWO*J2+TENTH)
IJ3=INT(TWO*J3+TENTH)
IF((IJ2.LE.IL1+IL3.AND.IJ2.GE.IABS(IL1-IL3)).AND.
$ (IJ3.LE.IL1+IL2.AND.IJ3.GE.IABS(IL1-IL2))) GO TO 11
DO 12 I=1,IVAL
12 D6J(I)=ZERO
RETURN
C
11 INMID=(IVAL+3)/2
SGNV=J2+J3+L2+L3
SGN=ONE
ISIGN=INT(SGNV+TENTH)
IF(MOD(ISIGN,2).NE.0) SGN=-ONE
D6J(1)=HALF
C
C UPWARD RECURSION.
C
IF(IVAL.EQ.1) GO TO 40
JJP1=J1MIN*(J1MIN+ONE)
F1=F(J1MIN,JJP1)
J1=J1MIN+ONE
J1S=J1*J1
E2=E(J1S)
IF(J1MIN.LT.TENTH) GO TO 15
D6J(2)=-F1*D6J(1)/(E2*J1MIN)
GO TO 16
15 D6J(2)=-HALF*(LLP2+JJP2-LLP1)*D6J(1)/SQRT(JJP2*LLP2)
16 SCALE=D6J(2)
IF(IVAL.EQ.2) GO TO 40
DO 21 IJ2=3,INMID
JJP1=J1*(J1+ONE)
F1=F(J1,JJP1)
J1=J1+ONE
E1=E2
J1S=J1*J1
E2=E(J1S)
21 D6J(IJ2)=-(F1*D6J(IJ2-1)+J1*E1*D6J(IJ2-2))/(E2*(J1-ONE))
SCALE=D6J(INMID)
IEXC=5
IF(ABS(SCALE).GT.CONST) GO TO 18
INMID=INMID-1
SCALE=D6J(INMID)
IEXC=3
GO TO 30
18 IF(IVAL.EQ.3) GO TO 40
C
C DOWNWARD RECURSION.
C
30 D6J(IVAL)=HALF
J1=J1MAX
J1S=J1*J1
JJP1=J1*(J1+ONE)
F1=F(J1,JJP1)
E1=E(J1S)
D6J(IVAL-1)=-F1*D6J(IVAL)/(E1*(J1+ONE))
IEND=IVAL-INMID
IF(IVAL.LE.IEXC) GO TO 31
DO 32 IJ2=2,IEND
J1=J1-ONE
E2=E1
J1S=J1*J1
JJP1=J1*(J1+ONE)
E1=E(J1S)
F1=F(J1,JJP1)
32 D6J(IVAL-IJ2)=-(J1*E2*D6J(IVAL-IJ2+2)+F1*D6J(IVAL-IJ2+1))/
$ (E1*(J1+ONE))
C
C MATCH UPWARD AND DOWNWARD RECURSIVE RESULTS BY SCALING.
C
31 SCALE=SCALE/D6J(INMID)
DO 33 IJ2=INMID,IVAL
33 D6J(IJ2)=SCALE*D6J(IJ2)
C
C NORMALIZE RESULTS AND SET PHASE.
C
40 SUM=ZERO
DO 41 IJ2=1,IVAL
J1=J1MIN+DBLE(IJ2-1)
41 SUM=SUM+(TWO*J1+ONE)*D6J(IJ2)**2
RNORM=ONE/SQRT(SUM*(TWO*L1+ONE))
IF((SGN*D6J(IVAL)).LT.ZERO) RNORM=-RNORM
DO 42 IJ2=1,IVAL
42 D6J(IJ2)=D6J(IJ2)*RNORM
RETURN
END
FUNCTION THREEJ (J1,J2,J3)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C COMPUTATION OF SPECIAL WIGNER 3J COEFFICIENT WITH
C VANISHING PROJECTIONS. SEE EDMONDS, P. 50.
C
C THIS VERSION EVALUATES BINOM AND PARITY IN-LINE
C SHOULD IMPROVE EFFICIENCY, ESPECIALLY ON CRAY;
C ALSO GIVES IMPROVEMENT ON AMDAHL (SG: 20 DEC 92)
C
C STATEMENT FUNCTION FOR DELTA ASSOCIATED W/ RACAH AND SIXJ SYMBOLS
C DELTA(I,J,K)= SQRT(1.D0/ ( BINOM(I+J+K+1,I+J-K) *
C 1 BINOM(K+K+1,I-J+K) * DBLE(K+J-I+1) ) )
C
I1=J1+J2+J3
IF (I1-2*(I1/2).NE.0) GO TO 8
1 I2=J1-J2+J3
IF (I2) 8,2,2
2 I3=J1+J2-J3
IF (I3) 8,3,3
3 I4=-J1+J2+J3
IF (I4) 8,4,4
4 I5=I1/2
I6=I2/2
SIGN=1.D0
IF (I5-2*(I5/2).NE.0) SIGN=-SIGN
C 7 THREEJ=SIGN*DELTA(J1,J2,J3)*BINOM(I5,J1)*BINOM(J1,I6)
C B1,B2 ARE BINOM ASSOCIATED W/ DELTA
N=J1+J2+J3+1
M=J1+J2-J3
NM = N-M
MNM = MIN(NM,M)
IF(MNM.LE.0) THEN
B1=1.D0
ELSE
FN = N+1
F = 0.D0
B = 1.D0
DO 101 I = 1,MNM
F = F+1.D0
C = (FN-F)*B
101 B = C/F
B1 = B
ENDIF
N=J3+J3+1
M=J1-J2+J3
NM = N-M
MNM = MIN(NM,M)
IF(MNM.LE.0) THEN
B2=1.D0
ELSE
FN = N+1
F = 0.D0
B = 1.D0
DO 102 I = 1,MNM
F = F+1.D0
C = (FN-F)*B
102 B = C/F
B2 = B
ENDIF
DELTA=SQRT(1.D0/(B1*B2*(J3+J2-J1+1)))
C B3=BINOM(I5,J1), B4=BINOM(J1,I6)
N=I5
M=J1
NM = N-M
MNM = MIN(NM,M)
IF(MNM.LE.0) THEN
B3=1.D0
ELSE
FN = N+1
F = 0.D0
B = 1.D0
DO 103 I = 1,MNM
F = F+1.D0
C = (FN-F)*B
103 B = C/F
B3 = B
ENDIF
N=J1
M=I6
NM = N-M
MNM = MIN(NM,M)
IF(MNM.LE.0) THEN
B4=1.D0
ELSE
FN = N+1
F = 0.D0
B = 1.D0
DO 104 I = 1,MNM
F = F+1.D0
C = (FN-F)*B
104 B = C/F
B4 = B
ENDIF
THREEJ=SIGN*DELTA*B3*B4
RETURN
8 THREEJ=0.D0
RETURN
END
|