実行時内部エラーのご報告

 投稿者：mike  投稿日：2017年 2月15日(水)10時32分14秒
  下記のプログラムの実行中に内部エラーが発生しました。
EXTERNAL FUNCTION getZminApex(sc,xShift,yShift)
のデバグのためprintで変数の変化を見ながら実行中にエラーが発生しました。
立方体が１回転する前に発生しました。
実行環境は　十進BASIC 0.6.6.0 / macOS 10.7.5　/　2.7Ghz Intel Core i7
メモリ 4GB 1333Mhz DDR3　です。
ご検討ください。

!
! ========= steepest descent method 3D ==========
!
! eigenFunctionSD2D.bas
! Mitsuru Ikeuchi (C) Copyleft
!
! ver 0.0.0   2017.01.31  created
!
OPTION ARITHMETIC NATIVE
DECLARE EXTERNAL SUB sd3d.setInitialCondition,sd3d.SDiteration,sd3d.drawState
DECLARE EXTERNAL FUNCTION sd3d.iterationCount,sd3d.stateEnergy,INKEY$
LET Ax = -PI/24
LET Ay = -PI/12
LET drawMode = 0
!setInitialCondition(stateMax,kx0,ky0)
CALL setInitialCondition(10,2.0,2.0)
!
LET ist = 0
FOR it=1 TO 1000
   LET Ay=Ay +PI/180
   LET dump = 0.05
   CALL SDiteration(6, 2, dump) !(stateMax, iterMax, dump)
   CALL drawState(ist,Ax,Ay,drawMode) !(state,xRotateAngle,yRotateAngle,drawMode)
   LET S$=INKEY$
   IF S$="0" THEN LET ist = MOD(ist+1,6)
   IF S$="1" THEN LET drawMode = MOD(drawMode+1,3)
   IF S$="." THEN EXIT FOR
NEXT it
!
END


EXTERNAL FUNCTION INKEY$
OPTION ARITHMETIC NATIVE
SET ECHO "OFF"
LET S$=""
CHARACTER INPUT NOWAIT: S$
LET INKEY$=S$
END FUNCTION
!
! ---------- steepest descent method 3D ----------
!
! system Hamiltonian: H = -delta/2 + V(r) , delta r = div grad r
! eigen energy Ei, eigen function set { |i> }
!
! procedure : successive approximation
! (i) trial function set { |0>,|1>,..,|i>,.. }
! (2) energy of |i> : ei = <i|H|i>/<i|i>
! (3) steepest gradient direction (H-ei)|i>
! (4) next generation : |i(next)> = |i> - dumpingFactor*(H-ei)|i>
! (5) orthogonalization { |0>,|1>,..,|i>,.. }  (Gram-Schmidt)
! (6) goto (2)
!
MODULE sd3d
MODULE OPTION ARITHMETIC NATIVE
OPTION BASE 0
PUBLIC SUB setInitialCondition, SDiteration, drawState
PUBLIC FUNCTION iterationCount, stateEnergy
SHARE NUMERIC NNx, NNy, NNz, dx, dy, dz, iterCount
SHARE NUMERIC cosAx,sinAx,cosAy,sinAy,cx0,cy0,cz0 !--- 3D graphics
SHARE NUMERIC sdEnergy(10)           ! electron state energy
SHARE NUMERIC sdState(10,65,65,65)   ! electron states 0...20 0:ground state
SHARE NUMERIC wrk(65,65,65)          ! state work space in steepestDescent
SHARE NUMERIC vv(65,65,65)           ! external potential
SHARE NUMERIC boxApex(0 TO 7,0 TO 2) ! boxApex(i,j) i-th box apex j-th coordinate x,y,z
SHARE NUMERIC boxEdge(0 TO 11,0 TO 2)! boxEdge(i,j) i-th edge j=0,1:j-th apex, j=2:for marking
LET NNx = 32                         ! x-max number of sdState(,NNx,NNy,NNz)
LET NNy = 32                         ! y-max number of sdState(,NNx,NNy,NNz)
LET NNz = 32                         ! z-max number of sdState(,NNx,NNy,NNz)
LET dx = 0.5                         ! (au) x division
LET dy = 0.5                         ! (au) y division
LET dz = 0.5                         ! (au) y division
LET iterCount = 0                    ! sd iteration count
!
! ---------- set initial condition
!
EXTERNAL SUB setInitialCondition(stateMax,kx0,ky0)   !public
   DECLARE EXTERNAL SUB setInitialState,setPotential,setBox
   RANDOMIZE
   CALL setInitialState(stateMax)
   CALL setPotential(8,8,8)
   CALL setBox
   ! set window
   LET xMargin = 60
   LET yMargin = 120
   SET WINDOW -xMargin,500-xMargin,-yMargin,500-yMargin
END SUB
!
EXTERNAL SUB setInitialState(stateMax)
   RANDOMIZE
   FOR ist=0 TO stateMax-1
      FOR i=1 TO NNx-2
         FOR j=1 TO NNy-2
            FOR k=1 TO NNz-2
               LET sdState(ist,i,j,k) = RND-0.5
            NEXT k
         NEXT j
      NEXT i
      CALL normalizeState(ist)
   NEXT ist
END SUB
!
EXTERNAL SUB setPotential(x0,y0,z0)
   FOR i=0 TO NNx-1
      FOR j=0 TO NNy-1
         FOR k=0 TO NNz-1
            LET x = i*dx
            LET y = j*dy
            LET z = k*dz
            LET r = SQR((x-x0)*(x-x0)+(y-y0)*(y-y0)+(z-z0)*(z-z0))
            IF r<0.25 THEN LET r = 0.25
            LET vv(i,j,k) = -1/r
         NEXT k
      NEXT j
   NEXT i
END SUB

EXTERNAL SUB setPotential2(x0,y0,z0)
   FOR i=0 TO NNx-1
      FOR j=0 TO NNy-1
         FOR k=0 TO NNz-1
            LET x = i*dx
            LET y = j*dy
            LET z = k*dz
            LET vv(i,j,k) = 2.0*((x-x0)*(x-x0)+(y-y0)*(y-y0)+(z-z0)*(z-z0))
         NEXT k
      NEXT j
   NEXT i
