マーデルング定数に関しては，このコラムで４作目となる．まず，このシリーズを書くことになった発端の話から振り返ってみたい．

　

　「食塩はＮａ+とＣｌ-が交互に並んで立方格子というジャングルジムのような構造を形成するイオン結晶である．ＮａＣｌ結晶の場合，最近接イオン間の距離をｒとすると，Ｎａ+イオンのまわりには，ｒだけ離れて６個のＣｌ-イオンがあり，第２近接には√２ｒだけ離れて１２個のＮａ+イオン，第３近接には√３ｒだけ離れて８個のＣｌ-イオン，第４近接には２ｒだけ離れて６個のＮａ+イオンがある．第５近接には√５ｒだけ離れて２４個のＣｌ-イオン，・・・という具合になる．また，異符号のイオン間には引力が働き，同符号のイオン間には斥力が働くから，ＮａＣｌ結晶の結合エネルギーは，これらのイオン間のクーロンポテンシャルを加えて，

　　Ｕ＝-e^2/r[6-12/√2+8/√3-6/2+24/√5-･･･]

となる．無限級数で示される括弧のなかの和

　　6-12/√2+8/√3-6/2+24/√5-･･･

をマーデルング定数と呼ぶ．」

　

　「ＣｓＣｌ型構造では，１個のＣｓイオンに最も接近しているのは８個のＣｌイオンである．この最近接距離をｒとすると，第２近接は（２／√３）ｒの距離にある６個のＣｓイオン，第３近接は（２√２／√３）ｒの距離にある１２個のＣｓイオン，第４近接は（√１１／√３）ｒの距離にある２４個のＣｌイオン，さらに８個のＣｓイオンが２ｒの距離（第５近接）に存在する．

　したがって，ＣｓＣｌのマーデルング級数は，

　　8-6/(2/√3)-12/(2√2/√3)+24/(√11/√3)-8/4-･･･

になる．」

　

　一般に，無限級数が収束する場合，それぞれの項のなかで絶対値が小さいものほど寄与が小さくなるから，正にせよ負にせよ絶対値の大きいものから足しあげていくのが常套であり，第ｎ近接を求めて，原点からの距離の近い順に足し合わせてみるのは妥当なところであると考えられる．したがって，たいていの物性をテーマとする物理の本には，この無限級数が収束するのは当然のごとく書かれてある．

　

　この記述をみた瞬間，小生はこの無限級数は本当に収束するのかという疑問を持った．たとえば，ＣｓＣｌのマーデルング級数を実際に計算すると

　　Ｍ＝８−５．１９６−７．３４８＋１２．５３４−４−・・・

になり，はたせるかな，この方法では計算は収束しない．直観は正しかったのである．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

プログラム（１）

　

1000 '

1010 ' *** Madelung constant ***

1020 '

1030 DEF FNODD(X)=INSTR(STR$((X-1)/2),".")

1040 '

1050 DIMENSION=3

1060 PRINT "dimension=";DIMENSION

1070 INPUT "distance =";D

1080 '

1090 MAX=(2*D+1)*(2*D+1)*(2*D+1)

1100 DIM R(MAX),R$(MAX),RR(MAX)

1110 DIM N(MAX),P(MAX),M(MAX)

1120 DIM SL(MAX),SR(MAX)

1130 D2=D*D

1140 EPS=.00001

1150 '

1160 ' ** MAIN **

1170 '

1180 *MAIN:

1190 'GOSUB *NACL

1200 'GOSUB *NACL2

1210　GOSUB *CSCL

1220 '

1230 GOSUB *SUM

1240 '

1250 GOSUB *QUICK.SORT

1260 GOSUB *TIE.RANK

1270 '

1280 PRINT:PRINT " R^2"," No."," +"," -"," M. const."

1290 SUM=0

1300 FOR I=1 TO NR

1310　 SUM=SUM+(P(I)-M(I))/SQR(RR(I))

