C************************************************************************
C BAND PLOT PROGRAM *
C PROGRAMMED BY H.HARIMA *
C SPIN FERRO TAIOU 1989/10/30 *
C MODIFIED BY A.YANASE *
C 1994/02/10 *
C MODIFIED BY A.YANASE 2001/01/22 *
C FOR AYPLOT SYSTEM *
C 01:SPACE GROUP AND LATTICE CONSTANTS *
C 02:ENERGY OUTPUT FILE *
C 03:FORMAT FILE *
C IFILE:PS-file by AYPLOT *
C************************************************************************
IMPLICIT REAL*8 (A-H,O-Z)
COMMON/SCL/WX,WY,EO,EM,XM,YM,IPR,ISO
CHARACTER*40 D
C
CALL TSPPRP
WRITE(6,*) ' END OF TSPPRP'
CALL TSVECT
C
CALL STRUCT(NLCOMP,NSPIN,JMARK,IPOINT,IFILE)
WRITE(6,*) ' END OF STRUCT'
write(6,*) 'jmark=',jmark
C
CALL ENERED(NLCOMP)
WRITE(6,*) ' END OF ENERED'
C
IPRR=IPR
JD=0
IF(ISO.EQ.3) JD=1
CALL AXENER(JD,IPRR,NSPIN)
WRITE(6,*) ' END OF AXENR'
C
CALL AYPSTR(IFILE)
CALL AYORIG(50.0,40.0)
write(6,*) ' END OF ORIGIN'
CALL LINEWD(1.2)
C
CALL PLOTR
write(6,*) ' END OF PLOTR'
C
IF(NSPIN.NE.2) IUD=1
IF(NSPIN.EQ.2) IUD=2
CALL LINEWD(1.1)
CALL PLOTSC(IUD)
IF(JMARK.EQ.1) CALL MARKPL(IUD,IPOINT,IFILE,NSPIN)
IF(JMARK.EQ.2) CALL MARKCN(IUD,IPOINT,IFILE,NSPIN)
IF(NSPIN.EQ.3) THEN
CALL LINEWD(0.9)
CALL SETDAS(6,4,0)
CALL PLOTSC(2)
IF(JMARK.EQ.1) CALL MARKPL(2,IPOINT,IFILE,NSPIN)
IF(JMARK.EQ.2) CALL MARKCN(2,IPOINT,IFILE,NSPIN)
END IF
write(6,*) ' END OF PLOTSC'
CALL LINEWD(2.0)
READ(3,*,END=10)EF
20 CALL LINEWD(1.2)
CALL FERMIE(EF)
write(6,*) ' END OF FERMIE'
30 READ(3,'(A40)',err=10)D
CALL TITLE(D,28)
WRITE(6,*) ' END OF TITLE'
10 CALL AYPEND
STOP
END
SUBROUTINE TITLE(D,K)
C***********************************************************************
C WRITE TITLE *
C***********************************************************************
IMPLICIT REAL*8 (A-H,O-Z)
COMMON/SCL/WX,WY,EO,EM,XM,YMM,IPR,ISO
CHARACTER*40 D
REAL*4 YM
YM=YMM
CALL AWRITE(25,D,40,2.0,YM+5.0,0)
RETURN
END
SUBROUTINE FERMIE(EF)
C***********************************************************************
C WRITE FERMI ENERGY *
C***********************************************************************
IMPLICIT REAL*8 (A-H,O-Z)
PARAMETER (NAXMAX=20)
COMMON/ISP/IX(NAXMAX),RT(NAXMAX),PS(2,NAXMAX),ICON(NAXMAX),
& MK(3,NAXMAX),NAXM
COMMON/SCL/WX,WY,EO,EM,XM,YM,IPR,ISO
REAL*4 XO,XE,YF,XMM,YMM
YMM=YM
XMM=XM
YF=(EF-EO)*WY
C write(6,*) naxm,'fermie','yf=',yf,ymm
IF(YF.GT.YMM.OR.YF.LT.0.0) RETURN
XO=0.0
DO 10 N=1,NAXM
C write(6,*) n,icon(n),naxm
IF(ICON(N).EQ.1.OR.N.EQ.NAXM) THEN
XE=WX*PS(2,N)
C write(6,*) xo,yf,xe
CALL SETDAS(6,4,1)
CALL MOVETO(XO,YF)
CALL LINETO(XE,YF)
IF(N.NE.NAXM) XO=WX*PS(1,N+1)
ENDIF
10 CONTINUE
CALL SETDAS(0,0,0)
CALL AWRITE(20,'E',1,XMM+2.0,YF-2.5,0)
CALL AWRITE(12,'F',1,XMM+7.0,YF-3.5,0)
RETURN
END
SUBROUTINE TSPPRP
C***********************************************************************
C PREPERATION TO CALL TSPACE *
C***********************************************************************
IMPLICIT REAL*8 (A-H,O-Z)
COMMON/SCL/WX,WY,EO,EM,XM,YM,IPR,ISO
DIMENSION JB(2,3)
READ(1,*)
READ(1,*) IL,NGEN,INV
CALL TSPACE(IL)
DO 1 I=1,NGEN
READ(1,*) JA,((JB(J,K),J=1,2),K=1,3)
CALL TSGENR(JA,JB)
1 CONTINUE
CALL TSPGRP(INV)
IF(IPR.GE.3) CALL TSPGDS
READ(1,104) A,B,C
READ(1,104) CA,CB,CC
104 FORMAT(3D23.16)
CALL TSLATC(A,B,C,CA,CB,CC)
RETURN
END
SUBROUTINE TSVECT
IMPLICIT REAL*8(A-H,O-Z)
COMPLEX*16 CR
COMMON/SPG2/IL,NGNG,IIGG(48),JV(2,3,48)
COMMON/SPG3/KB(3),ICB,MG,JG(48),JK(3,48)
& ,IZ,IFA(48,48),IFC,IDOUB,MTRG,JTRG(4,48),ITRC(3,48)
COMMON/SPG4/NR,NH,ND(12),IR(7,48,12),CR(48,12),IATR(12)
COMMON/SPGSPL/NIRG(12),NDIRG(12),NNRG
INTEGER TC(24),TH(12)
DATA TC/
& 1,2,2,2,3,3,3,3,3,3,3,3,
& 2,2,2,2,2,2,4,4,4,4,4,4/
DATA TH/
& 1,6,3,2,3,6,2,2,2,2,2,2/
DIMENSION KBB(3)
PAI=4.D0*DATAN(1.D0)
L=1
KBB(1)=0
KBB(2)=0
KBB(3)=0
IC=1
CALL TSIREP(KBB,IC,0)
IW=0
JRCH=1
IF(IL.LE.0) JRCH=2
C
C CHARCTER TEST
C
DO 203 IIR=1,NR
WB=0.0D0
DO 204 IGG=1,MG
IG=IIGG(JG(IGG))
INV=1
IF(JRCH.EQ.1.AND.IG.GT.24) THEN
IG=IG-24
INV=-1
END IF
IF(JRCH.EQ.2.AND.IG.GT.12) THEN
IG=IG-12
INV=-1
END IF
CL=2*L+1
IF(IG.EQ.1) GO TO 201
IF(JRCH.EQ.1) XX=(PAI/TC(IG))
IF(JRCH.EQ.2) XX=(PAI/TH(IG))
CL=SIN((2*L+1)*XX)/SIN(XX)
201 CONTINUE
CL=CL*INV
C WRITE(6,693) IGG,IG,CL,INV
693 FORMAT(2I5,F8.4,I5)
IF(ABS(CL).LT.0.001D0) GO TO 204
WB=WB+CL*CR(IGG,IIR)
204 CONTINUE
WB=WB/MG
NI=WB+0.5D0
WRITE(6,690) L,IIR,ND(IIR),NI
690 FORMAT(16I5)
NIRG(IIR)=NI
NDIRG(IIR)=ND(IIR)
203 CONTINUE
NNRG=NR
RETURN
END
SUBROUTINE STRUCT(NLCOMP,NSPIN,JMARK,IPOINT,IFILE)
C***********************************************************************
C SELECTION OF POINTS TO BE PLOTTED *
C***********************************************************************
IMPLICIT REAL*8 (A-H,O-Z)
PARAMETER (NAXMAX=20)
COMMON/ISP/IX(NAXMAX),RT(NAXMAX),PS(2,NAXMAX),ICON(NAXMAX)
& ,MK(3,NAXMAX),NAXM
COMMON/JMPAX/INDJAX(NAXMAX),RTJMP(NAXMAX),JOPT
COMMON/KKEND/KK1M(3,NAXMAX),KK2M(3,NAXMAX)
& ,ICC1M(NAXMAX),ICC2M(NAXMAX)
COMMON/SCL/WX,WY,EO,EM,XM,YM,IPR,ISO
DIMENSION KK1(3),KK2(3),KK3(3),KTV(3),KG(3)
CHARACTER*2 CCR1,CCR2,CCR3
CHARACTER*4 CHARA4
CHARACTER*10 MAGNET
C
READ(3,150) MAGNET
150 format(A10)
CHARA4=MAGNET
IF(CHARA4.EQ.'NONM') ISO=1
IF(CHARA4.EQ.'MAGN') ISO=2
IF(CHARA4.EQ.'SPIN') ISO=3
READ(3,*) NLCOMP,NSPIN,IFILE
IF((ISO.EQ.1.OR.ISO.EQ.3).AND.NSPIN.NE.0) THEN
WRITE(6,*) ' FOR PARAMAGNETIC BANDS NSPIN MUST BE 0',NSPIN
STOP
END IF
IF(ISO.EQ.2.AND.(NSPIN.LE.0.OR.NSPIN.GT.3)) THEN
WRITE(6,*) ' FOR FERROMAGNETIC BANDS'
WRITE(6,*) ' NSPIN=1 ONLY FOR UP SPIN BANDS'
WRITE(6,*) ' NSPIN=2 ONLY FOR DOWN SPIN BANDS'
WRITE(6,*) ' NSPIN=3 FOR BOTH SPIN BANDS'
STOP
END IF
READ(3,*) IPR,JMARK,IPOINT,JOPT
READ(3,*) EO,EM,YM,XM
WY=YM/(EM-EO)
C
PS(1,1)=0.0
NK=0
READ(3,*) NAXM
IF(NAXM.GT.NAXMAX) THEN
WRITE(6,*) ' NAXMAX=',NAXMAX,' NAXM=',NAXM
STOP
END IF
KK3(1)=1
KK3(2)=2
KK3(3)=3
ICC3=9999999
DO 99 NAX=1,NAXM
READ(3,*) (KK1(J),J=1,3),ICC1,(KK2(J),J=1,3),ICC2
IF(IPR.GE.2) WRITE(6,600) (KK1(J),J=1,3),ICC1,(KK2(J),J=1,3),ICC2
600 FORMAT( '(',3I4,')/',I4,'---','(',3I4,')/',I4)
C
DO 75 J=1,3
IIX=J
IF(KK1(J)*ICC2.NE.KK2(J)*ICC1) GOTO 81
75 CONTINUE
write(6,*) 'KK1 AND KK2 ARE THE SAME POINT'
STOP
81 IX(NAX)=IIX
DO 71 J=1,3
KK1M(J,NAX)=KK1(J)
KK2M(J,NAX)=KK2(J)
71 CONTINUE
ICC1M(NAX)=ICC1
ICC2M(NAX)=ICC2
KTV(1)=KK2(1)*ICC1-KK1(1)*ICC2
KTV(2)=KK2(2)*ICC1-KK1(2)*ICC2
KTV(3)=KK2(3)*ICC1-KK1(3)*ICC2
ICC=ICC1*ICC2
CALL ZZZY37(KTV(1),KTV(2),KTV(3),ICC,WW)
RT(NAX)=WW
BLK=0.0
ICON(NAX)=0
IF(NAX.NE.1) THEN
CALL EQUIKK(KK1,ICC1,KK3,ICC3,KG,IND)
IF(IND.EQ.0) THEN
ICON(NAX-1)=1
BLK=0.3*ABS(RT(1))
ENDIF
PS(1,NAX)=PS(2,NAX-1)+BLK
ENDIF
PS(2,NAX)=RT(NAX)+PS(1,NAX)
C IF(KTV(IIX).LT.0.) RT(NAX)=-RT(NAX)
IF(IPR.GE.2) WRITE(6,601) PS(1,NAX),PS(2,NAX),RT(NAX)
601 FORMAT(3F10.5)
DO 34 J=1,3
KK3(J)=KK2(J)
34 CONTINUE
ICC3=ICC2
KX=KK1(1)
KY=KK1(2)
KZ=KK1(3)
ICC=ICC1
CALL KPNAME(KX,KY,KZ,ICC,CCR1,KG)
READ(CCR1,'(A2)')MR
MK(1,NAX)=MR
KX=KK2(1)
KY=KK2(2)
KZ=KK2(3)
ICC=ICC2
CALL KPNAME(KX,KY,KZ,ICC,CCR3,KG)
READ(CCR3,'(A2)')MR
MK(3,NAX)=MR
KX=KK1(1)*ICC2+KK2(1)*ICC1
KY=KK1(2)*ICC2+KK2(2)*ICC1
KZ=KK1(3)*ICC2+KK2(3)*ICC1
ICC=ICC1*ICC2*2
CALL KPNAME(KX,KY,KZ,ICC,CCR2,KG)
READ(CCR2,'(A2)')MR
MK(2,NAX)=MR
IF(IPR.GE.2) WRITE(6,602) CCR1,CCR2,CCR3
602 FORMAT(3(3X,A2))
99 CONTINUE
WX=XM/PS(2,NAXM)
RETURN
END
SUBROUTINE PLOTR
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
C***********************************************************************
C PLOT FRAMES AND SCALES *
C***********************************************************************
IMPLICIT REAL*8 (A-H,O-Z)
PARAMETER (NAXMAX=20)
COMMON/ISP/IX(NAXMAX),RT(NAXMAX),PS(2,NAXMAX),ICON(NAXMAX),
& MK(3,NAXMAX),NAXM
COMMON/JMPAX/INDJAX(NAXMAX),RTJMP(NAXMAX),JOPT
COMMON/KKEND/KK1M(3,NAXMAX),KK2M(3,NAXMAX)
& ,ICC1M(NAXMAX),ICC2M(NAXMAX)
COMMON/SCL/WX,WY,EO,EM,XM,YMM,IPR,ISO
REAL*4 X,Y,XO,XE,YM,E
WRITE(6,*) ' NAXM=',NAXM
YM=YMM
XO=0.
DO 10 N=1,NAXM
IF(IPR.GE.3) WRITE(6,600) N,PS(1,N),PS(2,N),ICON(N)
600 FORMAT(I4,2F10.6,I2)
IF(ICON(N).EQ.1.OR.N.EQ.NAXM) THEN
XE=WX*PS(2,N)
CALL MOVETO(XO,0.0)
CALL LINETO(XE,0.0)
CALL MOVETO(XE,YM)
CALL LINETO(XO,YM)
XO=WX*PS(1,N+1)
ENDIF
10 CONTINUE
DO 301 I=1,NAXM
JM=1
IF(ICON(I).EQ.1.OR.I.EQ.NAXM) JM=2
DO 302 J=1,JM
X=PS(J,I)*WX
CALL MOVETO(X,0.0)
CALL LINETO(X,YM)
302 CONTINUE
IF(JOPT.EQ.0) GO TO 301
KX1=KK1M(1,I)
KY1=KK1M(2,I)
KZ1=KK1M(3,I)
KX2=KK2M(1,I)
KY2=KK2M(2,I)
KZ2=KK2M(3,I)
IC1=ICC1M(I)
IC2=ICC2M(I)
C write(6,*) ic1,ic2
CALL TSKFBZ(KX1,KY1,KZ1,IC1,IND1)
CALL TSKFBZ(KX2,KY2,KZ2,IC2,IND2)
IF((IND1.EQ.0.AND.IND2.NE.0).OR.(IND2.EQ.0.AND.IND1.NE.0)) THEN
IAX=I
CALL AXJUMP(IAX,IND1,IND2)
IF(INDJAX(IAX).EQ.1) THEN
X=(PS(1,I)+RTJMP(I))*WX
C write(6,*) x
CALL MOVETO(X,0.0)
CALL LINETO(X,YM)
END IF
ELSE
INDJAX(IAX)=0
END IF
301 CONTINUE
c WRITE(6,*) ' 301'
CALL MOVETO(-7.0,0.0)
CALL LINETO(-10.0,0.0)
CALL LINETO(-10.0,YM)
CALL LINETO(-7.0,YM)
ESDV=0.1
IF(EM-EO.LT.0.2) ESDV=0.02
IF(EM-EO.LT.0.1) ESDV=0.01
E0A=(EO-0.00001)/ESDV
IF(E0A.GE.0.0) I0=INT(E0A)+1
IF(E0A.LT.0.0) I0=INT(E0A)
EMA=(EM+0.00001)/ESDV
IF(EMA.GE.0.0) IM=INT(EMA)
IF(EMA.LT.0.0) IM=INT(EMA)-1
c WRITE(6,*) IM,I0,ESDV,E0A,EMA
DO 100 I=1,IM-I0+1
E=(I0+I-1)*ESDV
Y=(E-EO)*WY
IF(ABS(Y).LT.1.0) GO TO 101
IF(ABS(Y-YM).LT.1.0) GO TO 101
CALL MOVETO(-10.0,Y)
CALL LINETO(-7.0,Y)
c WRITE(6,*) I,E,Y
101 IF(ESDV.GE.0.1)
&CALL FWRITE(18,E,4,1,-13.0,Y-4.5,90)
IF(ESDV.LT.0.1)
&CALL FWRITE(18,E,4,2,-13.0,Y-4.5,90)
100 CONTINUE
C write(6,*) ' btype'
CALL AWRITE(20,'Energy (Ry)',11,-25.0,YM/2.0,90)
C write(6,*) ' btype'
DO 110 N=1,NAXM
MMM=MK(1,N)
X=WX*PS(1,N)
CALL CTYPE(X,-10.0,20,MMM)
MMM=MK(2,N)
X=WX*(PS(1,N)+PS(2,N))/2.0
CALL CTYPE(X,-7.5,18,MMM)
IF(ICON(N).EQ.1.OR.N.EQ.NAXM) THEN
MMM=MK(3,N)
X=WX*PS(2,N)
CALL CTYPE(X,-10.0,20,MMM)
ENDIF
110 CONTINUE
RETURN
END
SUBROUTINE CTYPE(X,Y,IH,M)
C***********************************************************************
C WRITE NAMES OF POINTS AND AXES *
C***********************************************************************
IMPLICIT REAL*8 (A-H,O-Z)
CHARACTER*4 D
REAL*4 X,Y,H
WRITE(D(1:2),'(A2)')M
K=2
X=X-2.5
IF(D(2:2).EQ.' ')K=1
IF(D(1:2).EQ.'GM') THEN
CALL GWRITE(IH,'G',X,Y)
ELSEIF(D(1:2).EQ.'XD') THEN
CALL CWRITE(IH,'X',1,X,Y)
ELSEIF(D(1:2).EQ.'SM') THEN
CALL GWRITE(IH,'S',X,Y)
ELSEIF(D(1:2).EQ.'LD') THEN
CALL GWRITE(IH,'L',X,Y)
ELSEIF(D(1:2).EQ.'DT') THEN
CALL GWRITE(IH,'D',X,Y)
else if(D(1:2).eq.'TP') then
CALL AWRITE(IH,"T'",2,X,Y,0)
else if(D(1:2).eq.'SP') then
CALL AWRITE(IH,"S'",2,X,Y,0)
ELSEIF(D(1:2).EQ.'MX') THEN
CALL CWRITE(IH,'M',1,X,Y)
ELSEIF(D(1:2).EQ.'LX') THEN
CALL CWRITE(IH,'L',1,X,Y)
ELSE
CALL AWRITE(IH,D,K,X,Y,0)
ENDIF
RETURN
END
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
SUBROUTINE AXJUMP(IAX,IND1,IND2)
IMPLICIT REAL*8 (A-H,O-Z)
PARAMETER (NAXMAX=20)
COMMON/ISP/IX(NAXMAX),RT(NAXMAX),PS(2,NAXMAX),ICON(NAXMAX),
& MK(3,NAXMAX),NAXM
COMMON/JMPAX/INDJAX(NAXMAX),RTJMP(NAXMAX),JOPT
COMMON/KKEND/KK1M(3,NAXMAX),KK2M(3,NAXMAX)
& ,ICC1M(NAXMAX),ICC2M(NAXMAX)
COMMON/SCL/WX,WY,EO,EM,XM,YMM,IPR,ISO
REAL*8 KB1(3),KB2(3)
DO 11 K=1,3
KB1(K)=DBLE(KK1M(K,IAX))/ICC1M(IAX)
KB2(K)=DBLE(KK2M(K,IAX))/ICC2M(IAX)
11 CONTINUE
IF(IND1.NE.0) INDP=IND1
IF(IND2.NE.0) INDP=IND2
DO 12 ICOUNT=1,100
KX=((KB1(1)+KB2(1))*0.5D0)*1000000
KY=((KB1(2)+KB2(2))*0.5D0)*1000000
KZ=((KB1(3)+KB2(3))*0.5D0)*1000000
IC=1000000
CALL TSKFBZ(KX,KY,KZ,IC,IND)
C write(6,*) ' ind',ICOUNT,IND1,IND,IND2,INDP
IF(IND.GT.INDP) GO TO 13
IF((IND.EQ.0.AND.IND1.EQ.0).OR.(IND.NE.0.AND.IND1.NE.0)) THEN
KB1(1)=DBLE(KX)/IC
KB1(2)=DBLE(KY)/IC
KB1(3)=DBLE(KZ)/IC
IND1=IND
ELSE
KB2(1)=DBLE(KX)/IC
KB2(2)=DBLE(KY)/IC
KB2(3)=DBLE(KZ)/IC
IND2=IND
END IF
KXX=(KB1(1)-KB2(1))*1000000
KYY=(KB1(2)-KB2(2))*1000000
KZZ=(KB1(3)-KB2(3))*1000000
ICC=1000000
C write(6,*) ' ICC',icc
CALL ZZZY37(KXX,KYY,KZZ,ICC,SS)
C WRITE(6,*) IAX,KX,KY,KZ,IC,SS
IF(SS.LT.1.0D-5) GO TO 13
12 CONTINUE
C WRITE(6,*) IAX,KB1,IC1,KB2,IC2
WRITE(6,*) ' STOP IN AXJUMP'
STOP
13 KKX=KX*ICC1M(IAX)-KK1M(1,IAX)*IC
KKY=KY*ICC1M(IAX)-KK1M(2,IAX)*IC
KKZ=KZ*ICC1M(IAX)-KK1M(3,IAX)*IC
IIC=IC*ICC1M(IAX)
C write(6,*) ' IIC',iic
CALL ZZZY37(KKX,KKY,KKZ,IIC,SSS)
IF(SSS.GT.0.8D-1.AND.ABS(SSS-RT(IAX)).GT.0.8D-1) THEN
INDJAX(IAX)=1
RTJMP(IAX)=SSS
ELSE
INDJAX(IAX)=0
END IF
write(6,*) ' AXJUMP',KKX,KKY,KKZ,IIC,INDJAX(IAX),SSS
RETURN
END
SUBROUTINE PLOTSC(IUD)
C***********************************************************************
C PLOT BAND *
C***********************************************************************
IMPLICIT REAL*8 (A-H,O-Z)
PARAMETER (NAXMAX=20)
PARAMETER (MAXKPT=2000,MAXEIG=120)
PARAMETER (MAXFKP=500,MAXIRP=55)
PARAMETER (N65=65)
COMMON/AXE/EAX(N65,MAXEIG,MAXIRP),NKPT(MAXIRP)
& ,NLIN(MAXIRP),JRR(MAXIRP),JUD(MAXIRP),KPM
COMMON/ISP/IX(NAXMAX),RT(NAXMAX),PS(2,NAXMAX),ICON(NAXMAX),
& MK(3,NAXMAX),NAXM
COMMON/SCL/WX,WY,EO,EM,XM,YMM,IPR,ISO
REAL*4 YM,XC,XA,YA,XB,YB
CHARACTER*4 DDD
YM=YMM
DO 110 K=1,KPM
C write(6,*) iud,k,jud(k)
IF(IUD.NE.JUD(K)) GO TO 110
NP=NKPT(K)
WRITE(DDD(1:2),'(A2)') MK(2,JRR(K))
IF(IPR.GE.1) WRITE(6,600) K,NP,DDD,NLIN(K)
600 FORMAT(2I5,1X,A2,I4)
IF(NLIN(K).LE.0) GO TO 110
DO 100 J=1,NLIN(K)
IF(IPR.GE.4) WRITE(6,*) ' BNAD NO=',J
XA=WX*PS(1,JRR(K))
YA=(EAX(1,J,K)-EO)*WY
IF(YA.GT.0.0.AND.YA.LT.YM) CALL MOVETO(XA,YA)
C IF(IPR.GE.2) WRITE(6,*) '2',XA,YA
DO 180 I=2,NP
C
XB=XA
YB=YA
XA=PS(1,JRR(K))+(RT(JRR(K))*(I-1))/(NP-1)
XA=XA*WX
YA=(EAX(I,J,K)-EO)*WY
C write(6,*) xa,ya,eax(i,j,k)
IF(YB.GT.YM) GO TO 291
IF(YA.LT.YM) GO TO 295
XC=XB+(XA-XB)*((YM-YB)/(YA-YB))
CALL LINETO(XC,YM)
GO TO 180
291 IF(YA.GT.YM) GO TO 180
XC=XB+(XA-XB)*((YB-YM)/(YB-YA))
CALL MOVETO(XC,YM)
GO TO 299
295 IF(YB.LT.0.0) GO TO 296
IF(YA.GT.0.0) GO TO 299
XC=XB+(XA-XB)*(YB/(YB-YA))
CALL LINETO(XC,0.0)
GO TO 180
296 IF(YA.LT.0.0) GO TO 180
XC=XB+(XA-XB)*(-YB/(YA-YB))
CALL MOVETO(XC,0.0)
299 CALL LINETO(XA,YA)
C WRITE(6,*) '1',XA,YA
180 CONTINUE
100 CONTINUE
110 CONTINUE
RETURN
END
SUBROUTINE ENERED(NLCOMP)
PARAMETER (MAXKPT=2000,MAXEIG=120)
IMPLICIT REAL*8 (A-H,O-Z)
CHARACTER*2 MRR,MRRM
COMMON/SCL/WX,WY,EO,EM,XM,YM,IPR,ISO
COMMON/ENR/KXM(3,MAXKPT),ICM(MAXKPT),IUDM(MAXKPT)
& ,MRRM(MAXKPT),MRNM(MAXKPT),MWEIM(MAXKPT),NSTM(MAXKPT)
& ,NEIGM(MAXKPT),EIGM(MAXEIG,MAXKPT),NKPINT
COMMON/MRK/NSIT(MAXKPT),MSIT(3,MAXKPT),XIMS(3,MAXKPT)
DIMENSION KX(3),kxx(3),EIG(MAXEIG),KBB(3,10),KG(3,10),itcr(12)
dimension ndegg(12),jtrr(12),ipart(12)
REWIND 2
C
N=0
ifrag=0
40 READ(2,150,END=41) KK,(KX(J),J=1,3),IC,IUD,MRR,MRN,MWEI,NST,NEIG
IF(NEIG.GT.MAXEIG) THEN
WRITE(6,*) ' MAXEIG=',MAXEIG,' NEIG=',NEIG
STOP
END IF
IF(NEIG.GT.0) READ(2,*) (EIG(J),J=1,NEIG)
150 FORMAT(I3,2X,3I3,1X,2I3,2X,A2,I2,3I3)
C 151 FORMAT(8F8.4)
IF(NEIG.EQ.0) GO TO 40
IF(NLCOMP.NE.0) THEN
NFACT=(NLCOMP-1)/2+1
DO 50 K=1,NEIG*NFACT
50 READ(2,*)
END IF
N=N+1
if(ifrag.eq.1) go to 40
IF(N.GT.MAXKPT) THEN
ifrag=1
go to 40
end if
CALL NEAREC(KX,IC,KBB,KG,NG)
DO 10 J=1,3
kxx(j)=kbb(j,1)
KXM(J,N)=KBB(J,1)
10 CONTINUE
icx=ic
jdob=0
if(iso.eq.3) jdob=1
call corres(kxx,icx,kx,ic,jdob,itcr,nrr,ind)
call tsirep(kxx,icx,jdob)
call dgtrst(jdubb,nrr,mmgg,nstrr,ndegg,jtrr,ipart)
C write(6,*) kxx,kx,ic,(itcr(j),j=1,nrr)
if(ind.eq.0) then
write(6,*) ' kx=',kx,ic,' kxx=',kxx,icx
write(6,*) ' stop at corres in ENERED'
stop
end if
ICM(N)=IC
C write(6,*) n,(kxm(j,n),j=1,3),icm(n),KX,IC
IUDM(N)=IUD
MRRM(N)=MRR
MRNM(N)=itcr(MRN)
if(ipart(itcr(mrn)).ne.0.and.ipart(itcr(mrn)).lt.itcr(mrn)) then
mrnm(n)=ipart(itcr(mrn))
end if
MWEIM(N)=MWEI
NSTM(N)=NST
NEIGM(N)=NEIG
DO 20 J=1,NEIG
EIGM(J,N)=EIG(J)
20 CONTINUE
GO TO 40
41 CONTINUE
if(ifrag.eq.1) then
WRITE(6,*) ' MAXKPT=',MAXKPT,' NKPINT=',N
STOP
END IF
NKPINT=N
DO 1 I=1,NKPINT
NSIT(I)=0
1 CONTINUE
RETURN
END
SUBROUTINE MARKPL(IUD,IPOINT,IFILE,NSPIN)
IMPLICIT REAL*8 (A-H,O-Z)
PARAMETER (MAXKPT=2000,MAXEIG=120)
PARAMETER (NAXMAX=20)
CHARACTER*2 MRRM
COMMON/ISP/IX(NAXMAX),RT(NAXMAX),PS(2,NAXMAX),ICON(NAXMAX),
& MK(3,NAXMAX),NAXM
COMMON/ENR/KXM(3,MAXKPT),ICM(MAXKPT),IUDM(MAXKPT)
& ,MRRM(MAXKPT),MRNM(MAXKPT),MWEIM(MAXKPT),NSTM(MAXKPT)
& ,NEIGM(MAXKPT),EIGM(MAXEIG,MAXKPT),NKPINT
COMMON/SCL/WX,WY,EO,EM,XMM,YMM,IPR,ISO
COMMON/MRK/NSIT(MAXKPT),MSIT(3,MAXKPT),XIMS(3,MAXKPT)
REAL*4 XM,YM,XA,YA
CHARACTER*1 CMARK(12)
DATA CMARK/'1','2','3','4','5','6','7','8','9','A','B','C'/
C
XM=XMM
YM=YMM
IF(IUD.EQ.1.OR.NSPIN.NE.3) write(IFILE,*) '.1 setgray'
IF(IUD.EQ.2.AND.NSPIN.EQ.3) write(IFILE,*) '.5 setgray'
IF(IUD.EQ.1.OR.NSPIN.NE.3) IPOIN=IPOINT
IF(IUD.EQ.2.AND.NSPIN.EQ.3) IPOIN=IPOINT*1.3
DO 10 I=1,NKPINT
IF(IUD.NE.IUDM(I)) GO TO 10
C write(6,*) i,mrrm(i),mrnm(i),nsit(i)
IF(NSIT(I).EQ.0) GO TO 10
DO 11 N=1,NSIT(I)
XA=(PS(1,MSIT(N,I))+RT(MSIT(N,I))*XIMS(N,I))*WX
DO 12 KK=1,NEIGM(I)
YA=(EIGM(KK,I)-EO)*WY
IF(YA.LT.YM.AND.YA.GT.0.0) THEN
CALL NRMARK(IPOIN,CMARK(MRNM(I)),XA,YA)
END IF
12 CONTINUE
11 CONTINUE
10 CONTINUE
write(IFILE,*) '0 setgray'
RETURN
END
SUBROUTINE MARKCN(IUD,IPOINT,IFILE,NSPIN)
IMPLICIT REAL*8 (A-H,O-Z)
PARAMETER (MAXKPT=2000,MAXEIG=120)
PARAMETER (NAXMAX=20)
CHARACTER*2 MRRM
COMMON/ISP/IX(NAXMAX),RT(NAXMAX),PS(2,NAXMAX),ICON(NAXMAX),
& MK(3,NAXMAX),NAXM
COMMON/ENR/KXM(3,MAXKPT),ICM(MAXKPT),IUDM(MAXKPT)
& ,MRRM(MAXKPT),MRNM(MAXKPT),MWEIM(MAXKPT),NSTM(MAXKPT)
& ,NEIGM(MAXKPT),EIGM(MAXEIG,MAXKPT),NKPINT
COMMON/SCL/WX,WY,EO,EM,XMM,YMM,IPR,ISO
COMMON/MRK/NSIT(MAXKPT),MSIT(3,MAXKPT),XIMS(3,MAXKPT)
REAL*4 XM,YM,XA,YA
C
XM=XMM
YM=YMM
IF(IUD.EQ.1.OR.NSPIN.NE.3) write(IFILE,*) '.1 setgray'
IF(IUD.EQ.2.AND.NSPIN.EQ.3) write(IFILE,*) '.5 setgray'
IF(IUD.EQ.1.OR.NSPIN.NE.3) IPOIN=IPOINT
IF(IUD.EQ.2.AND.NSPIN.EQ.3) IPOIN=IPOINT*1.3
DO 10 I=1,NKPINT
IF(IUD.NE.IUDM(I)) GO TO 10
C write(6,*) i,mrrm(i),mrnm(i),nsit(i)
IF(NSIT(I).EQ.0) GO TO 10
DO 11 N=1,NSIT(I)
XA=(PS(1,MSIT(N,I))+RT(MSIT(N,I))*XIMS(N,I))*WX
DO 12 KK=1,NEIGM(I)
YA=(EIGM(KK,I)-EO)*WY
IF(YA.LT.YM.AND.YA.GT.0.0) THEN
CALL CNMARK(IPOIN,MRNM(I),XA,YA)
END IF
12 CONTINUE
11 CONTINUE
10 CONTINUE
write(IFILE,*) '0 setgray'
RETURN
END
SUBROUTINE AXENER(JDOUB,IPR,NSPIN)
IMPLICIT REAL*8 (A-H,O-Z)
PARAMETER (NAXMAX=20)
PARAMETER (MAXKPT=2000,MAXEIG=120)
PARAMETER (MAXFKP=500,MAXIRP=55)
PARAMETER (N65=65)
CHARACTER*2 MRRM
COMMON/AXE/EAX(N65,MAXEIG,MAXIRP),NKPT(MAXIRP)
& ,NLIN(MAXIRP),JRR(MAXIRP),JUD(MAXIRP),NAXEN
COMMON/ISP/IX(NAXMAX),RT(NAXMAX),PS(2,NAXMAX),ICON(NAXMAX),
& MK(3,NAXMAX),NAXM
COMMON/KKEND/KK1M(3,NAXMAX),KK2M(3,NAXMAX)
& ,ICC1M(NAXMAX),ICC2M(NAXMAX)
COMMON/ENR/KXM(3,MAXKPT),ICM(MAXKPT),IUDM(MAXKPT)
& ,MRRM(MAXKPT),MRNM(MAXKPT),MWEIM(MAXKPT),NSTM(MAXKPT)
& ,NEIGM(MAXKPT),EIGM(MAXEIG,MAXKPT),NKPINT
DIMENSION KB(3),JTR(12),IPA(12),ND(12)
DIMENSION NRECPT(3,20)
JD=JDOUB
NAXEN=0
DO 11 I=1,NAXM
II=I
IF(IPR.GE.2) WRITE(6,601) II
601 FORMAT(' RECFND',I3)
CALL RECFND(II,NRECPT,NRPT,IPR)
IF(IPR.GE.2) WRITE(6,602) II,NRPT
602 FORMAT(' KPFIEN',I3,' NUMBER OF RECPR. LAT. POINTS=',I3)
CALL KPFIEN(II,NRECPT,NRPT,IPR)
DO 21 K=1,3
KB(K)=KK1M(K,I)*ICC2M(I)+KK2M(K,I)*ICC1M(I)
21 CONTINUE
IC=ICC1M(I)*ICC2M(I)*2
CALL TSIREP(KB,IC,JD)
CALL DGTRMD(JDO,NR,NH,NSTR,ND,MTR,JTR,IPA)
IF(IPR.GE.2) WRITE(6,603) KB,IC
603 FORMAT(' CMPTBL (',3I3,')/',I4)
CALL CMPTBL(JD,II,KB,IC,NR,NH,MTR,IPR)
DO 22 IR=1,NR
IF(JTR(IR).EQ.0.AND.IPA(IR).LT.IR) GO TO 22
IIR=IR
NAXEN=NAXEN+1
IF(NAXEN.GT.MAXIRP) THEN
WRITE(6,*) ' MAXIRP=',MAXIRP,' NAXEN=',NAXEN
STOP
END IF
IF(NSPIN.NE.2) IUD=1
IF(NSPIN.EQ.2) IUD=2
CALL INTPOR(JD,II,KB,IC,NH,IIR,IUD,IPR)
JUD(NAXEN)=IUD
IF(NSPIN.EQ.3) THEN
NAXEN=NAXEN+1
IF(NAXEN.GT.MAXIRP) THEN
WRITE(6,*) ' MAXIRP=',MAXIRP,' NAXEN=',NAXEN
STOP
END IF
IUD=2
CALL INTPOR(JD,II,KB,IC,NH,IIR,IUD,IPR)
JUD(NAXEN)=IUD
END IF
22 CONTINUE
11 CONTINUE
RETURN
END
SUBROUTINE RECFND(II,NRECPT,NRPT,IPR)
IMPLICIT REAL*8 (A-H,O-Z)
PARAMETER (NAXMAX=20)
C
COMMON/KKEND/KK1M(3,NAXMAX),KK2M(3,NAXMAX)
& ,ICC1M(NAXMAX),ICC2M(NAXMAX)
C
DIMENSION KB(3),NRECPT(3,20),KG(3,10),KBB(3,10)
JJ=0
DO 10 J=1,4
IF(J.EQ.1) THEN
DO 11 K=1,3
KB(K)=2*KK1M(K,II)*ICC2M(II)-KK2M(K,II)*ICC1M(II)
11 CONTINUE
ICC=ICC2M(II)*ICC1M(II)
ELSE IF(J.EQ.2) THEN
DO 12 K=1,3
KB(K)=KK1M(K,II)
12 CONTINUE
ICC=ICC1M(II)
ELSE IF(J.EQ.3) THEN
DO 13 K=1,3
KB(K)=KK2M(K,II)
13 CONTINUE
ICC=ICC2M(II)
ELSE IF(J.EQ.4) THEN
DO 14 K=1,3
KB(K)=2*KK2M(K,II)*ICC1M(II)-KK1M(K,II)*ICC2M(II)
14 CONTINUE
ICC=ICC1M(II)*ICC2M(II)
END IF
CALL NEAREC(KB,ICC,KBB,KG,NNG)
DO 16 N=1,NNG
IF(JJ.EQ.0) GO TO 17
DO 18 IJ=1,JJ
DO 19 K=1,3
IF(KG(K,N).NE.NRECPT(K,IJ)) GO TO 18
19 CONTINUE
GO TO 16
18 CONTINUE
17 JJ=JJ+1
DO 15 K=1,3
NRECPT(K,JJ)=KG(K,N)
15 CONTINUE
IF(IPR.GE.4) WRITE(6,601) II,J,JJ,(NRECPT(K,JJ),K=1,3)
601 FORMAT(6I5)
16 CONTINUE
10 CONTINUE
NRPT=JJ
RETURN
END
SUBROUTINE KPFIEN(II,NRECPT,NRPT,IPR)
IMPLICIT REAL*8 (A-H,O-Z)
PARAMETER (NAXMAX=20)
PARAMETER (MAXKPT=2000,MAXEIG=120)
PARAMETER (MAXFKP=500,MAXIRP=55)
CHARACTER*2 MRRM
C
COMPLEX*16 CW
COMMON/SPG2/IL,NG,IG(48),JV(2,3,48)
COMMON/STK/KS(3,48),JS(48),NS,ICBB,CW(48,12)
C
COMMON/KKEND/KK1M(3,NAXMAX),KK2M(3,NAXMAX)
& ,ICC1M(NAXMAX),ICC2M(NAXMAX)
COMMON/ISP/IX(NAXMAX),RT(NAXMAX),PS(2,NAXMAX),ICON(NAXMAX),
& MK(3,NAXMAX),NAXM
COMMON/ENR/KXM(3,MAXKPT),ICM(MAXKPT),IUDM(MAXKPT)
& ,MRRM(MAXKPT),MRNM(MAXKPT),MWEIM(MAXKPT),NSTM(MAXKPT)
& ,NEIGM(MAXKPT),EIGM(MAXEIG,MAXKPT),NKPINT
COMMON/AXP/IFKP(MAXFKP),JSKP(MAXFKP),KKB(3,MAXFKP)
& ,ICCB(MAXFKP),ICPTBL(12,MAXFKP),INDM(MAXFKP),NPOINT
C
DIMENSION KB(3),KBC(3)
DIMENSION NRECPT(3,20),JIND(MAXKPT),IMM(MAXKPT)
C
CALL ADDINV(NTLATC)
X1=DBLE(KK1M(1,II))/ICC1M(II)
Y1=DBLE(KK1M(2,II))/ICC1M(II)
Z1=DBLE(KK1M(3,II))/ICC1M(II)
DX=DBLE(KK2M(1,II))/ICC2M(II)-X1
DY=DBLE(KK2M(2,II))/ICC2M(II)-Y1
DZ=DBLE(KK2M(3,II))/ICC2M(II)-Z1
C write(6,*) x1,y1,z1
C write(6,*) dx,dy,dz
D6=1.0D-6
NPOINT=0
DO 10 I=1,NKPINT
JIND(I)=0
10 CONTINUE
DO 1 I=1,NKPINT
IF(JIND(I).NE.0) GO TO 1
KB(1)=KXM(1,I)
KB(2)=KXM(2,I)
KB(3)=KXM(3,I)
IC=ICM(I)
c write(6,*) 'kb',kb,ic
NJ=1
IMM(1)=I
IF(I.LE.NKPINT) THEN
DO 4 JI=I+1,NKPINT
IF(JIND(JI).NE.0) GO TO 4
DO 5 K=1,3
IF(KB(K)*ICM(JI).NE.KXM(K,JI)*IC) GO TO 4
5 CONTINUE
NJ=NJ+1
IMM(NJ)=JI
4 CONTINUE
END IF
C-----------------------------------------
C STAR OF K IS OBTAINED
CALL TSIREP(KB,IC,0)
CALL ZZZY38
C------------------------------------------
DO 2 J=1,NS
c write(6,*) j,(ks(k,j),k=1,3),icbb
DO 3 KGG=1,NRPT
XK=DBLE(KS(1,J))/ICBB+DBLE(NRECPT(1,KGG))
YK=DBLE(KS(2,J))/ICBB+DBLE(NRECPT(2,KGG))
ZK=DBLE(KS(3,J))/ICBB+DBLE(NRECPT(3,KGG))
c write(6,*) xk,yk,zk
IF(DABS(DX).LE.D6.AND.DABS(XK-X1).GT.D6) GO TO 3
IF(DABS(DY).LE.D6.AND.DABS(YK-Y1).GT.D6) GO TO 3
IF(DABS(DZ).LE.D6.AND.DABS(ZK-Z1).GT.D6) GO TO 3
IF(DABS(DX).GT.D6) THEN
IF(DABS(DY).GT.D6) THEN
IF(DABS((XK-X1)/DX-(YK-Y1)/DY).GT.D6) GO TO 3
END IF
IF(DABS(DZ).GT.D6) THEN
IF(DABS((XK-X1)/DX-(ZK-Z1)/DZ).GT.D6) GO TO 3
END IF
END IF
IF(DABS(DY).GT.D6.AND.DABS(DZ).GT.D6) THEN
IF(DABS((YK-Y1)/DY-(ZK-Z1)/DZ).GT.D6) GO TO 3
END IF
IF(IX(II).EQ.1.AND.
& ((XK-X1)/DX.LT.-1.0D0.OR.(XK-X1)/DX.GT.2.0D0)) GO TO 3
IF(IX(II).EQ.2.AND.
& ((YK-Y1)/DY.LT.-1.0D0.OR.(YK-Y1)/DY.GT.2.0D0)) GO TO 3
IF(IX(II).EQ.3.AND.
& ((ZK-Z1)/DZ.LT.-1.0D0.OR.(ZK-Z1)/DZ.GT.2.0D0)) GO TO 3
KBC(1)=KS(1,J)+NRECPT(1,KGG)*ICBB
KBC(2)=KS(2,J)+NRECPT(2,KGG)*ICBB
KBC(3)=KS(3,J)+NRECPT(3,KGG)*ICBB
IF(NPOINT.EQ.0) GO TO 20
DO 21 N=1,NPOINT
DO 22 K=1,3
IF(KBC(K)*ICCB(N).NE.KKB(K,N)*ICBB) GO TO 21
22 CONTINUE
IF(MRNM(IFKP(N)).EQ.MRNM(I)) GO TO 3
21 CONTINUE
20 DO 6 JI=1,NJ
JIND(IMM(JI))=I
NPOINT=NPOINT+1
IF(NPOINT.GT.MAXFKP) THEN
WRITE(6,*) ' MAXFKP=',MAXFKP,' NPOINT=',NPOINT
STOP
END IF
IFKP(NPOINT)=IMM(JI)
JSKP(NPOINT)=JS(J)
KKB(1,NPOINT)=KBC(1)
KKB(2,NPOINT)=KBC(2)
KKB(3,NPOINT)=KBC(3)
ICCB(NPOINT)=ICBB
IF(IPR.GE.4) THEN
N=NPOINT
WRITE(6,601) N,KBC,ICBB,IFKP(N)
& ,KB,IC,MRRM(IFKP(N)),MRNM(IFKP(N))
& ,j,kgg,(KS(k,J),K=1,3)
601 FORMAT(I3,'(',3I4,')/',I4,I3,
& '(',3I3,')/',I3,2X,A2,I3,' J=',I3,' KGG='
& ,I3,'KS=',3i4)
END IF
6 CONTINUE
c write(6,*) xk,yk,zk
c write(6,*) npoint,kbc,icbb,i,j,kgg
3 CONTINUE
2 CONTINUE
1 CONTINUE
CALL REMINV
RETURN
END
SUBROUTINE CMPTBL(JD,II,KBB,ICC,NRT,NH,MTR,IPR)
IMPLICIT REAL*8 (A-H,O-Z)
PARAMETER (NAXMAX=20)
PARAMETER (MAXKPT=2000,MAXEIG=120)
PARAMETER (MAXFKP=500,MAXIRP=55)
PARAMETER (N65=65)
CHARACTER*2 MRRM
COMMON/AXE/EAX(N65,MAXEIG,MAXIRP),NKPT(MAXIRP)
& ,NLIN(MAXIRP),JRR(MAXIRP),JUD(MAXIRP),NAXEN
COMMON/AXP/IFKP(MAXFKP),JSKP(MAXFKP),KKB(3,MAXFKP)
& ,ICCB(MAXFKP),ICPTBL(12,MAXFKP),INDM(MAXFKP),NPOINT
COMMON/ENR/KXM(3,MAXKPT),ICM(MAXKPT),IUDM(MAXKPT)
& ,MRRM(MAXKPT),MRNM(MAXKPT),MWEIM(MAXKPT),NSTM(MAXKPT)
& ,NEIGM(MAXKPT),EIGM(MAXEIG,MAXKPT),NKPINT
C
DIMENSION KB(3),KBX(3),KBB(3)
DIMENSION ITCR(12),ND(12),JTR(12),IPA(12)
DIMENSION NDX(12),JTRX(12),IPAX(12)
DIMENSION ICP(12,12),ICP0(12,12)
C
JDOB=JD
DO 11 I=1,NPOINT
INDM(I)=0
11 CONTINUE
DO 1 I=1,NPOINT
IF(INDM(I).NE.0) GO TO 1
KB(1)=KKB(1,I)
KB(2)=KKB(2,I)
KB(3)=KKB(3,I)
IC=ICCB(I)
KBX(1)=KXM(1,IFKP(I))
KBX(2)=KXM(2,IFKP(I))
KBX(3)=KXM(3,IFKP(I))
ICX=ICM(IFKP(I))
CALL CORRES(KB,IC,KBX,ICX,JDOB,ITCR,NRR,IND)
IF(IND.EQ.0) THEN
write(6,*) ' kb=',kb,ic,' kbx=',kbx,icx
WRITE(6,*) ' STOP AT CORRES IN CMPTBL'
STOP
END IF
CALL TSIREP(KBX,ICX,JDOB)
CALL DGTRMD(JDO,NRX,NHX,NSTRX,NDX,MTRX,JTRX,IPAX)
CALL TSIREP(KB,IC,JDOB)
CALL DGTRMD(JDO,NRR,NHH,NSTR,ND,MTRR,JTR,IPA)
DO 4 IR1=1,NRX
IF(JTR(IR1).EQ.0.AND.IR1.LT.IPAX(IR1).AND.
& ITCR(IR1).GT.IPA(ITCR(IR1))) THEN
IW1=ITCR(IR1)
ITCR(IR1)=ITCR(IPAX(IR1))
ITCR(IPAX(IR1))=IW1
END IF
4 CONTINUE
C write(6,*) ' cor',kb,ic,kbx,icx,(itcr(L),L=1,nrr)
IF(NH.NE.NHH.OR.MTR.NE.MTRR) THEN
CALL COMPAT(KBB,ICC,NRTT,KB,IC,NRR,JDOB,ICP,INDC)
C write(6,*) kbb,icc,kb,ic
C do 1000 i1=1,nrt
C write(6,*) ' icp',(icp(L,i1),L=1,nrr)
C 1000 continue
IF(INDC.NE.0) THEN
WRITE(6,*) 'STOP AT COMPAT IN COMTBL'
STOP
END IF
CALL CMPTRV(KBB,ICC,KB,IC,JDOB,ICP,ICP0)
C do 1001 i1=1,nrt
C write(6,*) ' icp0',(icp0(L,i1),L=1,nrr)
C 1001 continue
DO 2 IIR=1,NRT
ICPTBL(IIR,I)=ICP0(ITCR(MRNM(IFKP(I))),IIR)
2 CONTINUE
INDM(I)=1
IF(IPR.GE.4) WRITE(6,601) I,INDM(I),(ICPTBL(K,I),K=1,NRT)
601 FORMAT(14I3)
IF(I.EQ.NPOINT) GO TO 1
DO 20 J=I+1,NPOINT
DO 22 K=1,3
IF(KB(K)*ICCB(J).NE.KKB(K,J)*IC) GO TO 1
22 CONTINUE
DO 23 IIR=1,NRT
ICPTBL(IIR,J)=ICP0(ITCR(MRNM(IFKP(J))),IIR)
23 CONTINUE
INDM(J)=1
IF(IPR.GE.4)
& WRITE(6,601) J,INDM(J),(ICPTBL(K,J),K=1,NRT)
20 CONTINUE
ELSE
DO 3 IIR=1,NRT
ICPTBL(IIR,I)=0
3 CONTINUE
ICPTBL(ITCR(MRNM(IFKP(I))),I)=1
INDM(I)=-1
IF(IPR.GE.4)
& WRITE(6,601) I,INDM(I),(ICPTBL(K,I),K=1,NRT)
IF(I.EQ.NPOINT) GO TO 1
DO 30 J=I+1,NPOINT
DO 32 K=1,3
IF(KB(K)*ICCB(J).NE.KKB(K,J)*IC) GO TO 1
32 CONTINUE
DO 33 IIR=1,NRT
ICPTBL(IIR,J)=0
33 CONTINUE
ICPTBL(ITCR(MRNM(IFKP(J))),J)=1
INDM(J)=-1
IF(IPR.GE.4)
& WRITE(6,601) J,INDM(J),(ICPTBL(K,J),K=1,NRT)
30 CONTINUE
END IF
1 CONTINUE
RETURN
END
SUBROUTINE INTPOR(JD,II,KBB,ICC,NH,IIR,IUD,IPR)
IMPLICIT REAL*8 (A-H,O-Z)
PARAMETER (NAXMAX=20)
PARAMETER (MAXKPT=2000,MAXEIG=120)
PARAMETER (MAXFKP=500,MAXIRP=55)
PARAMETER (N65=65)
CHARACTER*2 MRRM,MRR
COMMON/SPGSPL/NIRG(12),NDIRG(12),NNRG
COMMON/AXE/EAX(N65,MAXEIG,MAXIRP),NKPT(MAXIRP)
& ,NLIN(MAXIRP),JRR(MAXIRP),JUD(MAXIRP),NAXEN
COMMON/AXP/IFKP(MAXFKP),JSKP(MAXFKP),KKB(3,MAXFKP)
& ,ICCB(MAXFKP),ICPTBL(12,MAXFKP),INDM(MAXFKP),NPOINT
COMMON/ISP/IX(NAXMAX),RT(NAXMAX),PS(2,NAXMAX),ICON(NAXMAX),
& MK(3,NAXMAX),NAXM
COMMON/KKEND/KK1M(3,NAXMAX),KK2M(3,NAXMAX)
& ,ICC1M(NAXMAX),ICC2M(NAXMAX)
COMMON/ENR/KXM(3,MAXKPT),ICM(MAXKPT),IUDM(MAXKPT)
& ,MRRM(MAXKPT),MRNM(MAXKPT),MWEIM(MAXKPT),NSTM(MAXKPT)
& ,NEIGM(MAXKPT),EIGM(MAXEIG,MAXKPT),NKPINT
COMMON/MRK/NSIT(MAXKPT),MSIT(3,MAXKPT),XIMS(3,MAXKPT)
C
DIMENSION KBB(3),IRNM(MAXEIG,MAXFKP),J00LR(MAXEIG,2)
DIMENSION XX(N65),YY(N65),X(MAXFKP),YM(MAXEIG,MAXFKP)
DIMENSION SM(MAXFKP),JMAX(MAXFKP)
DIMENSION IXM(MAXFKP),Y(MAXFKP),XXI(MAXFKP)
DIMENSION NL(MAXEIG),NR(MAXEIG)
C
D6=1.0D-6
IF(IPR.GE.2) THEN
KX=KBB(1)
KY=KBB(2)
KZ=KBB(3)
CALL KPNAME(KX,KY,KZ,ICC,MRR,KG)
WRITE(6,600) MRR,IIR
600 FORMAT(' INTERPORATION START FOR ',A2,I2)
END IF
JRR(NAXEN)=II
JDOB=JD
WIDE=DBLE(KK2M(IX(II),II))/ICC2M(II)
& -DBLE(KK1M(IX(II),II))/ICC1M(II)
NXI=0
DO 1 I=1,NPOINT
IF(IUD.NE.IUDM(IFKP(I))) GO TO 1
XI=(DBLE(KKB(IX(II),I))/ICCB(I)
& -DBLE(KK1M(IX(II),II))/ICC1M(II))/WIDE
C write(6,*) i,xi,nxi
IF(IIR.NE.1) GO TO 13
IF(XI+D6.LT.0.0) GO TO 13
IF(XI-D6.GT.1.0D0) GO TO 13
IF(XI+D6.LT.1.0D0.OR.II.EQ.NAXM.OR.ICON(II).EQ.1) THEN
NSIT(IFKP(I))=NSIT(IFKP(I))+1
IF(NSIT(IFKP(I)).GT.3) THEN
WRITE(6,*) ' NSIT MUST BE LESS THAN 3',NSIT(IFKP(I))
STOP
END IF
MSIT(NSIT(IFKP(I)),IFKP(I))=II
XIMS(NSIT(IFKP(I)),IFKP(I))=XI
END IF
C write(6,*) i,ifkp(i),nsit(ifkp(i)),ii,naxm,xi,'13up'
13 CONTINUE
IF(ICPTBL(IIR,I).EQ.0) GO TO 1
IF(NXI.EQ.0) GO TO 11
DO 15 NI=1,NXI
IF(DABS(XI-X(NI)).LT.1.0D-6) GO TO 12
15 CONTINUE
11 NXI=NXI+1
JMAX(NXI)=0
X(NXI)=XI
NI=NXI
12 CONTINUE
DO 5 L=1,NEIGM(IFKP(I))
DO 2 K=1,ICPTBL(IIR,I)
IF(JMAX(NI).EQ.0) GO TO 51
DO 52 JJ=1,JMAX(NI)
J1=JMAX(NI)-JJ+1
IF(YM(J1,NI).LT.EIGM(L,IFKP(I))) GO TO 53
YM(J1+1,NI)=YM(J1,NI)
IRNM(J1+1,NI)=IRNM(J1,NI)
52 CONTINUE
51 J1=0
53 YM(J1+1,NI)=EIGM(L,IFKP(I))
IRNM(J1+1,NI)=MRNM(IFKP(I))
JMAX(NI)=JMAX(NI)+1
2 CONTINUE
5 CONTINUE
1 CONTINUE
DO 6 I=1,NXI
IXM(I)=I
6 CONTINUE
XXI(1)=X(1)
DO 61 I=2,NXI
DO 62 J=1,I-1
JJ=I-J
IF(XXI(JJ).LT.X(I)) GO TO 63
XXI(JJ+1)=XXI(JJ)
IXM(JJ+1)=IXM(JJ)
62 CONTINUE
JJ=0
63 JJ=JJ+1
XXI(JJ)=X(I)
IXM(JJ)=I
61 CONTINUE
C write(6,*) ' ixm',(ixm(i),i=1,nxi)
C write(6,*) ' xxi',(xxi(i),i=1,nxi)
C write(6,*) II,IIR
CALL CHGORD(XXI,IXM,NXI,YM,JMAX,JL,NL,JR,NR)
NLIN(NAXEN)=0
NKPT(NAXEN)=N65
JJMAX=0
DO 21 I=1,NXI
IF(JMAX(I).GT.JJMAX) JJMAX=JMAX(I)
21 CONTINUE
CALL ENDTRM(II,XXI,IXM,NXI,IRNM,JMAX,JJMAX,J00LR)
C write(6,*) 'jjmax=',jjmax
C write(6,*) (jmax(ixm(i)),i=1,nxi)
DO 70 J=1,JJMAX
NNN=0
DO 71 I=1,NXI
IF(JMAX(IXM(I)).LT.J) GO TO 71
IF(XXI(I).LT.0.0D0) THEN
IF(J00LR(J,1).EQ.0.OR.J00LR(J,1).EQ.9) GO TO 71
NNN=NNN+1
X(NNN)=XXI(I)
IF(JL.EQ.1) THEN
Y(NNN)=YM(NL(J),IXM(I))
ELSE
Y(NNN)=YM(J,IXM(I))
END IF
ELSE IF(XXI(I).GT.1.0D0) THEN
IF(J00LR(J,2).EQ.0.OR.J00LR(J,2).EQ.9) GO TO 71
NNN=NNN+1
X(NNN)=XXI(I)
IF(JR.EQ.1) THEN
Y(NNN)=YM(NR(J),IXM(I))
ELSE
Y(NNN)=YM(J,IXM(I))
END IF
ELSE
NNN=NNN+1
X(NNN)=XXI(I)
Y(NNN)=YM(J,IXM(I))
END IF
71 CONTINUE
IF(NNN.LT.3) GO TO 70
IF(IPR.GE.4) THEN
WRITE(6,*) ' BAND NO=',J
WRITE(6,666) (X(I),I=1,NNN)
WRITE(6,666) (Y(I),I=1,NNN)
666 FORMAT(10F7.4)
END IF
IF(J00LR(J,1).NE.0) THEN
C1=4.0*((Y(2)-Y(1))/(X(2)-X(1)))
AMU1=2.0D0
ELSE
C1=0.0D0
AMU1=0.0D0
END IF
IF(J00LR(J,2).NE.0) THEN
CN=4.0*((Y(NNN)-Y(NNN-1))/(X(NNN)-X(NNN-1)))
ALMN=2.0D0
ELSE
CN=0.0D0
ALMN=0.0D0
END IF
INIT=0
ICO=0
DO 100 I=1,N65
DD=DBLE(I-1)/(N65-1)
IF(DD.LT.X(1)) GO TO 100
IF(DD.GT.X(NNN)) GO TO 101
IF(INIT.EQ.0) INIT=I
ICO=ICO+1
XX(ICO)=DD
100 CONTINUE
101 CONTINUE
IF(ICO.LE.3) GO TO 70
NLIN(NAXEN)=J
CALL S3N(X,Y,SM,XX,YY,NNN,ICO,C1,CN,AMU1,ALMN)
DO 72 I=1,N65
EAX(I,J,NAXEN)=99.0
72 CONTINUE
DO 75 I=1,ICO
EAX(INIT+I-1,J,NAXEN)=YY(I)
75 CONTINUE
70 CONTINUE
RETURN
END
SUBROUTINE CHGORD(XXI,IXM,NXI,YM,JMAX,JL,NL,JR,NR)
IMPLICIT REAL*8 (A-H,O-Z)
PARAMETER (MAXEIG=120)
PARAMETER (MAXFKP=500)
C
DIMENSION KKB(3),KBB(3),KG(3),NL(MAXEIG),NR(MAXEIG)
DIMENSION YM(MAXEIG,MAXFKP),XXI(MAXFKP),IXM(MAXFKP),JMAX(MAXFKP)
C
D6=1.0D-6
DO 1 IX=1,NXI
IF(ABS(XXI(IX)).LT.D6) GO TO 2
1 CONTINUE
JL=0
GO TO 31
2 I=IXM(IX)
JL=1
NS=1
DO 3 J=2,JMAX(I)
IF(YM(J,I)-YM(J-1,I).GT.D6) THEN
DO 4 K=1,NS
NL(J-NS-1+K)=J-K
4 CONTINUE
NS=1
ELSE
NS=NS+1
END IF
3 CONTINUE
DO 5 K=1,NS
NL(JMAX(I)-NS+K)=JMAX(I)+1-K
5 CONTINUE
C WRITE(6,*) IX,I,(NL(K),K=1,JMAX(I))
31 CONTINUE
DO 11 IX=1,NXI
IF(ABS(XXI(IX)-1.0D0).LT.D6) GO TO 12
11 CONTINUE
JR=0
GO TO 32
12 I=IXM(IX)
JR=1
NS=1
DO 13 J=2,JMAX(I)
IF(YM(J,I)-YM(J-1,I).GT.D6) THEN
DO 14 K=1,NS
NR(J-NS-1+K)=J-K
14 CONTINUE
NS=1
ELSE
NS=NS+1
END IF
13 CONTINUE
DO 15 K=1,NS
NR(JMAX(I)-NS+K)=JMAX(I)+1-K
15 CONTINUE
C WRITE(6,*) IX,I,(NR(K),K=1,JMAX(I))
32 CONTINUE
RETURN
END
SUBROUTINE ENDTRM(IAX,XXI,IXM,NXI,IRNM,JMAX,JJMAX,J00LR)
IMPLICIT REAL*8 (A-H,O-Z)
COMPLEX*16 DB(6,3,6,5)
PARAMETER (NAXMAX=20)
PARAMETER (MAXKPT=2000,MAXEIG=120)
PARAMETER (MAXFKP=500,MAXIRP=55)
COMMON/SPGSPL/NIRG(12),NDIRG(12),NNRG
COMMON/ISP/IX(NAXMAX),RT(NAXMAX),PS(2,NAXMAX),ICON(NAXMAX)
& ,MK(3,NAXMAX),NAXM
COMMON/KKEND/KK1M(3,NAXMAX),KK2M(3,NAXMAX)
& ,ICC1M(NAXMAX),ICC2M(NAXMAX)
COMMON/SCL/WX,WY,EO,EM,XM,YM,IPR,ISO
DIMENSION IXM(MAXFKP),XXI(MAXFKP),IRNM(MAXEIG,MAXFKP)
& ,J00LR(MAXEIG,2),JMAX(MAXFKP)
DIMENSION KB(3),KBI(3),KBO(3),KG(3),JNDLR(12,2)
DIMENSION NDES(12),JTRS(12),IPAS(12),ITCR(12)
D6=1.0D-6
INDUB=0
IF(ISO.EQ.3) INDUB=1
DO 13 JLR=1,2
IF(JLR.EQ.1) THEN
DO 1 I=1,NXI
IF(ABS(XXI(I)).LT.D6) GO TO 2
1 CONTINUE
DO 31 J=1,JJMAX
J00LR(J,1)=999
31 CONTINUE
GO TO 13
2 ILEFT=IXM(I)
DO 12 K=1,3
KBI(K)=KK1M(K,IAX)*ICC2M(IAX)*9+KK2M(K,IAX)*ICC1M(IAX)
KBO(K)=KK1M(K,IAX)*ICC2M(IAX)*11-KK2M(K,IAX)*ICC1M(IAX)
12 CONTINUE
ELSE
DO 11 I=1,NXI
IF(ABS(XXI(I)-1.0D0).LT.D6) GO TO 42
11 CONTINUE
DO 32 J=1,JJMAX
J00LR(J,2)=999
32 CONTINUE
GO TO 13
42 IRIGHT=IXM(I)
DO 14 K=1,3
KBI(K)=KK1M(K,IAX)*ICC2M(IAX)+KK2M(K,IAX)*ICC1M(IAX)*9
KBO(K)=-KK1M(K,IAX)*ICC2M(IAX)+KK2M(K,IAX)*ICC1M(IAX)*11
14 CONTINUE
END IF
ICC=10*ICC1M(IAX)*ICC2M(IAX)
ICO=ICC
C WRITE(6,*) KB,ICC
C WRITE(6,*) KBO,ICO
CALL EQUIKK(KBI,ICC,KBO,ICO,KG,INDEQK)
C WRITE(6,*) KG,IND
CALL CORRES(KBI,ICC,KBO,ICO,INDUB,ITCR,NRRR,IND)
C WRITE(6,*) ' CORRES',(ITCR(I),I=1,NRRR)
IF(JLR.EQ.1) THEN
KX=KK1M(1,IAX)
KY=KK1M(2,IAX)
KZ=KK1M(3,IAX)
ICC=ICC1M(IAX)
ELSE
KX=KK2M(1,IAX)
KY=KK2M(2,IAX)
KZ=KK2M(3,IAX)
ICC=ICC2M(IAX)
END IF
KB(1)=KX
KB(2)=KY
KB(3)=KZ
CALL TSIREP(KB,ICC,INDUB)
CALL DGTRST(JDUB,NRS,MMG,NSTR,NDES,JTRS,IPAS)
DO 23 IR1=1,NRS
IF(INDEQK.NE.0) THEN
JNDLR(IR1,JLR)=0
DO 24 IRG=1,NNRG
IF(NIRG(IRG).NE.0) THEN
CALL DRCGCF(IR1,IRG,IR1,DB,NI)
IF(NI.NE.0) JNDLR(IR1,JLR)=JNDLR(IR1,JLR)+NI
IF(IPAS(IR1).NE.0) THEN
IR2=IPAS(IR1)
CALL DRCGCF(IR1,IRG,IR2,DB,NI)
IF(NI.NE.0) JNDLR(IR1,JLR)=JNDLR(IR1,JLR)+NI
END IF
END IF
24 CONTINUE
ELSE
JNDLR(IR1,JLR)=99
END IF
23 CONTINUE
DO 33 J=1,JJMAX
IF(JLR.EQ.1) THEN
IF(JMAX(ILEFT).GE.J) J00LR(J,1)=JNDLR(IRNM(J,ILEFT),1)
IF(JMAX(ILEFT).LT.J) J00LR(J,1)=9
ELSE
IF(JMAX(IRIGHT).GE.J) J00LR(J,2)=JNDLR(IRNM(J,IRIGHT),2)
IF(JMAX(IRIGHT).LT.J) J00LR(J,2)=9
END IF
33 CONTINUE
C write(6,*) jjmax,jlr
C write(6,*) (j00lr(j,jlr),j=1,jjmax)
C if(jlr.eq.1) write(6,*) (irnm(j,ileft),j=1,jmax(ILEFT))
C if(jlr.eq.2) write(6,*) (irnm(j,iright),j=1,jmax(IRIGHT))
13 CONTINUE
RETURN
END
C SUBROUTINE DRCGCF ====*====3====*====4====*====5====*====6====*====7
C
C
C
C
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
SUBROUTINE DRCGCF(JR1,JRG,JR2,DB,NI)
IMPLICIT REAL*8 (A-H,O-Z)
COMPLEX*16 DB,CR,CR1,CR2,CRG,WA
COMPLEX*16 WWD1,WWD2,WWDG
COMMON/SPG2/IL,NG,IG(48),JV(2,3,48)
COMMON/SPG3/KB(3),ICB,MG,JG(48),JK(3,48)
& ,IZ,IFA(48,48),IFC,IDOUB,MTRG,JTRG(4,48),ITRC(3,48)
COMMON/SPG4/NR,NH,ND(12),IR(7,48,12),CR(48,12),IATR(12)
SAVE /SPG2/,/SPG3/,/SPG4/
DIMENSION KKB(3),KBB(3),INDA(6,3),INDB(6,3)
DIMENSION IR1(7,48),IR2(7,48),IRG(7,48),CR1(48),CR2(48),CRG(48)
DIMENSION WWD1(48,6),WWD2(48,6),WWDG(48,3),WA(6,3)
DIMENSION DB(6,3,6,5)
DIMENSION JGG(48),JGK(48)
DIMENSION NTAB(48)
DO 10 K=1,3
KBB(K)=KB(K)
10 CONTINUE
IBB=ICB
IDOUBB=IDOUB
KKB(1)=0
KKB(2)=0
KKB(3)=0
ICC=1
CALL TSIREP(KKB,ICC,0)
CALL CHKNIR(JRG,NR)
NHG=NH
NDG=ND(JRG)
DO 16 J=1,NHG
JGG(J)=JG(J)
CRG(J)=CR(J,JRG)
DO 17 K=1,7
IRG(K,J)=IR(K,J,JRG)
17 CONTINUE
16 CONTINUE
DO 11 K=1,3
KKB(K)=KBB(K)
11 CONTINUE
ICC=IBB
IDB=IDOUBB
CALL TSIREP(KKB,ICC,IDB)
CALL CHKNIR(JR1,NR)
NHK=NH
ND1=ND(JR1)
DO 12 J=1,NHK
JGK(J)=JG(J)
CR1(J)=CR(J,JR1)
DO 13 K=1,7
IR1(K,J)=IR(K,J,JR1)
13 CONTINUE
12 CONTINUE
CALL CHKNIR(JR2,NR)
ND2=ND(JR2)
DO 14 J=1,NHK
CR2(J)=CR(J,JR2)
DO 15 K=1,7
IR2(K,J)=IR(K,J,JR2)
15 CONTINUE
14 CONTINUE
C WRITE(6,*) NHK,NHG
C WRITE(6,*) (JGK(I),I=1,NHK)
C WRITE(6,*) (JGG(I),I=1,NHG)
CALL SUBGRP(NHK,JGK,NHG,JGG,NTAB,IND)
C WRITE(6,*) (NTAB(I),I=1,NHK)
C WRITE(6,*) IND
IF(IND.NE.0) GO TO 99
C
C CHARACTER TEST
C
CW=0
DO 23 I=1,NHK
CW=CW+DCONJG(CR1(I))*CR2(I)*CRG(NTAB(I))
C WRITE(6,*) I,NTAB(I),CW
C WRITE(6,600) CR1(I),CR2(I),CRG(NTAB(I))
C 600 FORMAT(3(' (',2F10.5,')'))
23 CONTINUE
CW=CW/NHK
NI=CW+0.5
C WRITE(6,*) ' NI=',NI
IF(NI.GT.5) THEN
WRITE(6,*) ' DIMENSION OF DB IS NOT ENUGH NI=',NI
WRITE(6,*) ' STOP AT 23 IN DIRSEL '
STOP
END IF
IF(NI.EQ.0) GO TO 80
DO 81 NA=1,ND2
DO 81 NBS=1,NDG
INDA(NA,NBS)=0
81 CONTINUE
C WRITE(6,*) ND1,ND2,NDG,NHK,NHG
DO 82 NNI=1,NI
C
C PROJECTION OPERATOR
C
C PIVOT IS SET AT THE COMPONENT OF SINGLE REPRESENTATIO
C WHICH IS NOT APPERED IN THE PREVIOUS PROJECTION FOR THE
C FIRST COMPONENT OF DUBLE REPRESENTATION
C
DO 83 IB=1,ND2
DO 831 IBS=1,NDG
IF(INDA(IB,IBS).EQ.1) GO TO 831
NB=IB
NBS=IBS
C---------------------------------------------------
C MATRIX ELEMENTS OF SINGLE REPRESENTATION
CALL GETMXE(IR2,ND2,NHK,NB,WWD2)
CALL GETMXE(IRG,NDG,NHG,NBS,WWDG)
C----------------------------------------------------
DO 84 N=1,ND1
C----------------------------------------------------
C SECOND SUFFIX OF PROJECTION OPERATOR IS FIXED AT 1
CALL GETMXE(IR1,ND1,NHK,1,WWD1)
C-----------------------------------------------------
DO 31 NA=1,ND2
DO 32 NAS=1,NDG
WA(NA,NAS)=0.0
INDB(NA,NAS)=0
32 CONTINUE
31 CONTINUE
DO 33 J=1,NHK
IF(ABS(WWD1(J,N)).GT.1.0D-7) THEN
DO 34 NA=1,ND2
IF(ABS(WWD2(J,NA)).LT.1.0D-7) GO TO 34
DO 35 NAS=1,NDG
IF(ABS(WWDG(NTAB(J),NAS)).LT.1.0D-7) GO TO 35
IF(N.EQ.1) INDB(NA,NAS)=1
WA(NA,NAS)=WA(NA,NAS)
& +DCONJG(WWD1(J,N))*WWD2(J,NA)*WWDG(NTAB(J),NAS)
35 CONTINUE
34 CONTINUE
END IF
33 CONTINUE
WC=0.0
DO 36 NA=1,ND2
DO 36 NAS=1,NDG
WC=WC+DCONJG(WA(NA,NAS))*WA(NA,NAS)
36 CONTINUE
C WRITE(6,*) 'NB=',NB,' NBS=',NBS,' N=',N,' NNI=',NNI,' WC=',WC
IF(WC.LT.1.0D-5) THEN
IF(N.EQ.1) GO TO 831
WRITE(6,*) ' THIS CAN NOT BE HAPPEN '
WRITE(6,*) ' STOP AT 36 IN DIRSEL '
STOP
END IF
DO 92 NA=1,ND2
DO 92 NAS=1,NDG
IF(INDB(NA,NAS).EQ.1) INDA(NA,NAS)=1
92 CONTINUE
WC=DSQRT(WC)
DO 37 NA=1,ND2
DO 37 NAS=1,NDG
DB(N,NAS,NA,NNI)=WA(NA,NAS)/WC
IF(ABS(DIMAG(DB(N,NAS,NA,NNI))).LT.1.0D-7) THEN
DB(N,NAS,NA,NNI)=DCMPLX(DBLE(DB(N,NAS,NA,NNI)),0.0D0)
END IF
IF(ABS(DBLE(DB(N,NAS,NA,NNI))).LT.1.0D-7) THEN
DB(N,NAS,NA,NNI)=DCMPLX(0.0D0,DIMAG(DB(N,NAS,NA,NNI)))
END IF
37 CONTINUE
84 CONTINUE
GO TO 82
831 CONTINUE
83 CONTINUE
WRITE(6,*) ' PIVOT SET ALGORISM IS NOT ENUGH '
WRITE(6,*) ' STOP AT 83 IN DIRSEL '
STOP
82 CONTINUE
80 CONTINUE
RETURN
99 WRITE(6,*) ' K-POINT GROUP IS NOT A SUBGROUP'
WRITE(6,*) ' STOP AT 99 IN DIRSEL'
STOP
END
SUBROUTINE S3N(X,Y,SM,XX,YY,N,NN,C1,CN,AMU1,ALMN)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION X(99),Y(99),SM(99),XX(200),YY(200)
DIMENSION H(99),ALM(99),AMU(99),C(99),P(99),Q(99),U(99)
N1=N-1
DO 110 I=2,N
H(I)=X(I)-X(I-1)
110 CONTINUE
DO 120 I=2,N1
ALM(I)=H(I+1)/(H(I)+H(I+1))
AMU(I)=1.0-ALM(I)
120 CONTINUE
DO 130 I=2,N1
C(I)=3.0*(ALM(I)*(Y(I)-Y(I-1))/H(I)+AMU(I)*(Y(I+1)-Y(I))/H(I+1))
130 CONTINUE
C(1)=C1
C(N)=CN
AMU(1)=AMU1
ALM(N)=ALMN
P(1)=2.0
Q(1)=-AMU(1)/P(1)
U(1)=C(1)/P(1)
DO 140 K=2,N
P(K)=ALM(K)*Q(K-1)+2.0
Q(K)=-AMU(K)/P(K)
U(K)=(C(K)-ALM(K)*U(K-1))/P(K)
140 CONTINUE
SM(N)=U(N)
DO 150 K=1,N1
K1=N1-K+1
SM(K1)=Q(K1)*SM(K1+1)+U(K1)
150 CONTINUE
DO 160 I=1,NN
XXI=XX(I)
DO 170 K=2,N
IF(XXI.GT.X(K)) GO TO 170
J1=K
GO TO 180
170 CONTINUE
180 J=J1-1
SMJ=SM(J)
SMJ1=SM(J1)
YJ=Y(J)
YJ1=Y(J1)
HJ1=H(J1)
XJ1=X(J1)-XXI
XJ=XXI-X(J)
HJ2=HJ1*HJ1
HJ3=HJ2*HJ1
YY(I)=SMJ*XJ1*XJ1*XJ/HJ2-SMJ1*XJ*XJ*XJ1/HJ2+YJ*XJ1*XJ1*(2.0*XJ+
& HJ1)/HJ3+YJ1*XJ*XJ*(2.0*XJ1+HJ1)/HJ3
160 CONTINUE
RETURN
END