|
|
下記のプログラムの実行中に内部エラーが発生しました。
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
|
|