1320　 PRINT RR(I),N(I),P(I),M(I),SUM

1330　 IF (I MOD 20)=0 THEN PRINT

1340 NEXT I

1350 '

1360 PRINT:PRINT "Madelung=";SUM

1370 PRINT "N-th neighbor=";NR

1380 END

1390 '

1400 *NACL:

1410 AA=1:BB=1:CC=1

1420 AB=0:BC=0:CA=0

1430 RETURN

1440 '

1450 *NACL2:

1460 AA=SQR(2)/2

1470 BB=AA:CC=AA

1480 AB=1/2

1490 BC=AB:CA=AB

1500 RETURN

1510 '

1520 *CSCL:

1530 AA=1:BB=1:CC=1

1540 AB=-1/3

1550 BC=AB:CA=AB

1560 RETURN

1570 '

1580 *SUM:

1590 NN=0

1600　 FOR X=-D TO D

1610　　 FOR Y=-D TO D

1620　　　 FOR Z=-D TO D

1630　　　　 PRINT "*";

1640　　　　 GOSUB *SUM3

1650　　　 NEXT Z

1660　　 NEXT Y

1670　 NEXT X

1680 PRINT

1690 RETURN

1700 '

1710 *SUM3:

1720 S1=X*X*AA*AA+Y*Y*BB*BB+Z*Z*CC*CC

1730 S2=2*(X*Y*AA*BB*AB+Y*Z*BB*CC*BC+Z*X*CC*AA*CA)

1740 R2=ABS(S1+S2):R=SQR(R2)

1750 IF R<=0 THEN 1820

1760 '

1770 IF R2>D2 THEN 1820

1780 NN=NN+1

1790 R(NN)=R2

1800 IF FNODD(X+Y+Z)=0 THEN R$(NN)="+" ELSE R$(NN)="-"

1810 '

1820 RETURN

1830 '

1840 ' *** QUICK SORT ***

1850 '

1860 *QUICK.SORT:

1870 PRINT:PRINT "ただいまソート中"

1880 '

1890 S=1:SL(1)=1:SR(1)=NN

1900 WHILE S<>0

1910　 L=SL(S):R=SR(S):S=S-1

1920　 WHILE L<R

1930　　 I=L:J=R:X=R((L+R)\2)

1940　　 WHILE I<=J

1950　　　 WHILE R(I)<X:I=I+1:WEND

1960　　　 WHILE R(J)>X:J=J-1:WEND

1970　　　 IF I<=J THEN SWAP R(I),R(J):SWAP R$(I),R$(J):I=I+1:J=J-1

1980　　 WEND

1990　　 IF I<R THEN S=S+1:SL(S)=I:SR(S)=R

2000　　 R=J

2010　 WEND

2020 WEND

2030 RETURN

2040 '

2050 ' *** TIE RANKING ***

2060 '

2070 *TIE.RANK:

2080 NR=0

2090 I=1

2100 WHILE I<=NN

2110　 L=1

2120　 IL=I+L

2130　 WHILE IL<=NN

2140　　'IF R(I)<R(IL) THEN 2190

2150　　 IF ABS((R(I)-R(IL))/R(IL))>EPS THEN 2190

2160　　 L=L+1

2170　　 IL=I+L

2180　 WEND

2190 '

2200　 P=0:M=0

2210　 FOR J=I TO I+L-1

2220　　 IF R$(J)="+" THEN P=P+1

2230　　 IF R$(J)="-" THEN M=M+1

2240　 NEXT J

2250　 NR=NR+1:RR(NR)=R(I):N(NR)=L:P(NR)=P:M(NR)=M:I=IL

2260 '

2270 WEND

2280 RETURN

　

　原点からの距離を２乗した値を配列Ｒに，引力・斥力の別を配列Ｒ＄に取り込んだ後，各データを原点からの距離の順に並べる．ソートでは最も平均的な効率がよく，最高速だといわれているクイックソートを採用した．

　

