IMPLICIT REAL*8 (A-H,O-Z) external func COMMON/DICUT/ICUT DIMENSION NRECPT(3,30),RECTAX(4,30) C REAL*8 XFU(3,2),WF(3) REAL*8 XX(3),XC(3),DR(3) real*8 func common/enrspl/EM(366145),PN(3),p(3),n(3),js INTEGER M(3),NN(3),NK(3) DIMENSION CO(3,2,100),RKG(4,20) DIMENSION JB(2,3) READ(1,*) READ(1,*) IL,NGEN,INV CALL TSPACE(IL) C CALL TSPHDS C CALL TSOPDS 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) CALL TSPGDS READ(1,*) A,B,C WRITE(6,601) A, B, C 601 FORMAT(' A=',F9.5,' B=',F9.5,' C=',F9.5) READ(1,*) CA,CB,CC CALL TSLATC(A,B,C,CA,CB,CC0) NB=6 DO 10 J=1,NB READ(31) IXMAX,NB,NNEE READ(31) (EM(K),K=1,NNEE) WRITE(6,*) IXMAX,NB,NNEE 10 CONTINUE PN(1)=1.0 PN(2)=1.0 PN(3)=1.0 n(1)=IXMAX-1 n(2)=n(1) n(3)=n(1) JS=3 DO 51 I=1,3 P(I)=Pn(I)/FLOAT(N(I)) 51 CONTINUE XFU(1,1)=0.0 XFU(2,1)=0.0 XFU(3,1)=-1.0D0 XFU(1,2)=2.0D0 XFU(2,2)=2.0D0 XFU(3,2)=1.0D0 SIZE=120.0 WF(1)=SIZE WF(2)=WF(1) WF(3)=WF(1) IEOH=1 AF=30.0 EF=0.41714 BT=0.0 GM=-20.0 ED=150.0 EL=1500.0 CALL AYPSTR(96) CALL AYORIG(30.0,30.0) CALL TPERSP(XFU,WF,AF,BT,GM,ED,EL) CALL TPTRAC(1) CALL TPSCHG(IEOH) CALL TPHILP(EF) CALL TPHILD(0.0,1,6,5,0) CALL LINEWD(1.6) CALL TPCIHL(3) CALL LINEWD(0.1) CALL TPCIHL(0) CALL TPHILD(0.0,0,0,0,0) DO 11 I=1,3 NN(I)=129 NK(I)=1 11 CONTINUE ICUT=0 CALL TPFXYZ(FUNC,EF,NN,NK,1) CALL LINEWD(1.3) ICUT=1 DR(1)=1.0D0 XC(1)=0.5D0 DR(2)=1.0D0 XC(2)=0.5D0 DR(3)=1.0D0 XC(3)=0.5D0 CALL TPSECT(FUNC,EF,DR,XC,0,1,257,257,4,4) CALL TPHILD(EF,1,4,4,0) CALL TPCIHL(2) CALL LINEWD(1.4) ICUT=2 DR(1)=1.0D0 XC(1)=0.5D0 DR(2)=1.0D0 XC(2)=0.5D0 DR(3)=-1.0D0 XC(3)=-0.5D0 CALL TPSECT(FUNC,EF,DR,XC,0,0,257,257,4,4) CALL LINEWD(1.5) ICUT=0 XC(1)=1.0 XC(2)=1.0 XC(3)=0.0 DR(1)=1.0 DR(2)=1.0 DR(3)=0.0D0 CALL TPSECT(FUNC,EF,DR,XC,0,0,129,129,4,4) CALL AYPEND STOP END FUNCTION FUNCD(XA) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION XA(3) WW=FUNC(XA) FUNCD=WW RETURN END FUNCTION LFUNC(XA) IMPLICIT REAL*8 (A-H,O-Z) COMMON/DICUT/ICUT DIMENSION XA(3) LFUNC=-1 IF(ICUT.EQ.0) THEN KX=1000000.0*XA(1) KY=1000000.0*XA(2) KZ=1000000.0*XA(3) CALL TSKFBZ(KX,KY,KZ,1000000,IND) IF(IND.NE.0.AND.XA(3).GT.0.0) LFUNC=1 ELSE IF(ICUT.EQ.1) THEN KX=1000005.0*XA(1) KY=1000005.0*XA(2) KZ=1000005.0*XA(3) CALL TSKFBZ(KX,KY,KZ,1000000,IND) IF(IND.NE.0.AND.XA(3).GT.0.0) LFUNC=1 IF(ABS(XA(2)-XA(1)).GT.1.0D0) LFUNC=1 IF(IND.EQ.0.AND.XA(3).LT.0.0) LFUNC=1 ELSE IF(ICUT.EQ.2) THEN IF(ABS(XA(2)-XA(1)).GT.1.0D0) LFUNC=1 IF(XA(3).GT.0.0) LFUNC=1 END IF RETURN END C FUNCTION TLINTP(XX,JS) C *** THREE DIMENSIONAL INTERPOLATION BY A LINEAR FORMULA C *** PERIODICITY AND F(X,Y,Z)=F(ABS(X),ABS(Y),ABS(Z)) C *** ARE ASSUMED. C *** IS=1 ORTHORHOMBIC SYMMETRY C *** IS=2 TETRAGONAL SYMMETRY C *** IS=3 CUBIC SYMMETRY REAL*8 function FUNC(XX) IMPLICIT REAL*8 (A-H,O-Z) REAL*8 KB C DIMENSION EM(NN),PM(3),N(3),XX(3),M(3),PN(3) common/enrspl/em(366145),pm(3),p(3),n(3),is DIMENSION IFAC(3) DIMENSION A(4),IPER(3,6),FF(2,2,2),F(8),KB(3) EQUIVALENCE(F(1),FF(1,1,1)) DATA IPER/1,2,3, 1,3,2, 3,1,2, 2,1,3, 2,3,1, 3,2,1/ DATA IFAC/1,2,4/ data niold,ipold/0,0/ save niold,ipold REAL*8 xx(3) C DO 53 I=1,3 C PM(I)=PN(I) C N(I)=M(I) C 53 CONTINUE C C *** KB IN THE FIRST ZONE *** C write(6,*) xx,ipold,niold DO 1 I=1,3 KB(I)=XX(I) 2 KB(I)=ABS(KB(I)) IF(KB(I).LE.PM(I)) GO TO 1 KB(I)=KB(I)-PM(I)*2.0 GO TO 2 1 CONTINUE IF(IS.EQ.1) GO TO 4 C *** KB(1).GE.KB(2).GE.KB(3) *** IF(KB(1).GE.KB(2)) GO TO 3 W=KB(1) KB(1)=KB(2) KB(2)=W 3 CONTINUE IF(IS.EQ.2) GO TO 4 IF(KB(2).GE.KB(3)) GO TO 4 W=KB(2) KB(2)=KB(3) KB(3)=W IF(KB(1).GE.KB(2)) GO TO 4 W=KB(1) KB(1)=KB(2) KB(2)=W C *** FIND THE BLOCK WHERE THE POINT IS *** 4 IX=KB(1)/P(1) IF(IX.EQ.N(1)) IX=N(1)-1 IY=KB(2)/P(2) IF(IY.EQ.N(2)) IY=N(2)-1 IZ=KB(3)/P(3) IF(IZ.EQ.N(3)) IZ=N(3)-1 IF(IS.EQ.1) NI=(IX*(N(2)+1)+IY)*(N(3)+1)+IZ+1 IF(IS.EQ.2) NI=((IX*(IX+1))/2+IY)*(N(3)+1)+IZ+1 IF(IS.EQ.3) NI=(IX*(IX+1)*(IX+2))/6 1 +(IY*(IY+1))/2+IZ+1 C *** IF THE BLOCK IS THE SAME, GO TO 40 *** IF(NI.EQ.NIOLD) GO TO 40 C *** PICK UP THE ENERGIES OF 8 POINTS *** NK=NI IF(IS.EQ.1) MX=(N(2)+1)*(N(3)+1) IF(IS.EQ.2) MX=(IX+1)*(N(3)+1) IF(IS.EQ.3) MX=((IX+1)*(IX+2))/2 DO 41 KX=1,2 DO 42 KY=1,2 DO 43 KZ=1,2 FF(KX,KY,KZ)=EM(NK) NK=NK+1 43 CONTINUE IF(IS.NE.3) NK=NK-2+(N(3)+1) IF(IS.EQ.3) NK=NK-2+IY+KY 42 CONTINUE IF(KX.LT.2) NK=NI+MX 41 CONTINUE 40 CONTINUE C *** FIND THE TETRAHEDRON *** X=KB(1)-P(1)*FLOAT(IX) Y=KB(2)-P(2)*FLOAT(IY) Z=KB(3)-P(3)*FLOAT(IZ) IF(X.GE.Y) GO TO 11 IF(Y.GE.Z) GO TO 12 IP=6 GO TO 13 12 IF(X.GE.Z) GO TO 24 IP=5 GO TO 13 24 IP=4 GO TO 13 11 IF(Y.GE.Z) GO TO 21 IF(X.GE.Z) GO TO 22 IP=3 GO TO 13 22 IP=2 GO TO 13 21 IP=1 C *** IF THE TETRAHEDRON IS THE SAME, GO TO 10 *** 13 IF(NIOLD.NE.NI) GO TO 14 IF(IPOLD.EQ.IP) GO TO 10 14 IPOLD=IP NIOLD=NI C *** DETERMINE THE COEFFICIENTS OF LINEAR FORM *** K1=IPER(1,IP) K2=IPER(2,IP) K3=IPER(3,IP) I1=IFAC(K1) I2=IFAC(K2) I3=IFAC(K3) A(4)=F(1) JA=I1+1 A(K1)=(F(JA)-F(1))/P(K1) JB=I1+I2+1 A(K2)=(F(JB)-F(JA))/P(K2) A(K3)=(F(8)-F(JB))/P(K3) 10 func=X*A(1)+Y*A(2)+Z*A(3)+A(4) RETURN END SUBROUTINE CLOCKM(ITIME) C C FOR UNIX C DIMENSION TERY(2) C C etime for AIX FORTRAN CALL ETIME_(TERY) C etime for SUN FORTARAN C CALL ETIME(TERY) ITIME=TERY(1)*1000 RETURN END SUBROUTINE PTIME(T) REAL*8 T CALL CLOCKM(IT) T=IT/3600000.0 RETURN END