Lanczos Diagonalization for Hubbard Ring
0010*#RUN: TIME=10 L=LIB/ASL7 0020C ------- * Dispersion Relation of the * ------- 0030C ------- * Extended Hubbard Hamiltonian * ------- 0040C ------- * * ------- 0050C ------- * H = - 1*( CiCj + CjCi ) * ------- 0060C ------- * + Ud* nn * ------- 0070C ------- * + Un( n+n )*( n+n ) * ------- 0080C ------- * * ------- 0090C ------- * 10 Site 364364 > 63504 * ------- 0100C ------- * 12 >>> 853776 * ------- 0110C ------- * 14 >>>>> 11778624 * ------- 0120C ------- * * ------- 0130C 0140 INTEGER NS, NK, NCF, CF(0:12870 ), DG(0:65536, 0:1) 0150 INTEGER NST(0:1), ST(0:12870,0:20,0:1), WG(0:12870, 0:1) 0160 INTEGER NE(0:1), SG(0:12870,0:20,0:1), S(0:12870,0:20,0:1) 0170 INTEGER CON(0:12870,0:20,0:1) 0180 INTEGER LD, LU(0:20), TD(0:99999, 0:2), TU(0:9999 ,0:20,0:2) 0190 REAL*8 UD, UN, MU, TIME1,TIME2, VEC(0:364364 ,0:3) 0200 REAL*8 BC(0:20), BS(0:20), NRM(0:364364), DIAG(0:364364) 0210C 0220 NS = 4 0230C NK = 0 0240 NE(0)= 2 0250 NE(1)= 2 0260C 0270 UD = 12.0D0 0280 UN = 6.0D0 0290 MU = 2.0D0*UN + UD/2.0D0 0300C 0310 WRITE(6,346) UD, UN, NE(0), NE(1) 0320 346 FORMAT(1H ,' UD= ',F10.7,' UN=',F10.7, ' Neu=',I3,' Ned=',I3 ) 0330C 0335C DO 90 I = 0 , NS/2 0340 DO 90 I = 0 , NS/2 0350 NK = I 0360C 0370 CALL CLOCK( TIME1 ) 0380 CALL CLOCK( TIME2 ) 0390 CALL ORTSET( NS, NK , BC, BS ) 0400 CALL CRESTA( NS, NE(0), NCF , CF ) 0410 CALL REDUCE( NS, NCF,CF, NST(0), ST(0,0,0) ) 0420 CALL CRESTA( NS, NE(1), NST(1), ST(0,0,1) ) 0430 CALL CREDG( NST(0), ST(0,0,0), DG(0,0) ) 0440 CALL CREDG( NST(1), ST(0,0,1), DG(0,1) ) 0450 CALL ROTATE( NS, NE , NST , ST, SG,WG, S, CON, DG ) 0460 CALL NORMAL( NS, BC , NST , WG, S, NRM ) 0470 CALL PRJOUT( NS, NST , CON,WG, S, NRM ) 0480 CALL DIAGON( NS, NE , NST , ST, UD,UN,MU, NRM, DIAG ) 0490 CALL TRANSF( NS, NE , NST , ST, SG,LD,LU,TD,TU, DG ) 0500 CALL LANTZ1( NS,BS,BC,NK,NST, S,CON, SG,LD,LU,TD,TU,DIAG,NRM,VEC) 0510C-WR- CALL NDIST( NS, NK , NST , ST, VEC ) 0520 CALL CDWORD( NS, NK , NST , ST, VEC ) 0530 90 CONTINUE 0540 CALL DTIME( TIME1, TIME2, ' IN ALL' ) 0550C 0560 STOP 0570 END 0580CCC 0590CC 0600C 0610 SUBROUTINE ORTSET( NS, NK, BC , BS ) 0620 INTEGER NS, NK 0630 REAL*8 PI, BC(0:20), BS(0:20) 0640C 0650 IF( MOD(NS,2) .NE. 0 .OR. NK .GT. NS/2 ) THEN 0660 WRITE(6,*)' Wrong Site Number ' 0670 STOP 0680 END IF 0690C 0700 PI = 4.0D0 * ATAN( 1.0D0 ) 0710 DO 10 I = 0, NS -1 0720 BC(I)= COS( 2.0D0*I*NK*PI/NS ) 0730 BS(I)= SIN( 2.0D0*I*NK*PI/NS ) 0740C-WR- WRITE(6,100) BC(I) , BS(I) 0750C 100 FORMAT(1H , 'COS= ', F10.5, ' SIN=' F10.5 ) 0760 10 CONTINUE 0770C 0780 RETURN 0790 END 0800CCC 0810CC 0820C 0830 SUBROUTINE CRESTA( NS, NE, NCONF, CONF ) 0840 INTEGER NS, NE, NCONF, CONF(0:12870) 0850C 0860 NCONF=0 0870 DO 20 I=0,2**NS-1 0880 IF( MOD(I ,2) +MOD(I/ 2** 1,2) +MOD(I/ 2** 2,2) 0890 & +MOD(I/ 2** 3,2) +MOD(I/ 2** 4,2) +MOD(I/ 2** 5,2) 0900 & +MOD(I/ 2** 6,2) +MOD(I/ 2** 7,2) +MOD(I/ 2** 8,2) 0910 & +MOD(I/ 2** 9,2) +MOD(I/ 2**10,2) +MOD(I/ 2**11,2) 0920 & +MOD(I/ 2**12,2) +MOD(I/ 2**13,2) +MOD(I/ 2**14,2) 0930 & +MOD(I/ 2**15,2) +MOD(I/ 2**16,2) +MOD(I/ 2**17,2) 0940 & +MOD(I/ 2**18,2) +MOD(I/ 2**19,2) .EQ. NE ) THEN 0950 CONF(NCONF)= I 0960 NCONF = NCONF+1 0970 END IF 0980 20 CONTINUE 0990C-WR- WRITE(6,100) NS,NE,NCONF 1000C 100 FORMAT(1H ,' NS=',I2,' NE=',I3,' NCONF=',I5) 1010C 1020 RETURN 1030 END 1040CCC 1050CC 1060C 1070 SUBROUTINE REDUCE( NS, NCF, CF , NST, ST ) 1080 INTEGER FLG, NS, NCF, CF(0:12870), NST, ST(0:12870) 1090C 1100 NST = 0 1110 DO 20 I = 0, NCF-1 1120 FLG = 0 1130 DO 10 J = 1, NS-1 1140 IF( MOD(CF(I)*2**J,2**NS) + CF(I)/2**(NS-J) .LT. CF(I) ) FLG = 1 1150 10 CONTINUE 1160 IF( FLG .EQ. 0 ) THEN 1170 ST(NST)= CF(I) 1180 NST = NST + 1 1190 END IF 1200 20 CONTINUE 1210C 1220C-WR- WRITE(6,100)NST 1230C 100 FORMAT(1H ,' Reduced Basis Set =', I5 ) 1240C 1250 RETURN 1260 END 1270CCC 1280CC 1290C 1300 SUBROUTINE CREDG( NST, ST , DG ) 1310 INTEGER NST, ST(0:12870), DG(0:65536) 1320C 1330 DO 10 I = 0, 65536 1340 DG(I)= -1 1350 10 CONTINUE 1360C 1370 DO 20 I = 0, NST - 1 1380 DG(ST(I))= I 1390 20 CONTINUE 1400C 1410 RETURN 1420 END 1430CCC 1440CC 1450C 1460 SUBROUTINE ROTATE( NS, NE, NST, ST, SG, WG, S, CON, DG ) 1470 INTEGER SP, NS, NE(0:1), NST(0:1), ST(0:12870,0:20,0:1) 1480 INTEGER FCT, S(0:12870,0:20,0:1), SG(0:12870,0:20,0:1) 1490 INTEGER CON(0:12870,0:20,0:1), WG(0:12870, 0:1) 1500 INTEGER DG(0:65536, 0:1) 1510C 1520 DO 90 SP = 0, 1 1530 FCT = (-1)**( NE(SP)-1 ) 1540C --- Initialization (m=0) --- 1550 DO 10 I = 0, NST(SP)-1 1560 SG(I,0,SP)= 1 1570 WG(I, SP)= 1 1580 S(I,0,SP)= 1 1590 CON(I,0,SP)= I 1600 10 CONTINUE 1610C 1620 DO 40 J = 1, NS -1 1630 DO 20 I = 0, NST(SP)-1 1640 ST(I,J,SP)= MOD(ST(I,J-1,SP)*2,2**NS) + ST(I,J-1,SP)/2**(NS-1) 1650 SG(I,J,SP)= SG(I,J-1,SP)*FCT** ( ST(I,J-1,SP)/2**(NS-1) ) 1660 S(I,J,SP)= 0 1670 20 CONTINUE 1680 DO 30 I = 0, NST(SP)-1 1690 CON(I,J,SP)= DG( ST(I,J,SP),SP ) 1700 IF( ST(I,J,SP) .EQ. ST(I,0,SP) ) THEN 1710 WG(I, SP) = WG(I, SP)+1 1720 S(I,J,SP) = SG(I,J,SP) 1730 END IF 1740 30 CONTINUE 1750 40 CONTINUE 1760C 1770 DO 50 I = 0, NST(SP)-1 1780C-WR- WRITE(6,100) WG(I,SP),(SG(I,J,SP),J=0,NS-1),(S(I,J,SP),J=0,NS-1) 1790C 100 FORMAT(1H ,'W= ',I2, 32I2 ) 1800 50 CONTINUE 1810C-WR- WRITE(6,*) 1820 90 CONTINUE 1830C 1840 RETURN 1850 END 1860CCC 1870CC 1880C 1890 SUBROUTINE NORMAL( NS, BC, NST, WG, S, NRM ) 1900 INTEGER NST(0:1), NS, WG(0:12870,0:1), S(0:12870,0:20,0:1) 1910 REAL*8 TP, BC(0:20), NRM(0:364364) 1920C 1930 DO 90 I = 0, NST(0)-1 1940 IF( WG(I,0) .EQ. 1 ) THEN 1950 DO 10 J = 0, NST(1)-1 1960 NRM( I*NST(1)+J ) = 1.0D0 1970 10 CONTINUE 1980 ELSE 1990 DO 40 J = 0, NST(1) -1 2000 TP = 0.0D0 2010 DO 30 K = 0, NS -1 2020 TP = TP + BC(K)*S(I,K,0)*S(J,K,1) 2030 30 CONTINUE 2040 IF( TP .GT. 1.0E-8 ) THEN 2050 NRM( I*NST(1)+J ) = 1.0D0 / SQRT(TP) 2060 ELSE 2070 NRM( I*NST(1)+J ) = 0.0D0 2080 END IF 2090 40 CONTINUE 2100 END IF 2110 90 CONTINUE 2120C 2130 RETURN 2140 END 2150CCC 2160CC 2170C 2180 SUBROUTINE PRJOUT( NS, NST, CON, WG, S, NRM ) 2190 INTEGER NS, NST(0:1), CON(0:12870,0:20,0:1), WG(0:12870,0:1) 2200 INTEGER S(0:12870,0:20,0:1) 2210 REAL*8 NRM(0:364364 ) 2220C 2230 DO 90 I = 0, NST(0)-1 2240 IF( WG(I,0) .NE. 1 ) THEN 2250 DO 80 J = 0, NST(1)-1 2260 IF( NRM( I*NST(1)+J ) .GT. 1.0E-8 ) THEN 2270 DO 70 K = 1, NS-1 2280 IF( S(I,K,0) .NE. 0 .AND. CON(J,K,1) .NE. J ) THEN 2290 NRM( I*NST(1) + CON(J,K,1) )= 0.0D0 2300 END IF 2310 70 CONTINUE 2320 END IF 2330 80 CONTINUE 2340 END IF 2350 90 CONTINUE 2360C 2370C-WR- WRITE(6,100) ( NRM(I), I=0, NST(0)*NST(1)-1 ) 2380C 100 FORMAT(1H ,'Normalization= ', F10.6 ) 2390C 2400 RETURN 2410 END 2420CCC 2430CC 2440C 2450 SUBROUTINE DTIME(TIME1,TIME2,ASCII) 2460 REAL*8 TIME1,TIME2 2470 CHARACTER*6 ASCII 2480 IF(TIME1.GE.TIME2) THEN 2490 CALL CLOCK(TIME2) 2500 WRITE(6,100) TIME2-TIME1,ASCII 2510 100 FORMAT(1H ,F15.8,' SEC WAS PASSED IN ',A10) 2520 ELSE 2530 CALL CLOCK(TIME1) 2540 WRITE(6,100) TIME1-TIME2,ASCII 2550 END IF 2560 RETURN 2570 END 2580CCC 2590CC 2600C 2610 SUBROUTINE DIAGON( NS, NE, NST, ST, UD, UN, MU, NRM, DIAG ) 2620 INTEGER NS, NE(0:1), NST(0:1), ST(0:12870,0:20,0:1) 2630 INTEGER T1, T2, T3, CNT, IS, ID 2640 REAL*8 UD, UN, MU, NRM(0:364364), DIAG(0:364364) 2650C 2660 DO 20 I = 0, NST(0)-1 2670 DO 10 J = 0, NST(1)-1 2680 T1 = ST(I,0,0) 2690 T2 = ST(J,0,1) 2700 CNT = I*NST(1)+J 2710 DIAG(CNT)= MU*( -NE(0) -NE(1) ) 2720C 2730 T3= IAND(T1,T2) 2740 DIAG(CNT)= DIAG(CNT) + UD*( 2750 & + MOD( T3 ,2)+ MOD( T3 / 2** 1,2) 2760 & + MOD( T3 / 2** 2,2)+ MOD( T3 / 2** 3,2) 2770 & + MOD( T3 / 2** 4,2)+ MOD( T3 / 2** 5,2) 2780 & + MOD( T3 / 2** 6,2)+ MOD( T3 / 2** 7,2) 2790 & + MOD( T3 / 2** 8,2)+ MOD( T3 / 2** 9,2) 2800 & + MOD( T3 / 2**10,2)+ MOD( T3 / 2**11,2) 2810 & + MOD( T3 / 2**12,2)+ MOD( T3 / 2**13,2) 2820 & + MOD( T3 / 2**14,2)+ MOD( T3 / 2**15,2) 2830 & + MOD( T3 / 2**16,2)+ MOD( T3 / 2**17,2) 2840 & + MOD( T3 / 2**18,2)+ MOD( T3 / 2**19,2) ) 2850C 2860 10 CONTINUE 2870 20 CONTINUE 2880C 2890 DO 50 IS = 0, NS-1 2900 IF( IS .NE. NS-1 ) THEN 2910 ID = IS + 1 2920 ELSE 2930 ID = 0 2940 END IF 2950C 2960 DO 40 I = 0, NST(0)-1 2970 DO 30 J = 0, NST(1)-1 2980 T1 = ST(I,0,0) 2990 T2 = ST(J,0,1) 3000 CNT = I*NST(1)+J 3010 DIAG(CNT)= DIAG(CNT) 3020 & + UN*( MOD(T1/2**IS,2) + MOD(T2/2**IS,2) ) 3030 & *( MOD(T1/2**ID,2) + MOD(T2/2**ID,2) ) 3040C 3050 30 CONTINUE 3060 40 CONTINUE 3070 50 CONTINUE 3080C 3090 DO 60 I = 0, NST(0)*NST(1)-1 3100 IF( NRM(I) .LT. 1.0E-8 ) DIAG(I)= 0.0D0 3110 60 CONTINUE 3120C 3130C-WR- WRITE(6,100) ( DIAG(I),I=0,NST(0)*NST(1)-1 ) 3140C 100 FORMAT(1H ,' Diagonal = ',F10.6 ) 3150C 3160 RETURN 3170 END 3180CCC 3190CC 3200C 3210 SUBROUTINE TRANSF(NS, NE, NST, ST, SG, LD, LU, TD, TU, DG ) 3220 INTEGER ITP, NS, NE(0:1), NST(0:1),LD, LU(0:20) 3230 INTEGER ST(0:12870,0:20,0:1), SG(0:12870,0:20,0:1) 3240 INTEGER TD(0:99999,0:2) , TU(0: 9999,0:20,0:2) 3250 INTEGER DG(0:65536,0:1), FL , MSK 3260C ---- Down Spin Part ---- 3270 LD = 0 3280 DO 20 I = 0, NS-1 3290 IF( I .NE. NS-1 ) THEN 3300 FL = 1 3310 MSK = 3* 2**I 3320 ELSE 3330 FL = (-1)**( NE(1)-1 ) 3340 MSK = 1+ 2**( NS -1 ) 3350 END IF 3360 DO 10 J = 0, NST(1)-1 3370 ITP = IEOR( MSK, ST(J,0,1) ) 3380 IF( DG(ITP,1) .NE. -1 ) THEN 3390 TD(LD ,0) = J 3400 TD(LD ,1) = DG(ITP,1) 3410 TD(LD ,2) =-FL 3420 LD = LD + 1 3430 END IF 3440 10 CONTINUE 3450 20 CONTINUE 3460C ---- Up Spin Part ---- 3470 DO 50 M = 0, NS-1 3480 LU(M)= 0 3490 DO 40 I = 0, NS-1 3500 IF( I .NE. NS-1 ) THEN 3510 FL = 1 3520 MSK = 3* 2**I 3530 ELSE 3540 FL = (-1)**( NE(0)-1 ) 3550 MSK = 1+ 2**( NS -1 ) 3560 END IF 3570 DO 30 J = 0, NST(0)-1 3580 ITP = IEOR( MSK, ST(J,M,0) ) 3590 IF( DG(ITP,0) .NE. -1 ) THEN 3600 TU(LU(M),M,0) = J 3610 TU(LU(M),M,1) = DG(ITP,0) 3620 TU(LU(M),M,2) =-FL*SG(J,M,0) 3630 LU(M) = LU(M) + 1 3640 END IF 3650 30 CONTINUE 3660 40 CONTINUE 3670 50 CONTINUE 3680C 3690C-WR- WRITE(6,100) LD 3700C-WR- WRITE(6,110)(LU(I),I=0,NS-1) 3710C 100 FORMAT(1H ,' Transfer Down=', I6 ) 3720C 110 FORMAT(1H ,' Transfer Up =',10I6 ) 3730 RETURN 3740 END 3750CCC 3760CC 3770C 3780 SUBROUTINE INPRO( NVEC, IDIS, ISRC, PROD, VEC ) 3790 INTEGER NVEC, IDIS, ISRC 3800 REAL*8 PROD, VEC(0:364364,0:3) 3810 PROD=0.0D0 3820 DO 10 I=0, NVEC-1 3830 PROD=PROD+VEC(I,IDIS)*VEC(I,ISRC) 3840 10 CONTINUE 3850 RETURN 3860 END 3870CCC 3880CC 3890C 3900 SUBROUTINE ADDER( NVEC, IDIS, ISRC1, ISRC2, A, B, VEC ) 3910 INTEGER NVEC, IDIS, ISRC1, ISRC2 3920 REAL*8 VEC(0:364364,0:3), A, B 3930 DO 10 I=0, NVEC-1 3940 VEC(I,IDIS)=A*VEC(I,ISRC1)+B*VEC(I,ISRC2) 3950 10 CONTINUE 3960 RETURN 3970 END 3980CCC 3990CC 4000C 4010 SUBROUTINE DIVER( NVEC, IDIS, ISRC, D, VEC) 4020 INTEGER NVEC, IDIS, ISRC 4030 REAL*8 D, VEC(0:364364,0:3) 4040 DO 10 I=0, NVEC-1 4050 VEC(I,IDIS)=VEC(I,ISRC)/D 4060 10 CONTINUE 4070 RETURN 4080 END 4090CCC 4100CC 4110C 4120 SUBROUTINE VEC0(NVEC,VEC) 4130 INTEGER NVEC 4140 REAL*8 VEC(0:364364,0:3) 4150 DO 10 I=0,3 4160 DO 5 J=0,NVEC -1 4170 VEC(J,I)=0.0D0 4180 5 CONTINUE 4190 10 CONTINUE 4200 RETURN 4210 END 4220CCC 4230CC 4240C 4250 SUBROUTINE VMOVE( NVEC, IDIS, ISRC, VEC) 4260 INTEGER NVEC, IDIS, ISRC 4270 REAL*8 VEC(0:364364,0:3) 4280 DO 10 I=0, NVEC-1 4290 VEC(I,IDIS)=VEC(I,ISRC) 4300 10 CONTINUE 4310 RETURN 4320 END 4330CCC 4340CC 4350C 4360 SUBROUTINE VECINI( NS, NK, NST, NRM, VEC ) 4370 INTEGER NS, NK, NST(0:1), NV 4380 REAL*8 PROD, VEC(0:364364,0:3), NRM(0:364364) 4390C 4400 NV = NST(0)*NST(1) 4410 IF( NK .EQ. 0 .OR. NK .EQ. NS/2 ) THEN 4420 CALL VEC0(NV,VEC) 4430 DO 10 I = 0, NV - 1 4440 VEC(I,0) = (MOD(I,3)+MOD(I,5)+1)/SQRT(DBLE(NV))*NRM(I) 4450 10 CONTINUE 4460 ELSE 4470 CALL VEC0(NV*2,VEC) 4480 DO 20 I = 0, NV-1 4490 VEC( I,0) = (MOD(I,3)+MOD(I,5)+1)/SQRT(DBLE(NV))*NRM(I) 4500 VEC(NV+I,0) = (MOD(I,7)+MOD(I,5)+1)/SQRT(DBLE(NV))*NRM(I) 4510 20 CONTINUE 4520 NV = NV * 2 4530 END IF 4540 CALL INPRO(NV,0,0,PROD,VEC) 4550 PROD=SQRT(PROD) 4560 CALL DIVER(NV,0,0,PROD,VEC) 4570 RETURN 4580 END 4590CCC 4600CC 4610C 4620 SUBROUTINE MULT( ID, IS, NS, NK, NST, LU, LD, TU, TD, S, SG, CON, 4630 & BS, BC,NRM, DIAG, VEC ) 4640 INTEGER ID, IS, NS, NK, NST(0:1),LU(0:20), LD 4650 INTEGER S(0:12870,0:20,0:1), CON(0:12870,0:20,0:1) 4660 INTEGER SG(0:12870,0:20,0:1), TU(0:9999,0:20,0:2),TD(0:99999,0:2) 4670 REAL*8 BS( 0:20) , BC( 0:20) 4680 REAL*8 VEC(0:364364,0:3) ,NRM(0:364364) , DIAG(0:364364) 4690C 4700 CALL SCALE( ID, IS, NS, NK, NST, NRM, VEC ) 4710 CALL MULTUP( ID, IS, NS, NK, NST, LU, TU, CON, SG, BS,BC, VEC ) 4720 CALL MULTDW( ID, IS, NS, NK, NST, LD, TD, S, CON, SG, BS,BC, VEC ) 4730 CALL RESCAL( ID, IS, NS, NK, NST, NRM, VEC ) 4740 CALL MULDIA( ID, IS, NS, NK, NST, DIAG, VEC ) 4750C WRITE(6,*) 4760C WRITE(6,100)(VEC(I,IS),I=0,2*NST(0)*NST(1)-1) 4770C WRITE(6,100)(VEC(I,ID),I=0,2*NST(0)*NST(1)-1) 4780C 100 FORMAT( 1H ,8F9.3 ) 4790C 4800 RETURN 4810 END 4820CCC 4830CC 4840C 4850 SUBROUTINE SCALE( ID, IS, NS, NK, NST, NRM, VEC ) 4860 INTEGER NV, ID, IS, NS, NK, NST(0:1) 4870 REAL*8 NRM(0:364364) , VEC(0:364364, 0:3) 4880C 4890 NV = NST(0)*NST(1) 4900C 4910 IF( NK .EQ. 0 .OR. NK .EQ. NS/2 ) THEN 4920 DO 10 I = 0, NV - 1 4930 VEC(I,ID)= 0.0D0 4940 VEC(I,IS)= VEC(I,IS) * NRM(I) 4950 10 CONTINUE 4960 ELSE 4970 DO 20 I = 0, NV - 1 4980 VEC(I ,ID)= 0.0D0 4990 VEC(I+NV,ID)= 0.0D0 5000 VEC(I ,IS)= VEC(I ,IS) * NRM(I) 5010 VEC(I+NV,IS)= VEC(I+NV,IS) * NRM(I) 5020 20 CONTINUE 5030 END IF 5040C 5050 RETURN 5060 END 5070CCC 5080CC 5090C 5100 SUBROUTINE RESCAL( ID, IS, NS, NK, NST, NRM, VEC ) 5110 INTEGER NV, ID, IS, NS, NK, NST(0:1) 5120 REAL*8 NRM(0:364364) , VEC(0:364364, 0:3) 5130C 5140 NV = NST(0)*NST(1) 5150C 5160 IF( NK .EQ. 0 .OR. NK .EQ. NS/2 ) THEN 5170 DO 10 I = 0, NV -1 5180 VEC(I,ID)= VEC(I,ID) * NRM(I) 5190 IF( NRM(I) .GT. 1.0E-8 ) THEN 5200 VEC(I,IS)= VEC(I,IS) / NRM(I) 5210 END IF 5220 10 CONTINUE 5230 ELSE 5240 DO 20 I = 0, NV -1 5250 VEC(I ,ID)= VEC(I ,ID) * NRM(I) 5260 VEC(I+NV,ID)= VEC(I+NV,ID) * NRM(I) 5270 IF( NRM(I) .GT. 1.0E-8 ) THEN 5280 VEC(I ,IS)= VEC(I ,IS) / NRM(I) 5290 VEC(I+NV,IS)= VEC(I+NV,IS) / NRM(I) 5300 END IF 5310 20 CONTINUE 5320 END IF 5330C 5340 RETURN 5350 END 5360CCC 5370CC 5380C 5390 SUBROUTINE MULDIA( ID, IS, NS, NK, NST, DIAG, VEC ) 5400 INTEGER NV, ID, IS, NS, NK, NST(0:1) 5410 REAL*8 DIAG(0:364364), VEC(0:364364,0:3) 5420C 5430 NV = NST(0) * NST(1) 5440C 5450 IF( NK .EQ. 0 .OR. NK .EQ. NS/2 ) THEN 5460 DO 10 I = 0, NV -1 5470 VEC(I,ID)= VEC(I,ID) + DIAG(I)*VEC(I,IS) 5480 10 CONTINUE 5490 ELSE 5500 DO 20 I = 0, NV -1 5510 VEC(I ,ID)= VEC(I ,ID) + DIAG(I)*VEC(I ,IS) 5520 VEC(I+NV,ID)= VEC(I+NV,ID) + DIAG(I)*VEC(I+NV,IS) 5530 20 CONTINUE 5540 END IF 5550C 5560 RETURN 5570 END 5580CCC 5590CC 5600C 5610 SUBROUTINE MULTUP(ID,IS, NS,NK, NST, LU ,TU, CON,SG,BS,BC, VEC) 5620 INTEGER BX,BY, NV,ID,IS, NS,NK, NST(0:1) ,TU(0:9999 , 0:20, 0:2) 5630 INTEGER CON(0:12870,0:20,0:1), SG(0:12870,0:20,0:1) , LU(0:20) 5640 REAL*8 VEC(0:364364, 0:3), BS(0:20) , BC(0:20) , FCT 5650C 5660 NV = NST(0)*NST(1) 5670 IF( NK .EQ. 0 .OR. NK .EQ. NS/2 ) THEN 5680 DO 30 M = 0, NS -1 5690 DO 20 I = 0, LU(M) -1 5700 BX = TU(I,M,0) * NST(1) 5710 BY = TU(I,M,1) * NST(1) 5720 FCT = TU(I,M,2) * BC(M) 5730*VDIR NODEP 5740 DO 10 L = 0, NST(1)-1 5750 VEC( BY+CON(L,M,1), ID ) = VEC( BY+CON(L,M,1), ID ) 5760 & + SG(L,M,1) * FCT * VEC( BX+L , IS ) 5770 10 CONTINUE 5780 20 CONTINUE 5790 30 CONTINUE 5800 ELSE 5810 DO 80 M = 0, NS -1 5820 IF( ABS( BC(M) ) .GT. 1.0E-8 ) THEN 5830 DO 50 I = 0, LU(M) -1 5840 BX = TU(I,M,0) * NST(1) 5850 BY = TU(I,M,1) * NST(1) 5860 FCT = TU(I,M,2) * BC(M) 5870*VDIR NODEP 5880 DO 40 L = 0, NST(1)-1 5890 VEC( BY+CON(L,M,1),ID ) = VEC( BY+CON(L,M,1),ID) 5900 & + SG(L,M,1) * FCT * VEC( BX+L ,IS) 5910 VEC(NV+BY+CON(L,M,1),ID ) = VEC(NV+BY+CON(L,M,1),ID) 5920 & + SG(L,M,1) * FCT * VEC(NV+BX+L ,IS) 5930 40 CONTINUE 5940 50 CONTINUE 5950 END IF 5960 IF( ABS( BS(M) ) .GT. 1.0E-8 ) THEN 5970 DO 70 I = 0, LU(M) -1 5980 BX = TU(I,M,0) * NST(1) 5990 BY = TU(I,M,1) * NST(1) 6000 FCT = TU(I,M,2) * BS(M) 6010*VDIR NODEP 6020 DO 60 L = 0, NST(1)-1 6030C 6040 VEC( BY+CON(L,M,1),ID ) = VEC( BY+CON(L,M,1),ID) 6050 & + SG(L,M,1) * FCT * VEC(NV+BX+L ,IS) 6060 VEC(NV+BY+CON(L,M,1),ID ) = VEC(NV+BY+CON(L,M,1),ID) 6070 & - SG(L,M,1) * FCT * VEC( BX+L ,IS) 6080 60 CONTINUE 6090 70 CONTINUE 6100 END IF 6110 80 CONTINUE 6120 END IF 6130C 6140 RETURN 6150 END 6160CCC 6170CC 6180C 6190 SUBROUTINE MULTDW(ID,IS, NS,NK, NST, LD, TD, S, CON,SG,BS,BC, VEC) 6200 INTEGER BX,BY, NV,ID,IS, NS,NK, NST(0:1),LD,TD(0:99999, 0:2) 6210 INTEGER CON(0:12870,0:20,0:1), SG(0:12870,0:20,0:1) 6220 INTEGER S(0:12870,0:20,0:1) 6230 REAL*8 VEC(0:364364, 0:3), BS(0:20) , BC(0:20) , FCT , TP 6240C 6250 NV = NST(0)*NST(1) 6260 IF( NK .EQ. 0 .OR. NK .EQ. NS/2 ) THEN 6270 DO 20 I = 0, LD-1 6280 BX = TD(I,0) 6290 BY = TD(I,1) 6300 TP = TD(I,2) 6310*VDIR NODEP 6320 DO 10 J = 0, NST(0)-1 6330 VEC(J*NST(1)+BY,ID) = VEC(J*NST(1)+BY,ID) 6340 & + VEC(J*NST(1)+BX,IS) * TP 6350 10 CONTINUE 6360 20 CONTINUE 6370C 6380 DO 50 M = 1, NS-1 6390 DO 40 I = 0, NST(0)-1 6400 IF( S(I,M,0) .NE. 0 ) THEN 6410 FCT = S(I,M,0) * BC(M) 6420*VDIR NOVECTOR 6430 DO 30 J = 0, LD-1 6440 BX = CON( TD(J,0), NS-M, 1 ) 6450 BY = TD(J,1) 6460 TP = TD(J,2)* FCT* SG(TD(J,0), NS-M, 1 ) 6470 VEC(I*NST(1)+BY,ID) = VEC(I*NST(1)+BY,ID) 6480 & + VEC(I*NST(1)+BX,IS) * TP 6490 30 CONTINUE 6500 END IF 6510 40 CONTINUE 6520 50 CONTINUE 6530 ELSE 6540 DO 70 I = 0, LD-1 6550 BX = TD(I,0) 6560 BY = TD(I,1) 6570 TP = TD(I,2) 6580*VDIR NODEP 6590 DO 60 J = 0, NST(0)-1 6600 VEC( J*NST(1)+BY,ID) = VEC( J*NST(1)+BY,ID) 6610 & + VEC( J*NST(1)+BX,IS) * TP 6620 VEC(NV+J*NST(1)+BY,ID) = VEC(NV+J*NST(1)+BY,ID) 6630 & + VEC(NV+J*NST(1)+BX,IS) * TP 6640 60 CONTINUE 6650 70 CONTINUE 6660C 6670 DO 120 M = 1, NS-1 6680 IF( ABS( BC(M) ) .GT. 1.0E-8 ) THEN 6690 DO 90 I = 0, NST(0)-1 6700 IF( S(I,M,0) .NE. 0 ) THEN 6710 FCT = BC(M) * S(I,M,0) 6720*VDIR NOVECTOR 6730 DO 80 J = 0, LD-1 6740 BX = CON( TD(J,0), NS-M, 1 ) 6750 BY = TD(J,1) 6760 TP = TD(J,2)* FCT *SG(TD(J,0),NS-M,1) 6770 VEC( I*NST(1)+BY,ID) = VEC( I*NST(1)+BY,ID) 6780 & + TP *VEC( I*NST(1)+BX,IS) 6790 VEC(NV+I*NST(1)+BY,ID) = VEC(NV+I*NST(1)+BY,ID) 6800 & + TP *VEC(NV+I*NST(1)+BX,IS) 6810 80 CONTINUE 6820 END IF 6830 90 CONTINUE 6840 END IF 6850 IF( ABS( BS(M) ) .GT. 1.0E-8 ) THEN 6860 DO 110 I = 0, NST(0)-1 6870 IF( S(I,M,0) .NE. 0 ) THEN 6880 FCT = BS(M) * S(I,M,0) 6890*VDIR NOVECTOR 6900 DO 100 J = 0, LD-1 6910 BX = CON( TD(J,0), NS-M, 1 ) 6920 BY = TD(J,1) 6930 TP = TD(J,2)* FCT *SG(TD(J,0),NS-M,1) 6940 VEC( I*NST(1)+BY,ID) = VEC( I*NST(1)+BY,ID) 6950 & + TP *VEC(NV+I*NST(1)+BX,IS) 6960 VEC(NV+I*NST(1)+BY,ID) = VEC(NV+I*NST(1)+BY,ID) 6970 & - TP *VEC( I*NST(1)+BX,IS) 6980 100 CONTINUE 6990 END IF 7000 110 CONTINUE 7010 END IF 7020 120 CONTINUE 7030 END IF 7040C 7050 RETURN 7060 END 7070CCC 7080CC 7090C 7100 SUBROUTINE LANTZ1(NS,BS,BC,NK,NST,S,CON,SG,LD,LU, TD, TU, 7110 & DIAG,NRM,VEC) 7120 INTEGER NVEC,LK,LLANTZ,NS,NK,NST(0:1), LD,LU(0:20) 7130 INTEGER S(0:12870,0:20,0:1), SG(0:12870,0:20,0:1) 7140 INTEGER CON(0:12870,0:20,0:1), TU(0:9999 ,0:20,0:2) 7150 INTEGER TD(0:99999, 0:2) 7160 REAL*8 DIAG(0:364364), NRM(0:364364), VEC(0:364364,0:3) 7170 REAL*8 MEER(1000,10),TEER(1000),ENERGY,ALPHA(1000),BEETA(1000) 7180 REAL*8 BS(0:20) , BC(0:20),E2,E2P, ESP(1000) 7190C 7200 E2 = 0.0D0 7210C LLANTZ=20 7220 IF( NK .EQ. 0 .OR. NK .EQ. NS/2 ) THEN 7230 NVEC=NST(0)*NST(1) 7240 ELSE 7250 NVEC=NST(0)*NST(1)*2 7260 END IF 7270C 7280 IF( NVEC .GT. 364364 ) THEN 7290 WRITE(6,*) 'Too Large Matrix. Please Increase VEC Size ' 7300 STOP 7310 END IF 7320C 7330 CALL VEC0( NVEC, VEC ) 7340 CALL VECINI( NS,NK, NST, NRM, VEC ) 7350 CALL VMOVE( NVEC, 3 , 0 , VEC ) 7360C (I) 7370 CALL MULT(1,0,NS,NK,NST,LU,LD,TU,TD,S,SG,CON,BS,BC,NRM,DIAG,VEC) 7380C (II) 7390C ******** LK IS AN INDEX FOR ALPHA AND BEETA ******** 7400 LK=1 7410 LLANTZ=0 7420 11 CONTINUE 7430C ******** REPORT PER 5 ROTATION ******** 7440 E2P=E2 7450 DO 40 J=0,4 7460 CALL MACRO1(NS,NK,NST,LU,LD,TU,TD,S,SG,CON,BS,BC,NRM,DIAG,VEC, 7470 & TEER, ALPHA, BEETA, LK, 0 ) 7480 40 CONTINUE 7490 CALL MALANZ(LK,ALPHA,BEETA,ENERGY,E2,MEER,ESP) 7500 LLANTZ=LLANTZ+1 7510 IF(LLANTZ.LE.40.AND.ABS(E2P-E2).GE.1.0E-14) GOTO 11 7520 IF(LLANTZ.LE.10 ) GOTO 11 7530 LLANTZ=LLANTZ-1 7540 WRITE(6,123) ( ESP(IK),IK=1,10 ) 7550C-WR- WRITE(8,124) ( ESP(IK),IK=1,21 ) 7560 123 FORMAT( 1H ,10F7.3 ) 7570C 124 FORMAT( 1H ,7F10.5 ) 7580C 7590 DO 88 I=1,1000 7600 TEER(I)=MEER(I,1) 7610 88 CONTINUE 7620 DO 99 I=0,2 7630 DO 98 J=0,NVEC -1 7640 VEC(J,I)=0.0D0 7650 98 CONTINUE 7660 99 CONTINUE 7670 CALL VMOVE(NVEC, 0 , 3 , VEC) 7680 CALL ADDER(NVEC, 3 , 3 , 0 ,TEER(1), 0.0D0,VEC) 7690 CALL MULT(1,0,NS,NK,NST,LU,LD,TU,TD,S,SG,CON,BS,BC,NRM,DIAG,VEC) 7700 LK=1 7710 DO 150 I=0,LLANTZ 7720 DO 140 J=0,4 7730 CALL MACRO1(NS,NK,NST,LU,LD,TU,TD,S,SG,CON,BS,BC,NRM,DIAG,VEC, 7740 & TEER, ALPHA, BEETA, LK, 1 ) 7750 140 CONTINUE 7760 150 CONTINUE 7770 CALL INVERS(NS,NK,NST,LU,LD,TU,TD,S,SG,CON,BS,BC, 7780 & NRM,DIAG,ENERGY,VEC) 7790C-WR- WRITE(6,179) ENERGY 7800C 179 FORMAT(1H ,F20.15) 7810 RETURN 7820 END 7830CCC 7840CC 7850C 7860 SUBROUTINE MACRO1(NS,NK,NST,LU,LD,TU,TD,S,SG,CON,BS,BC,NRM, 7870 & DIAG, VEC , TEER, ALPHA, BEETA, LK, ISW ) 7880 INTEGER NVEC,LK,ISW,NS,NK,NST(0:1),LD,LU(0:20) 7890 INTEGER S(0:12870,0:20,0:1), SG(0:12870,0:20,0:1) 7900 INTEGER CON(0:12870,0:20,0:1), TU(0:9999 ,0:20,0:2) 7910 INTEGER TD(0:99999, 0:2) 7920 REAL*8 DIAG(0:364364),NRM(0:364364),VEC(0:364364,0:3) 7930 REAL*8 TEER(1000),ALPHA(1000),BEETA(1000),PROD,BS(0:20),BC(0:20) 7940C 7950 IF( NK .EQ. 0 .OR. NK .EQ. NS/2 ) THEN 7960 NVEC=NST(0)*NST(1) 7970 ELSE 7980 NVEC=NST(0)*NST(1)*2 7990 END IF 8000C 8010 CALL INPRO( NVEC, 0 , 1 , PROD,VEC) 8020 ALPHA(LK)=PROD 8030 PROD=-1.0*PROD 8040 CALL ADDER( NVEC, 2 , 1 , 0 ,1.0D0,PROD,VEC) 8050 CALL INPRO( NVEC, 2 , 2 , PROD,VEC) 8060 PROD=SQRT(PROD) 8070 BEETA(LK)=PROD 8080 CALL DIVER( NVEC, 2 , 2 , PROD,VEC) 8090 CALL MULT(1,2,NS,NK,NST,LU,LD,TU,TD,S,SG,CON,BS,BC,NRM,DIAG,VEC) 8100 PROD=-1.0*PROD 8110 CALL ADDER( NVEC, 1 , 1 , 0 ,1.0D0,PROD,VEC) 8120 CALL VMOVE( NVEC, 0 , 2 , VEC) 8130 LK=LK+1 8140 IF(ISW.EQ.1) CALL ADDER( NVEC,3,3,0,1.0D0,TEER(LK),VEC) 8150 RETURN 8160 END 8170CCC 8180CC 8190C 8200 SUBROUTINE MALANZ(LK,ALPHA,BEETA,ENERGY,E2,VE,E) 8210 INTEGER LK,N,M,LNV,ISW,IW1(1000),IERR 8220 REAL*8 ALPHA(1000),BEETA(1000),ENERGY,E2 8230 REAL*8 EPS,E(1000),VE(1000,10),W1(6000) 8240 N = LK-1 8250 EPS = -1.0 8260 M = MIN(N,10) 8270 LNV = 1000 8280 ISW = -1 8290 CALL DCSTSS(ALPHA,N,BEETA,EPS,E,M,VE,LNV,ISW,IW1,W1,IERR) 8300 IF(IERR.NE.0) THEN 8310 WRITE(6,100)IERR 8320 100 FORMAT(1H ,' ERROR NO. ',I10) 8330 END IF 8340 ENERGY=E(1) 8350 E2 =E(2) 8360C-WR- WRITE(6,150)(E(J),J=1,10) 8370C 150 FORMAT(1H ,'*',10F7.3 ) 8380C-WR- IF(LK.GE. 65)WRITE(6,200)(VE(J,1),J=N-1,N) 8390C 200 FORMAT(2E30.20) 8400 RETURN 8410 END 8420CCC 8430CC 8440C 8450 SUBROUTINE INVERS(NS,NK,NST,LU,LD,TU,TD,S,SG,CON,BS,BC,NRM, 8460 & DIAG,LAMBDA,VEC) 8470 INTEGER NVEC, LROT, NS,NK,NST(0:1),LU(0:20),LD 8480 INTEGER S(0:12870,0:20,0:1), SG(0:12870,0:20,0:1) 8490 INTEGER CON(0:12870,0:20,0:1), TU(0:9999 ,0:20,0:2) 8500 INTEGER TD(0:99999, 0:2) 8510 REAL*8 DIAG(0:364364),NRM(0:364364), DTP,DTR 8520 REAL*8 LAMBDA,VEC(0:364364,0:3), BS(0:20), BC(0:20) 8530C 8540 IF( NK .EQ. 0 .OR. NK .EQ. NS/2 ) THEN 8550 NVEC=NST(0)*NST(1) 8560 ELSE 8570 NVEC=NST(0)*NST(1)*2 8580 END IF 8590C 8600 CALL VMOVE( NVEC, 0 , 3 , VEC) 8610 LROT=0 8620 GOTO 11 8630 10 CONTINUE 8640 CALL CGMETH(NS,NK,NST,LU,LD,TU,TD,S,SG,CON,BS,BC,NRM, 8650 & DIAG,LAMBDA,VEC) 8660 11 CALL VMOVE( NVEC, 1 , 0 , VEC) 8670 CALL INPRO( NVEC, 1 , 1 ,DTP, VEC) 8680 DTP=SQRT(DTP) 8690 CALL DIVER( NVEC, 1 , 1 ,DTP, VEC) 8700 CALL MULT(2,1,NS,NK,NST,LU,LD,TU,TD,S,SG,CON,BS,BC,NRM,DIAG,VEC) 8710 CALL INPRO( NVEC, 2 , 2 ,DTP, VEC) 8720 DTP=SQRT(DTP) 8730 CALL DIVER( NVEC, 2 , 2 ,DTP, VEC) 8740C 8750 IF(LAMBDA.LT.0) THEN 8760 CALL ADDER( NVEC, 3 , 2 , 1 ,1.0D0,1.0D0, VEC) 8770 ELSE 8780 CALL ADDER( NVEC, 3 , 2 , 1 ,-1.0D0,1.0D0, VEC) 8790 END IF 8800C 8810 CALL INPRO( NVEC, 3 , 3 ,DTR, VEC) 8820C-WR- WRITE(6,100)LAMBDA,DTP,SQRT(DTR) 8830 LROT=LROT+1 8840C 100 FORMAT(1H ,'LAMBDA=',E20.10,' ? ',E20.10,' D ',E20.10) 8850 IF(LROT.GE.30) THEN 8860 WRITE(6,281) 8870 281 FORMAT(1H ,'DID NOT CONVERGE IN INVERS') 8880 STOP 8890 END IF 8900 IF(SQRT(DTR).GE.0.5E- 9) GOTO 10 8910C-WR- WRITE(6,345)( VEC(K,1),K=0,MIN(NVEC-1,35) ) 8920C 345 FORMAT(1H , 5E15.7) 8930 RETURN 8940 END 8950CCC 8960CC 8970C 8980 SUBROUTINE CGMETH(NS,NK,NST,LU,LD,TU,TD,S,SG,CON,BS,BC,NRM, 8990 & DIAG,LAMBDA,VEC) 9000 INTEGER NVEC,LROT, NS,NK,NST(0:1),LU(0:20),LD 9010 INTEGER S(0:12870,0:20,0:1), SG(0:12870,0:20,0:1) 9020 INTEGER CON(0:12870,0:20,0:1), TU(0:9999 ,0:20,0:2) 9030 INTEGER TD(0:99999, 0:2) 9040 REAL*8 DIAG(0:364364), NRM(0:364364) 9050 REAL*8 LAMBDA,VEC(0:364364,0:3), BS(0:20), BC(0:20) 9060 REAL*8 DB,DR1,DR2,ALPHA,BEETA,EPSRON,DTP 9070C 9080 IF( NK .EQ. 0 .OR. NK .EQ. NS/2 ) THEN 9090 NVEC=NST(0)*NST(1) 9100 ELSE 9110 NVEC=NST(0)*NST(1)*2 9120 END IF 9130C (I) 9140 CALL INPRO( NVEC, 1 , 1 , DB, VEC ) 9150 CALL VMOVE( NVEC, 2 , 1 , VEC ) 9160 DO 10 I=0,NVEC-1 9170 VEC(I,0)=0.0 9180 10 CONTINUE 9190C (II) 9200 EPSRON=1.0E-8 9210 LROT=0 9220 CALL INPRO( NVEC, 1 , 1 ,DR1, VEC) 9230 5 CALL MULT(3,2,NS,NK,NST,LU,LD,TU,TD,S,SG,CON,BS,BC,NRM,DIAG,VEC) 9240 DO 30 I=0,NVEC-1 9250 VEC(I,3)=VEC(I,3)-LAMBDA*VEC(I,2) 9260 30 CONTINUE 9270 CALL INPRO( NVEC, 3 , 2 ,DTP, VEC) 9280 ALPHA=DR1/DTP 9290 CALL ADDER( NVEC, 0 , 0 , 2 ,1.0D0, ALPHA, VEC) 9300 CALL ADDER( NVEC, 1 , 1 , 3 ,1.0D0,-ALPHA, VEC) 9310 CALL INPRO( NVEC, 1 , 1 ,DR2, VEC) 9320C-WR- WRITE(6,999) DR2/DB 9330C 999 FORMAT(1H ,E30.20) 9340 IF(DR2/DB.LE.EPSRON) GOTO 200 9350 BEETA=DR2/DR1 9360 DR1=DR2 9370 LROT=LROT+1 9380 CALL ADDER( NVEC, 2 , 1 , 2 ,1.0D0, BEETA,VEC) 9390 IF(LROT.GE. 20) THEN 9400 WRITE(6,132) 9410 132 FORMAT(1H ,'DID NOT CONVERGE IN CG') 9420 GOTO 154 9430 END IF 9440 GOTO 5 9450 200 CONTINUE 9460 154 CONTINUE 9470 WRITE(6,678) LROT 9480 678 FORMAT(1H ,I10,'ROT IN CG') 9490 RETURN 9500 END 9510CCC 9520CC 9530C 9540 SUBROUTINE NDIST( NS, NK, NST , ST, VEC ) 9550 INTEGER ADD, NS, NK, NST(0:1), ST(0:12870,0:20,0:1), NV 9560 REAL*8 A(0:20),B(0:20), VEC(0:364364,0:3) 9570C 9580 NV= NST(0)*NST(1) 9590 DO 40 IS= 0,NS-1 9600 A(IS) = 0 9610 B(IS) = 0 9620 DO 30 I = 0,NST(0)-1 9630 DO 20 J = 0,NST(1)-1 9640 ADD = I*NST(1)+J 9650 IF( NK .EQ. 0 .OR. NK .EQ. NS/2 ) THEN 9660 A(IS) = A(IS)+VEC(ADD,1)*VEC(ADD,1)*MOD(ST(I,0,0)/2**IS,2) 9670 B(IS) = B(IS)+VEC(ADD,1)*VEC(ADD,1)*MOD(ST(J,0,1)/2**IS,2) 9680 ELSE 9690C 9700 A(IS) = A(IS)+MOD(ST(I,0,0)/2**IS,2) 9710 & *(VEC(ADD,1)*VEC(ADD,1)+VEC(NV+ADD,1)*VEC(NV+ADD,1)) 9720 B(IS) = B(IS)+MOD(ST(J,0,1)/2**IS,2) 9730 & *(VEC(ADD,1)*VEC(ADD,1)+VEC(NV+ADD,1)*VEC(NV+ADD,1)) 9740 END IF 9750C 9760 20 CONTINUE 9770 30 CONTINUE 9780 40 CONTINUE 9790 CALL DISP(NS, A,'UPHOLE') 9800 CALL DISP(NS, B,'DWHOLE') 9810 RETURN 9820 END 9830CCC 9840CC 9850C 9860 SUBROUTINE CDWORD( NS, NK, NST , ST, VEC ) 9870 INTEGER IS,ID,ADD, NS, NK, NST(0:1), ST(0:12870,0:20,0:1), NV 9880 REAL*8 S,CDW, VEC(0:364364,0:3) 9890C 9900 NV= NST(0)*NST(1) 9910 CDW= 0.0D0 9920C 9930 DO 50 IS= 1,NS-1 9940 DO 40 ID= 0,IS-1 9950 S = (-1.0D0)**( IS+ID ) 9960 DO 30 I = 0,NST(0)-1 9970 DO 20 J = 0,NST(1)-1 9980 ADD = I*NST(1)+J 9990 IF( NK .EQ. 0 .OR. NK .EQ. NS/2 ) THEN 10000 CDW = CDW + VEC(ADD,1)*VEC(ADD,1) * S 10010 & *( MOD(ST(I,0,0)/2**IS,2) + MOD(ST(J,0,1)/2**IS,2) -1 ) 10020 & *( MOD(ST(I,0,0)/2**ID,2) + MOD(ST(J,0,1)/2**ID,2) -1 ) 10030 ELSE 10040 CDW = CDW + ( VEC( ADD,1)*VEC( ADD,1) 10050 & +VEC(NV+ADD,1)*VEC(NV+ADD,1) ) * S 10060 & *( MOD(ST(I,0,0)/2**IS,2) + MOD(ST(J,0,1)/2**IS,2) -1 ) 10070 & *( MOD(ST(I,0,0)/2**ID,2) + MOD(ST(J,0,1)/2**ID,2) -1 ) 10080 END IF 10090 20 CONTINUE 10100 30 CONTINUE 10110 40 CONTINUE 10120 50 CONTINUE 10130C 10140 CDW = 2.0D0 * CDW / ( NS*(NS-1) ) 10150 WRITE(6,100) CDW 10160 100 FORMAT(1H ,'CDW ORDER=',F12.5 ) 10170C 10180 RETURN 10190 END 10200CCC 10210CC 10220C 10230 SUBROUTINE DISP(NS, A,CH) 10240 INTEGER NS 10250 CHARACTER*6 CH 10260 REAL*8 A(0:20) 10270C 10280 WRITE(6, * ) CH 10290 WRITE(6,100)( A(I),I=0,NS-1 ) 10300 100 FORMAT(1H ,8F9.6) 10310C 10320 RETURN 10330 END