　次に，サブルーチン＊ｔｉｅ．ｒａｎｋで同一距離のイオン数を計算する．この計算では単調増加数列のプラトー（非増加部分）の長さを求めることになるが，手順を分析しそれを具体的な形に表していくという意味で，プログラム設計者の腕の見せどころとなるステップである．なお，ＢＡＳＩＣの単精度計算は１０進法換算で７桁程度であるが，表示桁数は６桁までであり，表示された数値が同じであっても，内部的には異なる数として認識される．コンピュータを使って計算するとよくこのような誤算の問題に遭遇するが，２１４０行の代わりに２１５０行のように相対誤差を計算すれば，これを回避し順位計算を正しく行うことができる．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

１）ＮａＣｌの場合

　ｄ＝４と入力すると第１４近接まで計算され，下のような計算結果が表示される．左のカラムから順に，距離２乗値，イオン数，引力イオン数，斥力イオン数，マーデルング定数の途中経過を表している．

　

　たとえば，第５近接までをとると，

　　Ｍ＝6-12/√2+8/√3-6/2+24/√5＝９．８６７　　（第５近接）

　また，最下段には第１４近接まで計算したマーデルング定数が表示される．

　　Ｍ＝−４．５０９　　（第１４近接）

　

R^2　　　　　 No.　　　　　 +　　　　　　 -　　　　　　 M. const.

1　　　　　　 6　　　　　　 6　　　　　　 0　　　　　　 6

2　　　　　　 12　　　　　　0　　　　　　 12　　　　　 -2.48528

3　　　　　　 8　　　　　　 8　　　　　　 0　　　　　　 2.13352

4　　　　　　 6　　　　　　 0　　　　　　 6　　　　　　-.86648

5　　　　　　 24　　　　　　24　　　　　　0　　　　　　 9.86665

6　　　　　　 24　　　　　　0　　　　　　 24　　　　　　.0686874

8　　　　　　 12　　　　　　0　　　　　　 12　　　　　 -4.17395

9　　　　　　 30　　　　　　30　　　　　　0　　　　　　 5.82605

10　　　　　　24　　　　　　0　　　　　　 24　　　　　 -1.76342

11　　　　　　24　　　　　　24　　　　　　0　　　　　　 5.47285

12　　　　　　8　　　　　　 0　　　　　　 8　　　　　　 3.16345

13　　　　　　24　　　　　　24　　　　　　0　　　　　　 9.81985

14　　　　　　48　　　　　　0　　　　　　 48　　　　　 -3.00868

16　　　　　　6　　　　　　 0　　　　　　 6　　　　　　-4.50868

　

Madelung=-4.50868

N-th neighbor= 14

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

２）ＣｓＣｌの場合

　同じく，ｄ＝４と入力すると第１７近接まで計算され，

　　Ｍ＝−２．９６８　　（第１７近接）

と表示される．第５近接までをとると，

　　Ｍ＝8-6/(2/√3)-12/(2√2/√3)+24/(√11/√3)-8/4

　　　＝３．９８９　　（第５近接）

　

R^2　　　　　 No.　　　　　 +　　　　　　 -　　　　　　 M. const.

1　　　　　　 8　　　　　　 8　　　　　　 0　　　　　　 8

1.33333　　　 6　　　　　　 0　　　　　　 6　　　　　　 2.80385

2.66667　　　 12　　　　　　0　　　　　　 12　　　　　 -4.54462

3.66667　　　 24　　　　　　24　　　　　　0　　　　　　 7.98897

4　　　　　　 8　　　　　　 0　　　　　　 8　　　　　　 3.98897

5.33333　　　 6　　　　　　 0　　　　　　 6　　　　　　 1.39089

6.33333　　　 24　　　　　　24　　　　　　0　　　　　　 10.9275

6.66667　　　 24　　　　　　0　　　　　　 24　　　　　　1.63237

8　　　　　　 24　　　　　　0　　　　　　 24　　　　　 -6.85292