END SUB

EXTERNAL SUB setBox
   IF boxApex(7,2)<>1 THEN
      DATA 0,0,0,  1,0,0,  0,1,0,  1,1,0,  0,0,1,  1,0,1,  0,1,1,  1,1,1
      MAT READ boxApex !(0 TO 7,0 TO 2)
      DATA 0,1,9,  0,2,9,  0,4,9,  1,3,9,  1,5,9,  2,3,9,  2,6,9,  3,7,9,   4,5,9,  4,6,9,  5,7,9,  6,7,9
      MAT READ boxEdge !(0 TO 11,0 TO 2)
   END IF
END SUB
!
! ---------- steepest descent iteration
!
EXTERNAL SUB SDiteration(stateMax, iterMax, dump) !public
   DECLARE EXTERNAL FUNCTION steepestDescent
   DECLARE EXTERNAL SUB GramSchmidt
   FOR i=0 TO iterMax-1
      FOR ist=0 TO stateMax-1
         LET sdEnergy(ist) = steepestDescent(ist, dump)
      NEXT ist
      CALL GramSchmidt(stateMax)
      LET iterCount = iterCount + 1
   NEXT i
END SUB
!
EXTERNAL FUNCTION steepestDescent(ist,dump)  !---  steepest descent method
   DECLARE EXTERNAL FUNCTION energyOfState
   DECLARE EXTERNAL SUB normalizeState
   LET h2 = 2*dx*dx
   LET ei = energyOfState(ist)               !---  E_ist = <ist|H|ist>
   FOR i=1 TO NNx-2                          !---  |wrk> = (H-E_ist)|ist>
      FOR j=1 TO NNy-2
         FOR k=1 TO NNz-2
            LET kesdState = (6*sdState(ist,i,j,k)-sdState(ist,i+1,j,k)-sdState(ist,i-1,j,k)&
&-sdState(ist,i,j+1,k)-sdState(ist,i,j-1,k)-sdState(ist,i,j,k+1)-sdState(ist,i,j,k-1))/h2
            LET wrk(i,j,k) = kesdState+(vv(i,j,k)-ei)*sdState(ist,i,j,k)
         NEXT k
      NEXT j
   NEXT i
   FOR i=1 TO NNx-2                          !---  |ist(next)> = |ist> - dump*|wrk>  ( norm(|ist(next)>) <>1 )
      FOR j=1 TO NNy-2
         FOR k=1 TO NNz-2
            LET sdState(ist,i,j,k) = sdState(ist,i,j,k)-dump*wrk(i,j,k)
         NEXT k
      NEXT j
   NEXT i
   CALL normalizeState(ist)
   LET steepestDescent = ei
END FUNCTION
!
EXTERNAL FUNCTION energyOfState(ist) !--- E_ist = <ist|H|ist>/<ist|ist>
   LET h2 = 2*dx*dx
   LET s = 0
   LET sn=0
   FOR i=1 TO NNx-2
      FOR j=1 TO NNy-2
         FOR k=1 TO NNz-2
            LET kesdState = (6*sdState(ist,i,j,k)-sdState(ist,i+1,j,k)-sdState(ist,i-1,j,k)&
&-sdState(ist,i,j+1,k)-sdState(ist,i,j-1,k)-sdState(ist,i,j,k+1)-sdState(ist,i,j,k-1))/h2
            LET s = s + sdState(ist,i,j,k)*(kesdState+vv(i,j,k)*sdState(ist,i,j,k))
            LET sn = sn+sdState(ist,i,j,k)*sdState(ist,i,j,k)
         NEXT k
      NEXT j
   next i
   LET energyOfState = s/sn
END FUNCTION
!
EXTERNAL SUB GramSchmidt(stateMax)  !--- Gram-Schmidt orthonormalization
   DECLARE EXTERNAL FUNCTION innerProduct
   DECLARE EXTERNAL SUB normalizeState
   CALL normalizeState(0)
   FOR istate=1 TO stateMax-1
      FOR ist=0 TO istate-1
         LET s = innerProduct(ist,istate)
         FOR i=1 TO NNx-2
            FOR j=1 TO NNy-2
               FOR k=1 TO NNz-2
                  LET sdState(istate,i,j,k) = sdState(istate,i,j,k) - s*sdState(ist,i,j,k)
               NEXT k
            NEXT j
         NEXT i
      NEXT ist
      CALL normalizeState(istate)
   NEXT iState
END SUB
!
EXTERNAL FUNCTION innerProduct(ist,jst)  !--- <ist|jst>
   LET s = 0
   FOR i=1 TO NNx-2
      FOR j=1 TO NNy-2
         FOR k=1 TO NNz-2
            LET s = s + sdState(ist,i,j,k)*sdState(jst,i,j,k)
         NEXT k
      NEXT j
   NEXT i
   LET innerProduct = s*dx*dy*dz
END FUNCTION
!
EXTERNAL SUB normalizeState(ist)
   LET s = 0
   FOR i=1 TO NNx-2
      FOR j=1 TO NNy-2
         FOR k=1 TO NNz-2
            LET s = s + sdState(ist,i,j,k)*sdState(ist,i,j,k)
         NEXT k
      NEXT j
   NEXT i
   LET a = SQR(1/(s*dx*dy*dz))
   FOR i=1 TO NNx-2
      FOR j=1 TO NNy-2
         FOR k=1 TO NNz-2
            LET sdState(ist,i,j,k) = a*sdState(ist,i,j,k)
         NEXT k
      NEXT j
   NEXT i
END SUB
!
! ---------- utility
!
EXTERNAL FUNCTION iterationCount
   LET iterationCount = iterCount
END FUNCTION
!
EXTERNAL FUNCTION stateEnergy(ist)
   LET stateEnergy = sdEnergy(ist)
END FUNCTION
!
! ---------- 3D graphics
!
EXTERNAL SUB setRotateXYParameters(angleX,angleY,xCenter,yCenter,zCenter)
   LET cosAx = COS(angleX)
   LET sinAx = SIN(angleX)
   LET cosAy = COS(angleY)
   LET sinAy = SIN(angleY)
   LET cx0 = xCenter
   LET cy0 = yCenter
   LET cz0 = zCenter
