[013]定常状態の電子状態を求めるプログラム（3次元）投稿者：mike 投稿日：2017年 4月 7日(金)09時39分20秒 |
|
定常状態の電子のエネルギーと電子状態を求める３次元のプログラム(013harmonicsSD3D.bas)を公開します。 定常状態の電子のエネルギーと電子状態を求める問題は、系のハミルトニアンをHとして H |i> = e_i |i> の最小固有値とその近くの低いエネルギー固有値と固有関数を求める問題になります。 最急降下法（steepest descent method）によって、この最小固有値付近の固有値と固有関数を求めます。 適当に与えた試行関数（ここでは乱数NNx*NNy*NNz個の配列）の組{|i>}から出発し、エネルギーei=<i|H|i>を求め、 エネルギーの減少する方向 (H-ei)|i> から 次世代の関数|i(next)> = |i> - damp(H-ei)|i>を求め、 関数の組を直交化することを繰り返す（逐次近似する）ことで小さい方からN個の固有値と固有関数を求めます。 本プログラムでは、 水素様(V=-1/r)、放物型(V=0.5*r^2)、井戸型(V=0(r<4) =16)のポテンシャルを選択できます。 表示の説明： 波動関数は赤色が正、青色が負の値になることを表します。緑色はポテンシャルです。 ポテンシャルが放物型のとき、理論的には固有エネルギーは3次元の場合、E0=1.5, E1=E2=E3=2.5,.. となります。 |1>,|2>,|3>はほぼ同じ固有値になる（縮退している）ため、向きが不定です。 問題点： ３Dグラフィックの表示に問題があります。３Dの波動関数を表すのに|i>の空間の各格子点における大きさを その点の値に比例した丸の大きさで表すのですが、点がたくさんあるため、奥行き方向にソートして、 奥から描いていくということをさぼっています。このため、前後関係が逆になる場合、変な表示になります。 格子点の各点の位置関係は規則的なため、これを利用して矛盾のない表示ができないかと試行していますが、 今のところ成功していません。何かいいアイデアがありましたら教えてください。 試験環境： 本プログラムは十進BASIC 6.6.2 / macOS 10.7.5, 十進BASIC Ver 7.8.0 / windows 10 でテストしました。 ---------------- ! ! ========= Harmonics - steepest descent method 3D ========== ! ! 013harmonicsSD3D.bas ! Mitsuru Ikeuchi (C) copyleft ! ! ver 0.0.1 2017.04.06 created ! OPTION ARITHMETIC NATIVE DECLARE EXTERNAL SUB sd3d.setInitialCondition,sd3d.SDiteration,sd3d.drawState DECLARE EXTERNAL FUNCTION INKEY$ LET Ax = PI/12 !rotate angle around x-axis LET Ay = -PI/6 !rotate angle around y-axis LET ddAy = PI/180 !Ay=Ay+ddAy LET vIndex = 1 !vIndex=0: V=1/r, 1: V=(1/2)r^2, 2:V=0(r<=4) or V=16(r>4) LET stateMax = 6 !state 0,1,...,stateMax-1 LET drawMode = 0 !drawMode=0:density3D, 1:grid |i> in (x,y,0) LET dispState = 0 !display state |i> CALL setInitialCondition(stateMax,vIndex) DO LET Ay = Ay + ddAy IF Ay>PI THEN LET Ay = Ay - 2*PI IF Ay<-PI THEN LET Ay = Ay + 2*PI CALL SDiteration(stateMax, 2) !2 - iteration in SDiteration() CALL drawState(dispState,Ax,Ay,drawMode,vIndex) LET S$=INKEY$ IF S$="0" OR S$="1" OR S$="2" OR S$="3" OR S$="4" OR S$="5" THEN LET dispState = VAL(S$) IF S$="W" OR S$="w" THEN LET ddAy = -PI/180 IF S$="E" OR S$="e" THEN LET ddAy = PI/180 IF S$="S" OR S$="s" THEN LET ddAy = 0 IF s$="/" THEN LET drawMode = MOD(drawMode+1,2) IF S$="7" OR S$="8" OR S$="9" THEN LET vIndex = VAL(S$)-7 CALL setInitialCondition(stateMax,vIndex) END IF IF S$="." THEN EXIT DO LOOP END EXTERNAL FUNCTION INKEY$ !--- from decimal BASIC library inkey$.bas OPTION ARITHMETIC NATIVE SET ECHO "OFF" LET S$="" CHARACTER INPUT NOWAIT: S$ LET INKEY$=S$ END FUNCTION ! ========== SD3D(steepest descent method 3D) module ========== ! ! steepest descent method ! system Hamiltonian: H = -delta/2 + V(r) , delta r = div grad r ! eigen energy set { 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> - dampingFactor*(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 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 xApex(0 TO 7),yApex(0 TO 7),zApex(0 TO 7) ! boxApex x- y- z-coordinate SHARE NUMERIC pxApex(0 TO 7),pyApex(0 TO 7),pzApex(0 TO 7) ! rotated boxApex x- y- z-coordinate 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,vIndex) !public DECLARE EXTERNAL SUB setInitialState,setPotential,setBox RANDOMIZE LET iterCount = 0 CALL setInitialState(stateMax) CALL setPotential(vIndex) CALL setBox ! set window SET WINDOW 0,500,0,500 END SUB EXTERNAL SUB setInitialState(stateMax) DECLARE EXTERNAL SUB normalizeState 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(vIndex) 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-8)*(x-8)+(y-8)*(y-8)+(z-8)*(z-8)) IF vIndex=0 THEN IF r<0.25 THEN LET r = 0.25 LET vv(i,j,k) = -1/r ELSEIF vIndex=1 THEN LET vv(i,j,k) = MIN(0.5*r*r,18) ELSEIF vIndex=2 THEN IF r<=4 THEN LET vv(i,j,k) = 0 ELSE LET vv(i,j,k) = 16 ELSE LET vv(i,j,k) = 0 END IF NEXT k NEXT j NEXT i END SUB EXTERNAL SUB setBox IF boxEdge(11,1)<>7 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 FOR i=0 TO 7 READ x,y,z LET xApex(i) = x*NNx*dx LET yApex(i) = y*NNy*dy LET zApex(i) = z*NNz*dz NEXT i 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) !public DECLARE EXTERNAL FUNCTION steepestDescent DECLARE EXTERNAL SUB GramSchmidt,sortState LET damp = 0.05 !damping factor in steepest descent FOR i=0 TO iterMax-1 FOR ist=0 TO stateMax-1 LET sdEnergy(ist) = steepestDescent(ist, damp) NEXT ist CALL GramSchmidt(stateMax) CALL sortState(stateMax) LET iterCount = iterCount + 1 NEXT i END SUB EXTERNAL FUNCTION steepestDescent(ist, damp) !--- 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> - damp*|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)-damp*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 SUB sortState(stateMax) FOR ist=stateMax-2 TO 0 STEP -1 IF sdEnergy(ist)>sdEnergy(ist+1)+0.00001 THEN FOR i=0 TO NNx-1 FOR j=0 TO NNy-1 FOR k=0 TO NNz-1 LET w = sdState(ist,i,j,k) LET sdState(ist,i,j,k) = sdState(ist+1,i,j,k) LET sdState(ist+1,i,j,k) = w NEXT k NEXT j NEXT i LET w = sdEnergy(ist) LET sdEnergy(ist) = sdEnergy(ist+1) LET sdEnergy(ist+1) = w END IF NEXT ist 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 ! ---------- 3D graphics aid EXTERNAL SUB setRotateXY(angleX,angleY) LET cosAx = COS(angleX) LET sinAx = SIN(angleX) LET cosAy = COS(angleY) LET sinAy = SIN(angleY) LET cx0 = 0.5*NNx*dx LET cy0 = 0.5*NNy*dy LET cz0 = 0.5*NNz*dz END SUB EXTERNAL SUB rotateXY !particles and box apex FOR i=0 TO 7 LET pxApex(i) = cosAy*(xApex(i)-cx0)+sinAy*(sinAx*(yApex(i)-cy0)+cosAx*(zApex(i)-cz0)) + cx0 LET pyApex(i) = cosAx*(yApex(i)-cy0)-sinAx*(zApex(i)-cz0) + cy0 LET pzApex(i) =-sinAy*(xApex(i)-cx0)+cosAy*(sinAx*(yApex(i)-cy0)+cosAx*(zApex(i)-cz0)) + cz0 NEXT i END SUB EXTERNAL SUB drawDisk3D(x,y,z,r,mag,xShift,yShift) LET x1 = cosAy*(x-cx0)+sinAy*(sinAx*(y-cy0)+cosAx*(z-cz0)) + cx0 LET y1 = cosAx*(y-cy0)-sinAx*(z-cz0) + cy0 LET z1 =-sinAy*(x-cx0)+cosAy*(sinAx*(y-cy0)+cosAx*(z-cz0)) + cz0 DRAW disk WITH SCALE(r)*SHIFT(x1*mag+xShift,y1*mag+yShift) END SUB EXTERNAL SUB plotLines3D(x,y,z,mag,xShift,yShift) LET x1 = cosAy*(x-cx0)+sinAy*(sinAx*(y-cy0)+cosAx*(z-cz0)) + cx0 LET y1 = cosAx*(y-cy0)-sinAx*(z-cz0) + cy0 LET z1 =-sinAy*(x-cx0)+cosAy*(sinAx*(y-cy0)+cosAx*(z-cz0)) + cz0 PLOT LINES: x1*mag+xShift, y1*mag+yShift; END SUB EXTERNAL SUB markFarEdge ! seek far apex --> iMin LET zMin = pzApex(0) LET iMin = 0 FOR i=1 TO 7 IF zMin>pzApex(i) THEN LET zMin = pzApex(i) LET iMin = i END IF NEXT i !mark far edge FOR iEdge = 0 TO 11 LET boxEdge(iEdge,2) = 0 IF (boxEdge(iEdge,0)=iMin OR boxEdge(iEdge,1)=iMin) THEN LET boxEdge(iEdge,2) = 1 NEXT iEdge END SUB EXTERNAL SUB plotNearEdge(mag,xp,yp,lineColor) DECLARE EXTERNAL SUB plotEdge FOR iEdge=0 TO 11 IF boxEdge(iEdge,2)=0 THEN !far edge mark = 0 CALL plotEdge(iEdge,mag,xp,yp,lineColor) END IF NEXT iEdge END SUB EXTERNAL SUB plotFarEdge(mag,xp,yp,lineColor) DECLARE EXTERNAL SUB plotEdge FOR iEdge=0 TO 11 IF boxEdge(iEdge,2)=1 THEN !far edge mark = 1 CALL plotEdge(iEdge,mag,xp,yp,lineColor) END IF NEXT iEdge END SUB EXTERNAL SUB plotEdge(iEdge,mag,xp,yp,lineColor) SET LINE COLOR lineColor FOR i=0 TO 1 LET iApex = boxEdge(iEdge,i) PLOT LINES: pxApex(iApex)*mag+xp, pyApex(iApex)*mag+yp; NEXT i PLOT LINES END SUB ! ---------- draw state EXTERNAL SUB drawState(ist,xRotateAngle,yRotateAngle,drawMode,vIndex) DECLARE EXTERNAL SUB drawDensity3D,setRotateXY,drawDensity3D,drawStateGrid SET DRAW MODE HIDDEN LET a$="H like harmonic quantumWell " CLEAR IF drawMode=0 THEN CALL drawDensity3D(ist,xRotateAngle,yRotateAngle,200/(NNx*dx),140,150) !!(ist,Ax,Ay,sc,xp,yp) ELSEIF drawMode=1 THEN CALL drawStateGrid(ist,200/(NNx*dx),0.1,140,200) !(ist,sc,zMag,xp,yp) END IF !--- caption SET TEXT HEIGHT 10 SET TEXT COLOR 1 ! black PLOT TEXT, AT 50, 70 ,USING "iterarion count =##### ":iterCount PLOT TEXT, AT 50, 55 ,USING "E0 =###.######### E3 =###.#########(au)":sdEnergy(0),sdEnergy(3) PLOT TEXT, AT 50, 40 ,USING "E1 =###.######### E4 =###.#########(au)":sdEnergy(1),sdEnergy(4) PLOT TEXT, AT 50, 25 ,USING "E2 =###.######### E5 =###.#########(au)":sdEnergy(2),sdEnergy(5) PLOT TEXT, AT 50, 10 :"steepest descent method 3D" PLOT TEXT, AT 50,470 :"[.] exit [/] change drawMode " PLOT TEXT, AT 50,455 :"[7] V=-1/r, [8] V=0.5r^2, [9] V=0(r<=4) =16(r>4)" PLOT TEXT, AT 50,440 :"[0],[1],...,[5] view state " PLOT TEXT, AT 50,425 :"[W] <<- rotate [S]top rotate ->> [E] " PLOT TEXT, AT 50,410 :"potential = "&a$(vIndex*12+1:vIndex*12+11) SET DRAW MODE EXPLICIT END SUB EXTERNAL SUB drawDensity3D(ist,xRotateAngle,yRotateAngle,sc,xp,yp) !----- drawMode=0 DECLARE EXTERNAL SUB setRotateXY,rotateXY,markFarEdge,plotFarEdge,drawDisk3D,plotNearEdge CALL setRotateXY(xRotateAngle,yRotateAngle) CALL rotateXY !rotateXY BOX CALL markFarEdge ! boxEdge(iEdge,2)=1:far side edge or 0:near side edge CALL plotFarEdge(sc,xp,yp,8) !8:gray 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*dx,j*dy,k*dz,ABS(psi)*10,sc,xp,yp) END IF NEXT k NEXT j NEXT i !SET AREA COLOR 1 !CALL drawDisk3D(0,0,0,5,sc,xp,yp) CALL plotNearEdge(sc,xp,yp,1) !1:black SET TEXT COLOR 1 SET TEXT HEIGHT 10 PLOT TEXT, AT 50,100:"|"&STR$(ist)&">" PLOT TEXT, AT 50, 85 ,USING "Ax =####.#(deg) Ay =####.#(deg)":xRotateAngle*180/PI,yRotateAngle*180/PI END SUB EXTERNAL SUB drawStateGrid(ist,sc,zMag,xp,yp) !--- drawMode=1 DECLARE EXTERNAL SUB setRotateXY,drawGridXY CALL setRotateXY(-PI/3,-PI/6) CALL drawGridXY(ist,sc,zMag,xp,yp) SET TEXT COLOR 1 SET TEXT HEIGHT 10 PLOT TEXT, AT 50,100:"|"&STR$(ist)&"> in (x,y,0)" 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*dx,j*dx,zMag*(psi + 2*vv(i,j,NNz/2)),sc,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*dx,j*dx,zMag*(psi + 2*vv(i,j,NNz/2)),sc,xShift,yShift) !(x,y,z) NEXT j PLOT LINES NEXT i END SUB END MODULE 
|