9　　　　　　 32　　　　　　32　　　　　　0　　　　　　 3.81375

10.6667　　　 12　　　　　　0　　　　　　 12　　　　　　.139516

11.6667　　　 48　　　　　　48　　　　　　0　　　　　　 14.1925

12　　　　　　30　　　　　　0　　　　　　 30　　　　　　5.53222

13.3333　　　 24　　　　　　0　　　　　　 24　　　　　 -1.04045

14.3333　　　 24　　　　　　24　　　　　　0　　　　　　 5.2988

14.6667　　　 24　　　　　　0　　　　　　 24　　　　　 -.967996

16　　　　　　8　　　　　　 0　　　　　　 8　　　　　　-2.968

Madelung=-2.968

N-th neighbor= 17

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

　計算結果は，いずれの場合も振動が激しく，正確なマーデルング定数：

　　α＝１．７４７５６５　　（ＮａＣｌ）

　　α＝１．７６２６７５　　（ＣｓＣｌ）

とは大きくかけ離れてしまった．

　

　ＮａＣｌの場合，以下，

　　ｄ＝５　→　Ｍ＝５．２４５８　　（第２２近接）

　　ｄ＝10　→　Ｍ＝−６．０３５０　（第８５近接）

　ＣｓＣｌの場合，

　　ｄ＝５　→　Ｍ＝２．００８５　　（第２６近接）

　　ｄ＝10　→　Ｍ＝３．３９７８　　（第１０２近接）

と続く．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

　データ数は（２ｄ＋１）^3個あり，ｄをこれ以上大きくすると，メモリに入りきらないほど巨大なデータを扱うことになる．データを格納する配列の大きさはいくらでも大きくできるわけではなく，１つの配列につき６４ＫＢという上限があり，単精度実数型配列の場合は１６３８３が限界値となっている．

　

　作業用配列の使用量を最小にするため，苦肉の策であるが，以下のプログラムでは同じ機能のサブルーチンを２度繰り返して実行している．ダサイし，いやらしいことだが，経済的なメモリ使用の必要上，やむを得まい．このあたりがＢＡＳＩＣプログラムの最大の泣きどころになっている．

　

　さらに，クイック・ソートはもっとも有効な内部ソート（オンメモリソート）の１つであるが，ソートする部分の両端をスタックに記憶しておく必要があり，大きな作業用配列が必要となる．そこで，メモリ節約のために，スタックを必要とせず，しかも最高速だといわれるクイックソート並かそれ以上の速度が出せるコムソートを用いることにした．

　

プログラム（２）

　

1000 '

1010 ' *** Madelung constant ***

1020 '

1030 DEF FNODD(X)=INSTR(STR$((X-1)/2),".")

1040 '

1050 DIMENSION=3

1060 PRINT "dimension=";DIMENSION

1070 INPUT "distance =";D

1080 '

1090 D2=D*D

1100 EPS=.00001

1110 '

1120 ' ** MAIN **

1130 '

1140 *MAIN:

1150 'GOSUB *NACL

1160 'GOSUB *NACL2

1170　GOSUB *CSCL

1180 '

1190 GOSUB *DIM1

1200 DIM R(NN),R$(NN)

1210 GOSUB *SUM

1220 'GOSUB *QUICK.SORT

1230 GOSUB *COMB.SORT

1240 '

1250 GOSUB *DIM2

1260 DIM RR(NR)

1270 DIM N(NR),P(NR),M(NR)

1280 GOSUB *TIE.RANK

1290 '

1300 PRINT:PRINT " R^2"," No."," +"," -"," M. const."

1310 SUM=0

1320 FOR I=1 TO NR

1330　 SUM=SUM+(P(I)-M(I))/SQR(RR(I))

1340　 PRINT RR(I),N(I),P(I),M(I),SUM

1350　 IF (I MOD 20)=0 THEN PRINT

1360 NEXT I

1370 '

1380 PRINT:PRINT "Madelung=";SUM