END SUB

EXTERNAL SUB plotLines3D(x,y,z,xShift,yShift)
   LET x1 = cosAy*(x-cx0)+sinAy*(sinAx*(y-cy0)+cosAx*(z-cz0))
   LET y1 = cosAx*(y-cy0)-sinAx*(z-cz0)
   LET z1 =-sinAy*(x-cx0)+cosAy*(sinAx*(y-cy0)+cosAx*(z-cz0))
   PLOT LINES: x1+cx0+xShift, y1+cy0+yShift; !z=z1+cz0
END SUB

EXTERNAL SUB drawDisk3D(x,y,z,r,xShift,yShift)
   LET x1 = cosAy*(x-cx0)+sinAy*(sinAx*(y-cy0)+cosAx*(z-cz0))
   LET y1 = cosAx*(y-cy0)-sinAx*(z-cz0)
   LET z1 =-sinAy*(x-cx0)+cosAy*(sinAx*(y-cy0)+cosAx*(z-cz0))
   DRAW disk WITH SCALE(r)*SHIFT(x1+cx0+xShift,y1+cy0+yShift)
END SUB

EXTERNAL SUB plotText3DAt(x,y,z,xShift,yShift,A$)
   LET x1 = cosAy*(x-cx0)+sinAy*(sinAx*(y-cy0)+cosAx*(z-cz0))
   LET y1 = cosAx*(y-cy0)-sinAx*(z-cz0)
   LET z1 =-sinAy*(x-cx0)+cosAy*(sinAx*(y-cy0)+cosAx*(z-cz0))
   PLOT TEXT,AT x1+cx0+xShift,y1+cy0+yShift:A$
END SUB

EXTERNAL FUNCTION getRotateZ(x,y,z,xShift,yShift)
   LET x1 = cosAy*(x-cx0)+sinAy*(sinAx*(y-cy0)+cosAx*(z-cz0))
   LET y1 = cosAx*(y-cy0)-sinAx*(z-cz0)
   LET z1 =-sinAy*(x-cx0)+cosAy*(sinAx*(y-cy0)+cosAx*(z-cz0))
   LET getRotateZ = z1+cz0
END FUNCTION
!
! ---------- drawState
!
EXTERNAL SUB drawState(ist,xRotateAngle,yRotateAngle,drawMode)
!DECLARE EXTERNAL SUB setRotateXYParameters,drawDensity3D,drawStateGrid,drawStateGridx6
   SET DRAW MODE HIDDEN
   CLEAR
   LET sc = 6
   IF drawMode=0 THEN
      CALL drawDensity3D(ist,sc, xRotateAngle,yRotateAngle)
   ELSEIF drawMode=1 THEN
      CALL setRotateXYParameters(-PI/3,-PI/6,(NNx/2)*sc,(NNy/2)*sc,(NNy/2)*sc)
      CALL drawStateGrid(ist,320/NNx,1.0,0,50) !(ist,sc,zMag,xShift,yShift)
   ELSEIF drawMode=2 THEN
      CALL setRotateXYParameters(-PI/3,-PI/6,(NNx/2)*sc,(NNy/2)*sc,(NNz/2)*sc)
      CALL drawStateGridx6
   ELSE
   !
   END IF
   !--- caption
   SET TEXT HEIGHT 10
   SET TEXT COLOR 1 ! black
   PLOT TEXT, AT 0,-50 ,USING "iterarion count =##### ":iterCount
   PLOT TEXT, AT 0,-65 ,USING "E0 =###.######### E3 =###.#########(au)":sdEnergy(0),sdEnergy(3)
   PLOT TEXT, AT 0,-80 ,USING "E1 =###.######### E4 =###.#########(au)":sdEnergy(1),sdEnergy(4)
   PLOT TEXT, AT 0,-95 ,USING "E2 =###.######### E5 =###.#########(au)":sdEnergy(2),sdEnergy(5)
   PLOT TEXT, AT 0,-110 :"steepest descent method 2D"
   SET DRAW MODE EXPLICIT
END SUB
!
EXTERNAL SUB drawDensity3D(ist,sc,xRotateAngle,yRotateAngle) !----- drawMode=0
   DECLARE EXTERNAL SUB setRotateXYParameters, drawDisk3D, plotEdge3D
   DECLARE EXTERNAL FUNCTION getRotateZ,getZminApex
   CALL setRotateXYParameters(xRotateAngle,yRotateAngle,NNx*sc/2,NNy*sc/2,NNx*sc/2)
   LET xOffset = 80
   LET yOffset = 50
   SET LINE COLOR 1
   FOR i=1 TO NNx-2
      FOR j=1 TO NNy-2
         FOR k=1 TO NNz-2
            LET psi = sdState(ist,i,j,k)
            LET psi2 = psi*psi
            IF psi2>0.001 THEN
               IF psi>0 THEN SET AREA COLOR 4 ELSE SET AREA COLOR 2
               CALL drawDisk3D(i*sc,j*sc,k*sc,ABS(psi)*10,xOffset,yOffset)
            END IF
         NEXT k
      NEXT j
   NEXT i
   CALL plotEdge3D(0,sc,xOffset,yOffset,1) !x-axis:black !(iEdge,xShift,yShift,lineColorIndex)
   CALL plotEdge3D(1,sc,xOffset,yOffset,2) !y-axis:blue
   CALL plotEdge3D(2,sc,xOffset,yOffset,3) !z-axis:green
   FOR iEdge=3 TO 11
      CALL plotEdge3D(iEdge,sc,xOffset,yOffset,8) !other edge:grey
   NEXT iEdge
   CALL plotText3DAt(0,0,NNz*sc,xOffset,yOffset,"Z"&STR$(getRotateZ(0,0,NNz*sc,xOffset,yOffset)))
   CALL plotText3DAt(0,0,0,xOffset,yOffset,"Z"&STR$(getRotateZ(0,0,0,xOffset,yOffset)))
   LET apx = getZminApex(sc,xOffset,yOffset)
   CALL drawDisk3D(boxApex(apx,0)*NNx*sc,boxApex(apx,1)*NNy*sc,boxApex(apx,2)*NNz*sc,10,xOffset,yOffset)
   SET TEXT COLOR 1
   SET TEXT HEIGHT 10
   PLOT TEXT, AT 0,-20:"|"&STR$(ist)&">"
   PLOT TEXT, AT 0,-35 ,USING "Ax =###.#(deg) Ay =###.#(deg)":MOD(xRotateAngle*180/PI,360),MOD(yRotateAngle*180/PI,360)