1390 PRINT "N-th neighbor=";NR

1400 END

1410 '

1420 *NACL:

1430 AA=1:BB=1:CC=1

1440 AB=0:BC=0:CA=0

1450 RETURN

1460 '

1470 *NACL2:

1480 AA=SQR(2)/2

1490 BB=AA:CC=AA

1500 AB=1/2

1510 BC=AB:CA=AB

1520 RETURN

1530 '

1540 *CSCL:

1550 AA=1:BB=1:CC=1

1560 AB=-1/3

1570 BC=AB:CA=AB

1580 RETURN

1590 '

1600 *SUM:

1610 NN=0

1620　 FOR X=-D TO D

1630　　 FOR Y=-D TO D

1640　　　 FOR Z=-D TO D

1650　　　　 PRINT "*";

1660　　　　 GOSUB *SUM3

1670　　　 NEXT Z

1680　　 NEXT Y

1690　 NEXT X

1700 PRINT

1710 RETURN

1720 '

1730 *SUM3:

1740 S1=X*X*AA*AA+Y*Y*BB*BB+Z*Z*CC*CC

1750 S2=2*(X*Y*AA*BB*AB+Y*Z*BB*CC*BC+Z*X*CC*AA*CA)

1760 R2=ABS(S1+S2):R=SQR(R2)

1770 IF R<=0 THEN 1840

1780 '

1790 IF R2>D2 THEN 1840

1800 NN=NN+1

1810 R(NN)=R2

1820 IF FNODD(X+Y+Z)=0 THEN R$(NN)="+" ELSE R$(NN)="-"

1830 '

1840 RETURN

1850 '

1860 ' *** QUICK SORT ***

1870 '

1880 *QUICK.SORT:

1890 PRINT:PRINT "ただいまソート中"

1900 '

1910 S=1:SL(1)=1:SR(1)=NN

1920 WHILE S<>0

1930　 L=SL(S):R=SR(S):S=S-1

1940　 WHILE L<R

1950　　 I=L:J=R:X=R((L+R)\2)

1960　　 WHILE I<=J

1970　　　 WHILE R(I)<X:I=I+1:WEND

1980　　　 WHILE R(J)>X:J=J-1:WEND

1990　　　 IF I<=J THEN SWAP R(I),R(J):SWAP R$(I),R$(J):I=I+1:J=J-1

2000　　 WEND

2010　　 IF I<R THEN S=S+1:SL(S)=I:SR(S)=R

2020　　 R=J

2030　 WEND

2040 WEND

2050 RETURN

2060 '

2070 ' *** TIE RANKING ***

2080 '

2090 *TIE.RANK:

2100 NR=0

2110 I=1

2120 WHILE I<=NN

2130　 L=1

2140　 IL=I+L

2150　 WHILE IL<=NN

2160　　 IF ABS((R(I)-R(IL))/R(IL))>EPS THEN 2200

2170　　 L=L+1

2180　　 IL=I+L

2190　 WEND

2200 '

2210　 P=0:M=0

2220　 FOR J=I TO I+L-1

2230　　 IF R$(J)="+" THEN P=P+1

2240　　 IF R$(J)="-" THEN M=M+1

2250　 NEXT J

2260　 NR=NR+1:RR(NR)=R(I):N(NR)=L:P(NR)=P:M(NR)=M:I=IL

2270 '

2280 WEND

2290 RETURN

2300 '

2310 ' *** dimension-1 ***

2320 '

2330 *DIM1:

2340 NN=0

2350　 FOR X=-D TO D

2360　　 FOR Y=-D TO D

2370　　　 FOR Z=-D TO D

2380　　　　 PRINT "*";

2390　　　　 GOSUB *DIM3

2400　　　 NEXT Z

2410　　 NEXT Y

2420　 NEXT X

2430 PRINT

2440 RETURN

2450 '

2460 *DIM3:

2470 S1=X*X*AA*AA+Y*Y*BB*BB+Z*Z*CC*CC

2480 S2=2*(X*Y*AA*BB*AB+Y*Z*BB*CC*BC+Z*X*CC*AA*CA)

2490 R2=ABS(S1+S2):R=SQR(R2)

2500 IF R<=0 THEN 2530

2510 IF R2>D2 THEN 2530

2520 NN=NN+1

2530 RETURN

2540 '

2550 ' *** dimension-2 ***

2560 '

2570 *DIM2:

2580 NR=0

2590 I=1

2600 WHILE I<=NN

2610　 L=1

2620　 IL=I+L

2630　 WHILE IL<=NN

2640　　 IF ABS((R(I)-R(IL))/R(IL))>EPS THEN 2680

2650　　 L=L+1

2660　　 IL=I+L

2670　 WEND

2680 '

2690　 NR=NR+1:I=IL

2700 '

2710 WEND

2720 RETURN

2730 '

2740 ' *** COMB SORT-11 ***

2750 '

2760 *COMB.SORT:

2770 PRINT:PRINT "ただいまソート中"

2780 '

2790 SF=1.3

2800 GAP=NN

2810 SW=1

2820 WHILE SW=1 OR GAP>1

2830　 IF INT(GAP/SF)>1 THEN GAP=INT(GAP/SF) ELSE GAP=1

2840　 IF GAP=9 OR GAP=10 THEN GAP=11

2850　 SW=0

2860　 FOR I=1 TO NN-GAP

2870　　 J=I+GAP

2880　　 IF R(I)>R(J) THEN SWAP R(I),R(J):SWAP R$(I),R$(J):SW=1

2890　 NEXT I

2900 WEND

2910 RETURN

　

　しかし，このプログラムを用いても，

１）ＮａＣｌの場合

　　ｄ＝15　→　Ｍ＝１０．０５３９　（第１８９近接）

２）ＣｓＣｌの場合

　　ｄ＝14　→　Ｍ＝５．５９２９　　（第１９８近接）

まで計算するのが限界であった．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

　さらに，以下のプログラムでは，メモリの節約効果を高めるように工夫してある．１０カ所の変更箇所には＊マークをつけておいた．

　

　主たる変更点を述べると，サブルーチン＊ＳＵＭでは，データの取り込み範囲をｘ≧０とし，サブルーチン＊ＳＵＭ３の因子ＭをＭ＝２としている．１８２０行でこの値を配列Ｌに取り込むが，これによって取り込みに要する時間もメモリ使用量も約１／２に節減される．

　

プログラム（３）

　

1000 '

1010 ' *** Madelung constant ***

1020 '

1030 DEF FNODD(X)=INSTR(STR$((X-1)/2),".")

1040 '

1050 DIMENSION=3

1060 PRINT "dimension=";DIMENSION

1070 INPUT "distance =";D

1080 '

1090 D2=D*D

1100 EPS=.00001

1110 '

1120 ' ** MAIN **

1130 '

1140 *MAIN:

1150 'GOSUB *NACL

1160 'GOSUB *NACL2

1170　GOSUB *CSCL

1180 '

1190 GOSUB *DIM1

1200 DIM R(NN),R$(NN),L(NN)　　　＊

1210 GOSUB *SUM

1220 'GOSUB *QUICK.SORT

1230 GOSUB *COMB.SORT

1240 '

1250 GOSUB *DIM2

1260 DIM RR(NR)

1270 DIM N(NR),P(NR),M(NR)

1280 GOSUB *TIE.RANK

1290 '

1300 PRINT:PRINT " R^2"," No."," +"," -"," M. const."

1310 SUM=0

1320 FOR I=1 TO NR

1330　 SUM=SUM+(P(I)-M(I))/SQR(RR(I))

1340　 PRINT RR(I),N(I),P(I),M(I),SUM

1350　 IF (I MOD 20)=0 THEN PRINT

1360 NEXT I

1370 '

1380 PRINT:PRINT "Madelung=";SUM