END SUB
!
EXTERNAL SUB plotEdge3D(iEdge,sc,xShift,yShift,lineColorIndex)
   DECLARE EXTERNAL SUB plotLines3D
   SET LINE COLOR lineColorIndex
   FOR j=0 TO 1
      LET apex = boxEdge(iEdge,j)
      LET x = boxApex(apex,0)*NNx*sc
      LET y = boxApex(apex,1)*NNy*sc
      LET z = boxApex(apex,2)*NNz*sc
      CALL plotLines3D(x,y,z,xShift,yShift)
   NEXT j
   PLOT LINES
END SUB

EXTERNAL FUNCTION getZminApex(sc,xShift,yShift)
   DECLARE EXTERNAL FUNCTION getRotateZ !(x,y,z,xShift,yShift)
   LET zmin = 1E99
   FOR iApex = 0 TO 7
      LET x = boxApex(iApex,0)*NNx*sc
      LET y = boxApex(iApex,1)*NNy*sc
      LET z = boxApex(iApex,2)*NNz*sc
      LET rotz = getRotateZ(x,y,z,xShift,yShift)
      PRINT USING "###   rotZ=####.##   x=###.##  y=###.##   z=###.##":iApex,rotz,x,y,z
      IF rotz<zmin THEN
         LET zmin = rotz
         LET zminApex = iApex
      END IF
   NEXT iApex
   LET getZminApex = zminApex
END FUNCTION
!
EXTERNAL SUB drawStateGrid(ist,sc,zMag,xShift,yShift)  !--- drawMode=1
   DECLARE EXTERNAL SUB drawGridXY
   CALL drawGridXY(ist,sc,zMag,xShift,yShift)
   SET TEXT COLOR 1
   SET TEXT HEIGHT 10
   PLOT TEXT, AT 0,-20:"|"&STR$(ist)&">"
END SUB
!
EXTERNAL SUB drawStateGridx6
   DECLARE EXTERNAL SUB drawGridXY
   CALL drawGridXY(0,4,0.5,-30,40) !(ist,sc,zMag,xShift,yShift)  !--- drawMode=2
   CALL drawGridXY(1,4,0.5,120,40)
   CALL drawGridXY(2,4,0.5,270,40)
   CALL drawGridXY(3,4,0.5,-30,190)
   CALL drawGridXY(4,4,0.5,120,190)
   CALL drawGridXY(5,4,0.5,270,190)
   SET TEXT HEIGHT 10
   PLOT TEXT, AT -30,-12:"|0>"
   PLOT TEXT, AT 120,-12:"|1>"
   PLOT TEXT, AT 270,-12:"|2>"
   PLOT TEXT, AT -30,138:"|3>"
   PLOT TEXT, AT 120,138:"|4>"
   PLOT TEXT, AT 270,138:"|5>"
END SUB
!
EXTERNAL SUB drawGridXY(ist,sc,zMag,xShift,yShift)
   DECLARE EXTERNAL SUB plotLines3D
   FOR j=0 TO NNy-1 STEP 1
      FOR i=0 TO NNx-1
         LET psi = sdState(ist,i,j,NNz/2)*200
         IF psi>1 THEN
            SET LINE COLOR 4 ! red
         ELSEIF psi<-1 THEN
            SET LINE COLOR 2 ! blue
         ELSE
            SET LINE COLOR 3 ! potential:green
         END IF
         CALL plotLines3D(i*sc,j*sc,zMag*(psi + 2*vv(i,j,NNz/2)),xShift,yShift) !(x,y,z)
      NEXT i
      PLOT LINES
   NEXT j
   FOR i=0 TO NNx-1 STEP 1
      FOR j=0 TO NNy-1
         LET psi = sdState(ist,i,j,NNz/2)*200
         IF psi>1 THEN
            SET LINE COLOR 4 ! red
         ELSEIF psi<-1 THEN
            SET LINE COLOR 2 ! blue
         ELSE
            SET LINE COLOR 3 ! potential:green
         END IF
         CALL plotLines3D(i*sc,j*sc,zMag*(psi + 2*vv(i,j,NNz/2)),xShift,yShift) !(x,y,z)
      NEXT j
      PLOT LINES
   NEXT i
END SUB
!
END MODULE
 		 

Re: 実行時内部エラーのご報告

 投稿者：SHIRAISHI Kazuo  投稿日：2017年 2月15日(水)17時31分51秒
  > No.4269[元記事へ]

mikeさんへのお返事です。

とりあえず手元にあった
OS 10.5.8
2Ghz Core 2 duo
1GB 667Mhz DDR2
のMAC BOOKだと
Iteration Count=2000
まで実行して普通に終了しました。
ただし、使用した十進BASICのバージョンは6.6.0です。
（計算に関係する部分は0.6.6.0から変わっていません）