1390 PRINT "N-th neighbor=";NR

1400 END

1410 '

1420 *NACL:

1430 AA=1:BB=1:CC=1

1440 AB=0:BC=0:CA=0

1450 RETURN

1460 '

1470 *NACL2:

1480 AA=SQR(2)/2

1490 BB=AA:CC=AA

1500 AB=1/2

1510 BC=AB:CA=AB

1520 RETURN

1530 '

1540 *CSCL:

1550 AA=1:BB=1:CC=1

1560 AB=-1/3

1570 BC=AB:CA=AB

1580 RETURN

1590 '

1600 *SUM:

1610 NN=0

1620　 FOR X=0 TO D　　　＊

1630　　 FOR Y=-D TO D

1640　　　 FOR Z=-D TO D

1650　　　　 PRINT "*";

1660　　　　 GOSUB *SUM3

1670　　　 NEXT Z

1680　　 NEXT Y

1690　 NEXT X

1700 PRINT

1710 RETURN

1720 '

1730 *SUM3:

1740 S1=X*X*AA*AA+Y*Y*BB*BB+Z*Z*CC*CC

1750 S2=2*(X*Y*AA*BB*AB+Y*Z*BB*CC*BC+Z*X*CC*AA*CA)

1760 R2=ABS(S1+S2):R=SQR(R2)

1770 IF R<=0 THEN 1850

1780 '

1790 IF R2>D2 THEN 1850

1800 M=2:IF X=0 THEN M=1　　　＊

1810 NN=NN+1

1820 R(NN)=R2:L(NN)=M　　　＊

1830 IF FNODD(X+Y+Z)=0 THEN R$(NN)="+" ELSE R$(NN)="-"

1840 '

1850 RETURN

1860 '

1870 ' *** QUICK SORT ***

1880 '

1890 *QUICK.SORT:

1900 PRINT:PRINT "ただいまソート中"

1910 '

1920 S=1:SL(1)=1:SR(1)=NN

1930 WHILE S<>0

1940　 L=SL(S):R=SR(S):S=S-1

1950　 WHILE L<R

1960　　 I=L:J=R:X=R((L+R)\2)

1970　　 WHILE I<=J

1980　　　 WHILE R(I)<X:I=I+1:WEND

1990　　　 WHILE R(J)>X:J=J-1:WEND

2000　　　 IF I<=J THEN SWAP R(I),R(J):SWAP R$(I),R$(J)　　　　　　　　　　　　　　　　　　　　　　　　:SWAP L(I),L(J):I=I+1:J=J-1　　　＊

2010　　 WEND

2020　　 IF I<R THEN S=S+1:SL(S)=I:SR(S)=R

2030　　 R=J

2040　 WEND

2050 WEND

2060 RETURN

2070 '

2080 ' *** TIE RANKING ***

2090 '

2100 *TIE.RANK:

2110 NR=0

2120 I=1

2130 WHILE I<=NN

2140　 L=1

2150　 IL=I+L

2160　 WHILE IL<=NN

2170　　 IF ABS((R(I)-R(IL))/R(IL))>EPS THEN 2210

2180　　 L=L+1

2190　　 IL=I+L

2200　 WEND

2210 '

2220　 P=0:M=0

2230　 FOR J=I TO I+L-1

2240　　 IF R$(J)="+" THEN P=P+L(J)　　　＊

2250　　 IF R$(J)="-" THEN M=M+L(J)　　　＊

2260　 NEXT J

2270　 NR=NR+1:RR(NR)=R(I):N(NR)=P+M:P(NR)=P:M(NR)=M:I=IL　　　＊

2280 '

2290 WEND

2300 RETURN

2310 '

2320 ' *** dimension-1 ***

2330 '

2340 *DIM1:

2350 NN=0

2360　 FOR X=0 TO D　　　＊

2370　　 FOR Y=-D TO D

2380　　　 FOR Z=-D TO D

2390　　　　 PRINT "*";