> 下記のプログラムの実行中に内部エラーが発生しました。
> EXTERNAL FUNCTION getZminApex(sc,xShift,yShift)
> のデバグのためprintで変数の変化を見ながら実行中にエラーが発生しました。
> 立方体が１回転する前に発生しました。
> 実行環境は　十進BASIC 0.6.6.0 / macOS 10.7.5　/　2.7Ghz Intel Core i7
> メモリ 4GB 1333Mhz DDR3　です。
> ご検討ください。
>
> !
> ! ========= steepest descent method 3D ==========
> !
> ! eigenFunctionSD2D.bas
> ! Mitsuru Ikeuchi (C) Copyleft
> !
> ! ver 0.0.0   2017.01.31  created
> !
> OPTION ARITHMETIC NATIVE
> DECLARE EXTERNAL SUB sd3d.setInitialCondition,sd3d.SDiteration,sd3d.drawState
> DECLARE EXTERNAL FUNCTION sd3d.iterationCount,sd3d.stateEnergy,INKEY$
> LET Ax = -PI/24
> LET Ay = -PI/12
> LET drawMode = 0
> !setInitialCondition(stateMax,kx0,ky0)
> CALL setInitialCondition(10,2.0,2.0)
> !
> LET ist = 0
> FOR it=1 TO 1000
>    LET Ay=Ay +PI/180
>    LET dump = 0.05
>    CALL SDiteration(6, 2, dump) !(stateMax, iterMax, dump)
>    CALL drawState(ist,Ax,Ay,drawMode) !(state,xRotateAngle,yRotateAngle,drawMode)
>    LET S$=INKEY$
>    IF S$="0" THEN LET ist = MOD(ist+1,6)
>    IF S$="1" THEN LET drawMode = MOD(drawMode+1,3)
>    IF S$="." THEN EXIT FOR
> NEXT it
> !
> END
>
>
> EXTERNAL FUNCTION INKEY$
> OPTION ARITHMETIC NATIVE
> SET ECHO "OFF"
> LET S$=""
> CHARACTER INPUT NOWAIT: S$
> LET INKEY$=S$
> END FUNCTION
> !
> ! ---------- steepest descent method 3D ----------
> !
> ! system Hamiltonian: H = -delta/2 + V(r) , delta r = div grad r
> ! eigen energy Ei, eigen function set { |i> }
> !
> ! procedure : successive approximation
> ! (i) trial function set { |0>,|1>,..,|i>,.. }
> ! (2) energy of |i> : ei = <i|H|i>/<i|i>
> ! (3) steepest gradient direction (H-ei)|i>
> ! (4) next generation : |i(next)> = |i> - dumpingFactor*(H-ei)|i>
> ! (5) orthogonalization { |0>,|1>,..,|i>,.. }  (Gram-Schmidt)
> ! (6) goto (2)
> !
> MODULE sd3d
> MODULE OPTION ARITHMETIC NATIVE
> OPTION BASE 0
> PUBLIC SUB setInitialCondition, SDiteration, drawState
> PUBLIC FUNCTION iterationCount, stateEnergy
> SHARE NUMERIC NNx, NNy, NNz, dx, dy, dz, iterCount
> SHARE NUMERIC cosAx,sinAx,cosAy,sinAy,cx0,cy0,cz0 !--- 3D graphics
> SHARE NUMERIC sdEnergy(10)           ! electron state energy
> SHARE NUMERIC sdState(10,65,65,65)   ! electron states 0...20 0:ground state
> SHARE NUMERIC wrk(65,65,65)          ! state work space in steepestDescent
> SHARE NUMERIC vv(65,65,65)           ! external potential
> SHARE NUMERIC boxApex(0 TO 7,0 TO 2) ! boxApex(i,j) i-th box apex j-th coordinate x,y,z
> SHARE NUMERIC boxEdge(0 TO 11,0 TO 2)! boxEdge(i,j) i-th edge j=0,1:j-th apex, j=2:for marking
> LET NNx = 32                         ! x-max number of sdState(,NNx,NNy,NNz)
> LET NNy = 32                         ! y-max number of sdState(,NNx,NNy,NNz)
> LET NNz = 32                         ! z-max number of sdState(,NNx,NNy,NNz)
> LET dx = 0.5                         ! (au) x division
> LET dy = 0.5                         ! (au) y division
> LET dz = 0.5                         ! (au) y division
> LET iterCount = 0                    ! sd iteration count
> !
> ! ---------- set initial condition
> !
> EXTERNAL SUB setInitialCondition(stateMax,kx0,ky0)   !public
>    DECLARE EXTERNAL SUB setInitialState,setPotential,setBox
>    RANDOMIZE
>    CALL setInitialState(stateMax)
>    CALL setPotential(8,8,8)
>    CALL setBox
>    ! set window
>    LET xMargin = 60
>    LET yMargin = 120
>    SET WINDOW -xMargin,500-xMargin,-yMargin,500-yMargin
> END SUB
> !
> EXTERNAL SUB setInitialState(stateMax)
>    RANDOMIZE
>    FOR ist=0 TO stateMax-1
>       FOR i=1 TO NNx-2
>          FOR j=1 TO NNy-2
>             FOR k=1 TO NNz-2
>                LET sdState(ist,i,j,k) = RND-0.5
>             NEXT k
>          NEXT j
>       NEXT i
>       CALL normalizeState(ist)
>    NEXT ist
> END SUB
> !
> EXTERNAL SUB setPotential(x0,y0,z0)
>    FOR i=0 TO NNx-1
>       FOR j=0 TO NNy-1
>          FOR k=0 TO NNz-1
>             LET x = i*dx
>             LET y = j*dy
>             LET z = k*dz
>             LET r = SQR((x-x0)*(x-x0)+(y-y0)*(y-y0)+(z-z0)*(z-z0))
>             IF r<0.25 THEN LET r = 0.25
>             LET vv(i,j,k) = -1/r
>          NEXT k
>       NEXT j
>    NEXT i
> END SUB
>
> EXTERNAL SUB setPotential2(x0,y0,z0)
>    FOR i=0 TO NNx-1
>       FOR j=0 TO NNy-1
>          FOR k=0 TO NNz-1
>             LET x = i*dx
>             LET y = j*dy
>             LET z = k*dz
>             LET vv(i,j,k) = 2.0*((x-x0)*(x-x0)+(y-y0)*(y-y0)+(z-z0)*(z-z0))
>          NEXT k
>       NEXT j
>    NEXT i
> END SUB
>
> EXTERNAL SUB setBox
>    IF boxApex(7,2)<>1 THEN
>       DATA 0,0,0,  1,0,0,  0,1,0,  1,1,0,  0,0,1,  1,0,1,  0,1,1,  1,1,1
>       MAT READ boxApex !(0 TO 7,0 TO 2)
>       DATA 0,1,9,  0,2,9,  0,4,9,  1,3,9,  1,5,9,  2,3,9,  2,6,9,  3,7,9,   4,5,9,  4,6,9,  5,7,9,  6,7,9
>       MAT READ boxEdge !(0 TO 11,0 TO 2)
>    END IF
> END SUB
> !
> ! ---------- steepest descent iteration
> !
> EXTERNAL SUB SDiteration(stateMax, iterMax, dump) !public
>    DECLARE EXTERNAL FUNCTION steepestDescent
>    DECLARE EXTERNAL SUB GramSchmidt
>    FOR i=0 TO iterMax-1
>       FOR ist=0 TO stateMax-1
>          LET sdEnergy(ist) = steepestDescent(ist, dump)
>       NEXT ist
>       CALL GramSchmidt(stateMax)
>       LET iterCount = iterCount + 1
>    NEXT i
> END SUB
> !
> EXTERNAL FUNCTION steepestDescent(ist,dump)  !---  steepest descent method
>    DECLARE EXTERNAL FUNCTION energyOfState
>    DECLARE EXTERNAL SUB normalizeState
>    LET h2 = 2*dx*dx
>    LET ei = energyOfState(ist)               !---  E_ist = <ist|H|ist>
>    FOR i=1 TO NNx-2                          !---  |wrk> = (H-E_ist)|ist>
>       FOR j=1 TO NNy-2
>          FOR k=1 TO NNz-2
>             LET kesdState = (6*sdState(ist,i,j,k)-sdState(ist,i+1,j,k)-sdState(ist,i-1,j,k)&
> &-sdState(ist,i,j+1,k)-sdState(ist,i,j-1,k)-sdState(ist,i,j,k+1)-sdState(ist,i,j,k-1))/h2
>             LET wrk(i,j,k) = kesdState+(vv(i,j,k)-ei)*sdState(ist,i,j,k)
>          NEXT k
>       NEXT j
>    NEXT i
>    FOR i=1 TO NNx-2                          !---  |ist(next)> = |ist> - dump*|wrk>  ( norm(|ist(next)>) <>1 )
>       FOR j=1 TO NNy-2
>          FOR k=1 TO NNz-2
>             LET sdState(ist,i,j,k) = sdState(ist,i,j,k)-dump*wrk(i,j,k)
>          NEXT k
>       NEXT j
>    NEXT i
>    CALL normalizeState(ist)
>    LET steepestDescent = ei
> END FUNCTION
> !
> EXTERNAL FUNCTION energyOfState(ist) !--- E_ist = <ist|H|ist>/<ist|ist>
>    LET h2 = 2*dx*dx
>    LET s = 0
>    LET sn=0
>    FOR i=1 TO NNx-2
>       FOR j=1 TO NNy-2
>          FOR k=1 TO NNz-2
>             LET kesdState = (6*sdState(ist,i,j,k)-sdState(ist,i+1,j,k)-sdState(ist,i-1,j,k)&
> &-sdState(ist,i,j+1,k)-sdState(ist,i,j-1,k)-sdState(ist,i,j,k+1)-sdState(ist,i,j,k-1))/h2
>             LET s = s + sdState(ist,i,j,k)*(kesdState+vv(i,j,k)*sdState(ist,i,j,k))
>             LET sn = sn+sdState(ist,i,j,k)*sdState(ist,i,j,k)
>          NEXT k
>       NEXT j
>    next i
>    LET energyOfState = s/sn
> END FUNCTION
> !
> EXTERNAL SUB GramSchmidt(stateMax)  !--- Gram-Schmidt orthonormalization
>    DECLARE EXTERNAL FUNCTION innerProduct
>    DECLARE EXTERNAL SUB normalizeState
>    CALL normalizeState(0)
>    FOR istate=1 TO stateMax-1
>       FOR ist=0 TO istate-1
>          LET s = innerProduct(ist,istate)
>          FOR i=1 TO NNx-2
>             FOR j=1 TO NNy-2
>                FOR k=1 TO NNz-2
>                   LET sdState(istate,i,j,k) = sdState(istate,i,j,k) - s*sdState(ist,i,j,k)
>                NEXT k
>             NEXT j
>          NEXT i
>       NEXT ist
>       CALL normalizeState(istate)
>    NEXT iState
> END SUB
> !
> EXTERNAL FUNCTION innerProduct(ist,jst)  !--- <ist|jst>
>    LET s = 0
>    FOR i=1 TO NNx-2
>       FOR j=1 TO NNy-2
>          FOR k=1 TO NNz-2
>             LET s = s + sdState(ist,i,j,k)*sdState(jst,i,j,k)
>          NEXT k
>       NEXT j
>    NEXT i
>    LET innerProduct = s*dx*dy*dz
> END FUNCTION
> !
> EXTERNAL SUB normalizeState(ist)
>    LET s = 0
>    FOR i=1 TO NNx-2
>       FOR j=1 TO NNy-2
>          FOR k=1 TO NNz-2
>             LET s = s + sdState(ist,i,j,k)*sdState(ist,i,j,k)
>          NEXT k
>       NEXT j
>    NEXT i
>    LET a = SQR(1/(s*dx*dy*dz))
>    FOR i=1 TO NNx-2
>       FOR j=1 TO NNy-2
>          FOR k=1 TO NNz-2
>             LET sdState(ist,i,j,k) = a*sdState(ist,i,j,k)
>          NEXT k
>       NEXT j
>    NEXT i
> END SUB
> !
> ! ---------- utility
> !
> EXTERNAL FUNCTION iterationCount
>    LET iterationCount = iterCount
> END FUNCTION
> !
> EXTERNAL FUNCTION stateEnergy(ist)
>    LET stateEnergy = sdEnergy(ist)
> END FUNCTION
> !
> ! ---------- 3D graphics
> !
> EXTERNAL SUB setRotateXYParameters(angleX,angleY,xCenter,yCenter,zCenter)
>    LET cosAx = COS(angleX)
>    LET sinAx = SIN(angleX)
>    LET cosAy = COS(angleY)
>    LET sinAy = SIN(angleY)
>    LET cx0 = xCenter
>    LET cy0 = yCenter
>    LET cz0 = zCenter
> END SUB
>
> EXTERNAL SUB plotLines3D(x,y,z,xShift,yShift)
>    LET x1 = cosAy*(x-cx0)+sinAy*(sinAx*(y-cy0)+cosAx*(z-cz0))
>    LET y1 = cosAx*(y-cy0)-sinAx*(z-cz0)
>    LET z1 =-sinAy*(x-cx0)+cosAy*(sinAx*(y-cy0)+cosAx*(z-cz0))
>    PLOT LINES: x1+cx0+xShift, y1+cy0+yShift; !z=z1+cz0
> END SUB
>
> EXTERNAL SUB drawDisk3D(x,y,z,r,xShift,yShift)
>    LET x1 = cosAy*(x-cx0)+sinAy*(sinAx*(y-cy0)+cosAx*(z-cz0))
>    LET y1 = cosAx*(y-cy0)-sinAx*(z-cz0)
>    LET z1 =-sinAy*(x-cx0)+cosAy*(sinAx*(y-cy0)+cosAx*(z-cz0))
>    DRAW disk WITH SCALE(r)*SHIFT(x1+cx0+xShift,y1+cy0+yShift)
> END SUB
>
> EXTERNAL SUB plotText3DAt(x,y,z,xShift,yShift,A$)
>    LET x1 = cosAy*(x-cx0)+sinAy*(sinAx*(y-cy0)+cosAx*(z-cz0))
>    LET y1 = cosAx*(y-cy0)-sinAx*(z-cz0)
>    LET z1 =-sinAy*(x-cx0)+cosAy*(sinAx*(y-cy0)+cosAx*(z-cz0))
>    PLOT TEXT,AT x1+cx0+xShift,y1+cy0+yShift:A$
> END SUB
>
> EXTERNAL FUNCTION getRotateZ(x,y,z,xShift,yShift)
>    LET x1 = cosAy*(x-cx0)+sinAy*(sinAx*(y-cy0)+cosAx*(z-cz0))
>    LET y1 = cosAx*(y-cy0)-sinAx*(z-cz0)
>    LET z1 =-sinAy*(x-cx0)+cosAy*(sinAx*(y-cy0)+cosAx*(z-cz0))
>    LET getRotateZ = z1+cz0
> END FUNCTION
> !
> ! ---------- drawState
> !
> EXTERNAL SUB drawState(ist,xRotateAngle,yRotateAngle,drawMode)
> !DECLARE EXTERNAL SUB setRotateXYParameters,drawDensity3D,drawStateGrid,drawStateGridx6
>    SET DRAW MODE HIDDEN
>    CLEAR
>    LET sc = 6
>    IF drawMode=0 THEN
>       CALL drawDensity3D(ist,sc, xRotateAngle,yRotateAngle)
>    ELSEIF drawMode=1 THEN
>       CALL setRotateXYParameters(-PI/3,-PI/6,(NNx/2)*sc,(NNy/2)*sc,(NNy/2)*sc)
>       CALL drawStateGrid(ist,320/NNx,1.0,0,50) !(ist,sc,zMag,xShift,yShift)
>    ELSEIF drawMode=2 THEN
>       CALL setRotateXYParameters(-PI/3,-PI/6,(NNx/2)*sc,(NNy/2)*sc,(NNz/2)*sc)
>       CALL drawStateGridx6
>    ELSE
>    !
>    END IF
>    !--- caption
>    SET TEXT HEIGHT 10
>    SET TEXT COLOR 1 ! black
>    PLOT TEXT, AT 0,-50 ,USING "iterarion count =##### ":iterCount
>    PLOT TEXT, AT 0,-65 ,USING "E0 =###.######### E3 =###.#########(au)":sdEnergy(0),sdEnergy(3)
>    PLOT TEXT, AT 0,-80 ,USING "E1 =###.######### E4 =###.#########(au)":sdEnergy(1),sdEnergy(4)
>    PLOT TEXT, AT 0,-95 ,USING "E2 =###.######### E5 =###.#########(au)":sdEnergy(2),sdEnergy(5)
>    PLOT TEXT, AT 0,-110 :"steepest descent method 2D"
>    SET DRAW MODE EXPLICIT
> END SUB
> !
> EXTERNAL SUB drawDensity3D(ist,sc,xRotateAngle,yRotateAngle) !----- drawMode=0
>    DECLARE EXTERNAL SUB setRotateXYParameters, drawDisk3D, plotEdge3D
>    DECLARE EXTERNAL FUNCTION getRotateZ,getZminApex
>    CALL setRotateXYParameters(xRotateAngle,yRotateAngle,NNx*sc/2,NNy*sc/2,NNx*sc/2)
>    LET xOffset = 80
>    LET yOffset = 50
>    SET LINE COLOR 1
>    FOR i=1 TO NNx-2
>       FOR j=1 TO NNy-2
>          FOR k=1 TO NNz-2
>             LET psi = sdState(ist,i,j,k)
>             LET psi2 = psi*psi
>             IF psi2>0.001 THEN
>                IF psi>0 THEN SET AREA COLOR 4 ELSE SET AREA COLOR 2
>                CALL drawDisk3D(i*sc,j*sc,k*sc,ABS(psi)*10,xOffset,yOffset)
>             END IF
>          NEXT k
>       NEXT j
>    NEXT i
>    CALL plotEdge3D(0,sc,xOffset,yOffset,1) !x-axis:black !(iEdge,xShift,yShift,lineColorIndex)
>    CALL plotEdge3D(1,sc,xOffset,yOffset,2) !y-axis:blue
>    CALL plotEdge3D(2,sc,xOffset,yOffset,3) !z-axis:green
>    FOR iEdge=3 TO 11
>       CALL plotEdge3D(iEdge,sc,xOffset,yOffset,8) !other edge:grey
>    NEXT iEdge
>    CALL plotText3DAt(0,0,NNz*sc,xOffset,yOffset,"Z"&STR$(getRotateZ(0,0,NNz*sc,xOffset,yOffset)))
>    CALL plotText3DAt(0,0,0,xOffset,yOffset,"Z"&STR$(getRotateZ(0,0,0,xOffset,yOffset)))
>    LET apx = getZminApex(sc,xOffset,yOffset)
>    CALL drawDisk3D(boxApex(apx,0)*NNx*sc,boxApex(apx,1)*NNy*sc,boxApex(apx,2)*NNz*sc,10,xOffset,yOffset)
>    SET TEXT COLOR 1
>    SET TEXT HEIGHT 10
>    PLOT TEXT, AT 0,-20:"|"&STR$(ist)&">"
>    PLOT TEXT, AT 0,-35 ,USING "Ax =###.#(deg) Ay =###.#(deg)":MOD(xRotateAngle*180/PI,360),MOD(yRotateAngle*180/PI,360)
> END SUB
> !
> EXTERNAL SUB plotEdge3D(iEdge,sc,xShift,yShift,lineColorIndex)
>    DECLARE EXTERNAL SUB plotLines3D
>    SET LINE COLOR lineColorIndex
>    FOR j=0 TO 1
>       LET apex = boxEdge(iEdge,j)
>       LET x = boxApex(apex,0)*NNx*sc
>       LET y = boxApex(apex,1)*NNy*sc
>       LET z = boxApex(apex,2)*NNz*sc
>       CALL plotLines3D(x,y,z,xShift,yShift)
>    NEXT j
>    PLOT LINES
> END SUB
>
> EXTERNAL FUNCTION getZminApex(sc,xShift,yShift)
>    DECLARE EXTERNAL FUNCTION getRotateZ !(x,y,z,xShift,yShift)
>    LET zmin = 1E99
>    FOR iApex = 0 TO 7
>       LET x = boxApex(iApex,0)*NNx*sc
>       LET y = boxApex(iApex,1)*NNy*sc
>       LET z = boxApex(iApex,2)*NNz*sc
>       LET rotz = getRotateZ(x,y,z,xShift,yShift)
>       PRINT USING "###   rotZ=####.##   x=###.##  y=###.##   z=###.##":iApex,rotz,x,y,z
>       IF rotz<zmin THEN
>          LET zmin = rotz
>          LET zminApex = iApex
>       END IF
>    NEXT iApex
>    LET getZminApex = zminApex
> END FUNCTION
> !
> EXTERNAL SUB drawStateGrid(ist,sc,zMag,xShift,yShift)  !--- drawMode=1
>    DECLARE EXTERNAL SUB drawGridXY
>    CALL drawGridXY(ist,sc,zMag,xShift,yShift)
>    SET TEXT COLOR 1
>    SET TEXT HEIGHT 10
>    PLOT TEXT, AT 0,-20:"|"&STR$(ist)&">"
> END SUB
> !
> EXTERNAL SUB drawStateGridx6
>    DECLARE EXTERNAL SUB drawGridXY
>    CALL drawGridXY(0,4,0.5,-30,40) !(ist,sc,zMag,xShift,yShift)  !--- drawMode=2
>    CALL drawGridXY(1,4,0.5,120,40)
>    CALL drawGridXY(2,4,0.5,270,40)
>    CALL drawGridXY(3,4,0.5,-30,190)
>    CALL drawGridXY(4,4,0.5,120,190)
>    CALL drawGridXY(5,4,0.5,270,190)
>    SET TEXT HEIGHT 10
>    PLOT TEXT, AT -30,-12:"|0>"
>    PLOT TEXT, AT 120,-12:"|1>"
>    PLOT TEXT, AT 270,-12:"|2>"
>    PLOT TEXT, AT -30,138:"|3>"
>    PLOT TEXT, AT 120,138:"|4>"
>    PLOT TEXT, AT 270,138:"|5>"
> END SUB
> !
> EXTERNAL SUB drawGridXY(ist,sc,zMag,xShift,yShift)
>    DECLARE EXTERNAL SUB plotLines3D
>    FOR j=0 TO NNy-1 STEP 1
>       FOR i=0 TO NNx-1
>          LET psi = sdState(ist,i,j,NNz/2)*200
>          IF psi>1 THEN
>             SET LINE COLOR 4 ! red
>          ELSEIF psi<-1 THEN
>             SET LINE COLOR 2 ! blue
>          ELSE
>             SET LINE COLOR 3 ! potential:green
>          END IF
>          CALL plotLines3D(i*sc,j*sc,zMag*(psi + 2*vv(i,j,NNz/2)),xShift,yShift) !(x,y,z)
>       NEXT i
>       PLOT LINES
>    NEXT j
>    FOR i=0 TO NNx-1 STEP 1
>       FOR j=0 TO NNy-1
>          LET psi = sdState(ist,i,j,NNz/2)*200
>          IF psi>1 THEN
>             SET LINE COLOR 4 ! red
>          ELSEIF psi<-1 THEN
>             SET LINE COLOR 2 ! blue
>          ELSE
>             SET LINE COLOR 3 ! potential:green
>          END IF
>          CALL plotLines3D(i*sc,j*sc,zMag*(psi + 2*vv(i,j,NNz/2)),xShift,yShift) !(x,y,z)
>       NEXT j
>       PLOT LINES
>    NEXT i
> END SUB
> !
> END MODULE
> 		 

Re: 実行時内部エラーのご報告

 投稿者：白石　和夫  投稿日：2017年 2月15日(水)17時51分56秒
  > No.4270[元記事へ]

SHIRAISHI Kazuoさんへのお返事です。

後日，OS 10.9で試してみます。
 		 

Re: 実行時内部エラーのご報告

 投稿者：mike  投稿日：2017年 2月15日(水)19時18分9秒
  白石　和夫さんへのお返事です。

白石先生、ご検討いただき、ありがとうございます。

BASIC処理系は終了していたのですが、再度立ち上げ、当プログラムを走らせたところ、
内部エラーが再現しました。（デバッグ前のプログラムは正常終了します。）

しかし、内部エラーが起きてから、macOSはそのまま使っていたことに気付き、
再起動後、当プログラムを走らせたところ、内部エラーは起こらず、
そのまま2000回までで正常終了しました。
再度。再起動し試みましたが、同様に正常終了しました。

macが一時的に異常であった可能性が高そうです。
お騒がわせして申し訳ありませんでした。
 		 

戻る