2400　　　　 GOSUB *DIM3

2410　　　 NEXT Z

2420　　 NEXT Y

2430　 NEXT X

2440 PRINT

2450 RETURN

2460 '

2470 *DIM3:

2480 S1=X*X*AA*AA+Y*Y*BB*BB+Z*Z*CC*CC

2490 S2=2*(X*Y*AA*BB*AB+Y*Z*BB*CC*BC+Z*X*CC*AA*CA)

2500 R2=ABS(S1+S2):R=SQR(R2)

2510 IF R<=0 THEN 2540

2520 IF R2>D2 THEN 2540

2530 NN=NN+1

2540 RETURN

2550 '

2560 ' *** dimension-2 ***

2570 '

2580 *DIM2:

2590 NR=0

2600 I=1

2610 WHILE I<=NN

2620　 L=1

2630　 IL=I+L

2640　 WHILE IL<=NN

2650　　 IF ABS((R(I)-R(IL))/R(IL))>EPS THEN 2690

2660　　 L=L+1

2670　　 IL=I+L

2680　 WEND

2690 '

2700　 NR=NR+1:I=IL

2710 '

2720 WEND

2730 RETURN

2740 '

2750 ' *** COMB SORT-11 ***

2760 '

2770 *COMB.SORT:

2780 PRINT:PRINT "ただいまソート中"

2790 '

2800 SF=1.3

2810 GAP=NN

2820 SW=1

2830 WHILE SW=1 OR GAP>1

2840　 IF INT(GAP/SF)>1 THEN GAP=INT(GAP/SF) ELSE GAP=1

2850　 IF GAP=9 OR GAP=10 THEN GAP=11

2860　 SW=0

2870　 FOR I=1 TO NN-GAP

2880　　 J=I+GAP

2890　　 IF R(I)>R(J) THEN SWAP R(I),R(J):SWAP R$(I),R$(J)　　　　　　　　　　　　　　　　　　　　　　　　:SWAP L(I),L(J):SW=1　　　＊

2900　 NEXT I

2910 WEND

2920 RETURN

　

　これにより，

１）ＮａＣｌの場合

　　ｄ＝19　→　Ｍ＝−１．１２９９　（第３０３近接）

２）ＣｓＣｌの場合

　　ｄ＝15　→　Ｍ＝１２．６２４　　（第２２６近接）

　　ｄ＝18　→　Ｍ＝２．６７５４　　（第３２６近接）

まで計算範囲が拡大している．配列の６４ＫＢの制約の中では，これが限界であろう．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

　今回のコラムでは，なかなか収束しない計算を展望してみた．途中，MSDOS/N88BASICによる制約のために，本題から一時脱線して，メモリ節減というアリ地獄にはまってしまった．確かにメモリの節減は大切なことであるが，いまはその辺はＯＳが面倒をみてくれて，プログラマはメモリサイズを気にしなくてもよい時代になっている．

　

　マーデルング級数の困難な点は，級数のなかの陽イオンの項と陰イオンの項が別々に発散することである．すなわち，∞−∞の不定形であり，条件つき収束のため，マーデルング級数の収束はきわめて遅く，プログラム的には収束しないという結論になってしまうことである．

　

　マーデルング定数の計算法に関しては，固体物理の教科書（例えば，キッテル「固体物理学入門」丸善・・・現在，第７版が出版されている）を参考にしたが，そこには，フーリエ変換を用いて計算する方法が紹介されている．特異点を打ち消し無限大になるのを抑えながら計算する非常に巧みな方法なのであろうが，小生にとってはまるでトリックでも見せられている感じで，あまり感心しなかった．結晶の規則的な配列は，そのような難解な方法を使わなくても，初等的な方法で十分計算可能と思われたからである．

　

　たとえ収束は遅くとも，いろいろな結晶構造のマーデルング定数を求めるプログラムを組み立ててみたい．おぼろげな概念であるが，これがこのシリーズを書くことになった動機であり，目的となっている．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