teacup. [ 掲示板 ] [ 掲示板作成 ] [ 有料掲示板 ] [ ブログ ]

 投稿者
  題名
  内容 入力補助画像・ファイル<IMG>タグが利用可能です。(詳細)
    
 URL
[ ケータイで使う ] [ BBSティッカー ] [ 書込み通知 ] [ 検索 ]

スレッド一覧

  1. スレッドが使えます(2)
  2. Paract BASIC(20)
  3. Amusement_Program(10)
  4. 改修予定のJIS非互換(3)
  5. 複数ページ長編プログラム(新規投稿)(16)
  6. 「十進BASIC第2掲示板」投稿記事リスト(17)
  7. Full BASIC互換ライブラリ(8)
  8. 「十進BASIC掲示板過去ログ」インデックス(トピック)(17)
  9. 人の色覚の数理(14)
  10. 「十進BASIC掲示板過去ログ」インデックス(ツリー)(91)
スレッド一覧(全10)  他のスレッドを探す  スレッド作成

*掲示板をお持ちでない方へ、まずは掲示板を作成しましょう。無料掲示板作成


Re: 変更の仕方

 投稿者:GAI  投稿日:2017年 9月29日(金)23時31分12秒
返信・引用
  > No.4402[元記事へ]

しばっちさんへのお返事です。


>
> 総数 2627625となりましたが、個数はあっているでしょうか?
>

プログラムありがとうございます。
計算上
16C4*12C4*8C4/4!=2627625
ですので正しく動いています。
 
 

Re: 変更の仕方

 投稿者:しばっち  投稿日:2017年 9月29日(金)22時36分12秒
返信・引用
  > No.4401[元記事へ]

GAIさんへのお返事です。

> これを異なる16個を4つずつ4組作るものに改変しようと試みていましたがいまいちやり方がわからなくなりました。
> どなたか変更するとどの様に変えればいいのか教えて下さい。
>

総数 2627625となりましたが、個数はあっているでしょうか?

FOR I1=2 TO 14
   FOR J1=I1+1 TO 15
      FOR K1=J1+1 TO 16
         LET W$="123456789ABCDEFG"
         LET A$=W$(1:1)&W$(I1:I1)&W$(J1:J1)&W$(K1:K1)
         LET W$(1:1)="*"
         LET W$(I1:I1)="*"
         LET W$(J1:J1)="*"
         LET W$(K1:K1)="*"
         CALL ERASE(W$)
         FOR I2=2 TO 10
            FOR J2=I2+1 TO 11
               FOR K2=J2+1 TO 12
                  LET X$=W$
                  LET A$=A$(1:4)&X$(1:1)&X$(I2:I2)&X$(J2:J2)&X$(K2:K2)
                  LET X$(1:1)="*"
                  LET X$(I2:I2)="*"
                  LET X$(J2:J2)="*"
                  LET X$(K2:K2)="*"
                  CALL ERASE(X$)
                  FOR I3=2 TO 6
                     FOR J3=I3+1 TO 7
                        FOR K3=J3+1 TO 8
                           LET Y$=X$
                           LET A$=A$(1:8)&Y$(1:1)&Y$(I3:I3)&Y$(J3:J3)&Y$(K3:K3)
                           LET Y$(1:1)="*"
                           LET Y$(I3:I3)="*"
                           LET Y$(J3:J3)="*"
                           LET Y$(K3:K3)="*"
                           CALL ERASE(Y$)
                           LET A$=A$&Y$
                           LET C=C+1
                           PRINT C;":";A$(1:4);"  ";A$(5:8);"  ";A$(9:12);"  ";A$(13:16) !' 2627625 まで
                           LET A$=A$(1:12)
                        NEXT K3
                     NEXT J3
                  NEXT I3
               NEXT K2
            NEXT J2
         NEXT I2
      NEXT K1
   NEXT J1
NEXT I1
END

EXTERNAL  SUB ERASE(A$)
DO
   LET FL=0
   FOR I=1 TO LEN(A$)
      IF A$(I:I)="*" THEN
         LET A$(I:I)=""
         LET FL=1
      END IF
   NEXT I
LOOP UNTIL FL=0
END SUB
 

変更の仕方

 投稿者:GAI  投稿日:2017年 9月29日(金)16時12分43秒
返信・引用
  以前異なる8個のものを2つずつ4組構成するプログラムで

DIM A$(105) !全パターン 例. 12345678なら、(1,2)(3,4)(5,6)(7,8)と読む
LET C=1
FOR i=2 TO 8 !1組目(1,i) ※comb(8,2)/4=7通り
   LET w$="12345678" !restore it

   LET A$(C)=w$(1:1)&w$(i:i) !その2つを抜き取って、連番を再構成する
   LET w$(i:i)=""
   LET w$(1:1)=""

   FOR j=2 TO 6 !2組目(1,j) ※comb(6,2)/3=5通り
      LET x$=w$ !restore it

      LET A$(C)=A$(C)(1:2)&x$(1:1)&x$(j:j)
      LET x$(j:j)=""
      LET x$(1:1)=""

      FOR k=2 TO 4 !3組目(1,k) ※comb(4,2)/2=3通り
         LET y$=x$ !restore it

         LET A$(C)=A$(C)(1:4)&y$(1:1)&y$(k:k)
         LET y$(k:k)=""
         LET y$(1:1)=""

         LET A$(C)=A$(C)&y$ !4組目 ※comb(2,2)/1=1通り


         PRINT C;": ";A$(C) !結果を表示する

         LET C=C+1
         IF C<=105 THEN LET A$(C)=A$(C-1)(1:4) !copy it
      NEXT k
   NEXT j
NEXT i
END

が紹介されていました。

これを異なる16個を4つずつ4組作るものに改変しようと試みていましたがいまいちやり方がわからなくなりました。
どなたか変更するとどの様に変えればいいのか教えて下さい。

 

[024]-続き

 投稿者:mike  投稿日:2017年 9月 6日(水)12時22分45秒
返信・引用
  プログラムの後半部分です

---------------------

! ---------- 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) !shift*xRotateAx*yRotateAy*(shift^-1)
   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

! ---------- drawState

EXTERNAL SUB drawState(drawMode)
   DECLARE EXTERNAL SUB setRotateXYParameters,drawState3D,draw3DPsix6,drawPsix6,drawRho
   DECLARE EXTERNAL FUNCTION totalEnergy
   SET DRAW MODE HIDDEN
   CALL setRotateXYParameters(-PI/6,-PI/12,NNx/2,NNy/2,0)
   CLEAR
   SET TEXT HEIGHT 8
   SET TEXT COLOR 1 ! black
   IF drawMode=0 THEN
      CALL drawRho(4,1.0,100,150)  !(sc,zMag,xShift,yShift)
      PLOT TEXT, AT 100, 150 :"electron density rho(x,y) and Veff(x,y)"
   ELSEIF drawMode=1 THEN
      CALL draw3DPsix6
   ELSEIF drawMode=2 THEN
      CALL drawPsix6
   END IF
   !--- caption
   PLOT TEXT, AT 50,115 ,USING "iteration count =######        error =##.##########":iterCount,iterationError
   PLOT TEXT, AT 50,100 ,USING "total energy =###.##########   N of electron =##":totalEnergy,numberOfElectron
   PLOT TEXT, AT 50, 85 ,USING "E0 =###.########## Occ =#.#####":sdEnergy(0),occupation(0)
   PLOT TEXT, AT 50, 70 ,USING "E1 =###.########## Occ =#.#####":sdEnergy(1),occupation(1)
   PLOT TEXT, AT 50, 55 ,USING "E2 =###.########## Occ =#.#####":sdEnergy(2),occupation(2)
   PLOT TEXT, AT 50, 40 ,USING "E3 =###.########## Occ =#.#####":sdEnergy(3),occupation(3)
   PLOT TEXT, AT 50, 25 ,USING "E4 =###.########## Occ =#.#####":sdEnergy(4),occupation(4)
   PLOT TEXT, AT 50, 10 :"RS-DFT - Local Density Approximation 2D"
   SET TEXT HEIGHT 10
   PLOT TEXT, AT 50,470 :"[esc] exit            [D] chage dispMode"
   PLOT TEXT, AT 50,455 :"[1] ... [6] number of electron"
   PLOT TEXT, AT 50,440 :"[7] hermonics k*x^2   [8] quantum well"
   SET DRAW MODE EXPLICIT
END SUB

EXTERNAL SUB draw3DPsix6
   DECLARE EXTERNAL SUB drawState3D
   CALL drawState3D(0,2,0.5, 30,120) !(ist,sc,zMag,xShift,yShift)
   CALL drawState3D(1,2,0.5,180,120)
   CALL drawState3D(2,2,0.5,330,120)
   CALL drawState3D(3,2,0.5, 30,270)
   CALL drawState3D(4,2,0.5,180,270)
   CALL drawState3D(5,2,0.5,330,270)
   SET TEXT HEIGHT 10
   PLOT TEXT, AT  30,140:"|0>"
   PLOT TEXT, AT 180,140:"|1>"
   PLOT TEXT, AT 330,140:"|2>"
   PLOT TEXT, AT  30,290:"|3>"
   PLOT TEXT, AT 180,290:"|4>"
   PLOT TEXT, AT 330,290:"|5>"
END SUB

EXTERNAL SUB drawRho(sc,zMag,xShift,yShift)
   DECLARE EXTERNAL SUB plotLines3D
   FOR j=0 TO NNy-1 STEP 1
      FOR i=0 TO NNx-1
         LET psi = rho(i,j)*2000
         IF psi>1 THEN
            SET LINE COLOR 13
         ELSE
            SET LINE COLOR 3 ! potential:green
         END IF
         CALL plotLines3D(i*sc,j*sc,zMag*(psi + vv(i,j)),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 = rho(i,j)*2000
         IF psi>1 THEN
            SET LINE COLOR 13
         ELSE
            SET LINE COLOR 3 ! potential:green
         END IF
         CALL plotLines3D(i*sc,j*sc,zMag*(psi + vv(i,j)),xShift,yShift) !(x,y,z)
      NEXT j
      PLOT LINES
   NEXT i
END SUB

EXTERNAL SUB drawState3D(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)*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 + vv(i,j)),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)*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 + vv(i,j)),xShift,yShift) !(x,y,z)
      NEXT j
      PLOT LINES
   NEXT i
END SUB

EXTERNAL SUB drawPsix6
   DECLARE EXTERNAL SUB drawPsi
   CALL drawPsi(0, 30,135) !(ist,xPos,yPos)
   CALL drawPsi(1,180,135)
   CALL drawPsi(2,330,135)
   CALL drawPsi(3, 30,285)
   CALL drawPsi(4,180,285)
   CALL drawPsi(5,330,285)
   SET TEXT HEIGHT 10
   PLOT TEXT, AT  30,140:"|0>"
   PLOT TEXT, AT 180,140:"|1>"
   PLOT TEXT, AT 330,140:"|2>"
   PLOT TEXT, AT  30,290:"|3>"
   PLOT TEXT, AT 180,290:"|4>"
   PLOT TEXT, AT 330,290:"|5>"
END SUB

EXTERNAL SUB drawPsi(ist,xPos,yPos)  !drawMode=0
   MAT psicol = vcol !--- psicol(,) <- Vcol(,) setted SUB initDraw
   FOR i=0 TO NNx-1 !--- set psicol(,) for MAT PLOT CELLS
      FOR j=0 TO NNy-1
         LET psi = sdState(ist,i,j)
         IF abs(psi)>0.002 THEN
            IF psi>=0 THEN
               LET col = psi*10
               IF col>1 THEN LET col = 1
               IF col<0 THEN LET col = 0
               LET psicol(i,j) = 40+INT(col*50)
            ELSE
               LET col = -psi*10
               IF col>1 THEN LET col = 1
               IF col<0 THEN LET col = 0
               LET psicol(i,j) = 100+INT(col*50)
            END IF
         END IF
      NEXT j
   NEXT i
   MAT PLOT CELLS,IN xPos,yPos; xPos+128,yPos+128 :psicol
END SUB

END MODULE
 

[024]多電子系の電子密度とエネルギー(2次元)

 投稿者:mike  投稿日:2017年 9月 6日(水)12時20分34秒
返信・引用
   定常状態の多電子系の電子密度とエネルギーを求める2次元のプログラム(024sampleLDA2D.bas)を公開します。
1次元定常状態の多電子系の電子密度とエネルギーを求める2次元のプログラム(022sampleLDA1D.bas)の2次元版です。

表示の説明:
オリーブ色は電子密度を表します。[D]を押すことで表示を変更できます。
[1] ~ [6] を押すと、系の中の電子の個数が変更できます。
表示のEnは電子軌道のエネルギー、Occは軌道の占有率です。

試験環境:
本プログラムは十進BASIC 6.6.3.3 / macOS 10.7.5 でテストしました。


プログラムが500行を超えたため、2つに分割して投稿します。


便宜のため、下記のURLでgoogle driveの一部を公開設定しました。
https://drive.google.com/drive/folders/0B6IBSbrPZRAVbFgxQU1LeE0yLUk
公開期限はありませんので、いつでもダウンロードできると思います。
本プログラムは024sampleLDA2Dv1.basです。
今までに投稿したプログラム([001]~[023])もあります。


---------------

!
! ========= RSDFT - local density approximation 2D ==========
!
! 024sampleLDA2D.bas
!   Copyright(C) 2017 Mitsuru Ikeuchi
!   Released under the MIT license ( https://opensource.org/licenses/MIT )
!
! ver 0.0.1   2017.09.06  created
!
OPTION ARITHMETIC NATIVE
DECLARE EXTERNAL SUB lda2d.setInitialCondition,lda2d.iterateLDA,lda2d.drawState,lda2d.setNumberOfElectron
DECLARE EXTERNAL FUNCTION INKEY$

LET stateMax = 10  !state 0,1,...,stateMax-1
LET vIndex = 0     !0:harmonic potential,  1:quantum well
LET nElectron = 4
LET iterMax = 5    !5 = iteration in iterateLDA
LET dispMode = 0   !0:Vext+rho, 1:Vext+rho+orbit, 2:rho+Vext+Veff+Vh+Vx+Vc
CALL setInitialCondition(stateMax,vIndex,nElectron)

DO
   CALL iterateLDA(stateMax, iterMax)
   CALL drawState(dispMode)
   LET S$=INKEY$
   IF S$="D" OR S$="d" THEN
      LET dispMode = mod(dispMode+1,3)
   ELSEIF S$="1" OR S$="2" OR S$="3" OR S$="4" OR S$="5" OR S$="6" THEN
      LET nElectron = VAL(s$)
      CALL setNumberOfElectron(nElectron)
   ELSEIF S$="7" THEN
      LET vIndex = 0
      CALL setInitialCondition(stateMax,vIndex,nElectron)
   ELSEIF S$="8" THEN
      LET vIndex = 1
      CALL setInitialCondition(stateMax,vIndex,nElectron)
   END IF
LOOP UNTIL S$=chr$(27)
END

EXTERNAL FUNCTION INKEY$
OPTION ARITHMETIC NATIVE
SET ECHO "OFF"
LET S$=""
CHARACTER INPUT NOWAIT: S$
LET INKEY$=S$
END FUNCTION

! ---------- RSDFT - local density approximation 2D ----------
!
! - real space density functional theory - local density approximation
! - solve Kohn-Sham equation - successive approximation
! - Vxc : LDA(local density approximation)
!         J. P. Perdew and A. Zunger; Phys. Rev., B23, 5048 (1981)
!
! procedure
! (1) given: trial |i>, occupation(i)
!
! (2) set electron density rho
!     rho <-- |i>, occupation[(i), mixing rho(iter-1)
!
! (3) set effective potential
!      Veff = Vext + Vh + Vx + Vc
!      Vh <-- rho (Poisson eq. ,SOR iteration)
!      Vx,Vc <-- rho (LDA:Perdew-Zunger)
!
! (4) solve Kohn-Sham equation (successive approximation)
!      |i> steepest descent method: |i(next)> = |i> - dump{H-E}|i>
!      E(i) <-- <i|H|i>
!      {|0>,..,|i>,..,|N>} orthogonallization : Gram-Schmidt
!
! (5) sort state
!      sort orbit by E(i)
!
! (6) set occupation
!      occupation(i) <-- E(i)
!
! (7) goto (2)
!
MODULE lda2d
MODULE OPTION ARITHMETIC NATIVE
OPTION BASE 0
PUBLIC SUB setInitialCondition, setNumberOfElectron, iterateLDA, drawState
SHARE NUMERIC NNx,NNy, dx,dy, lz, iterCount
SHARE NUMERIC numberOfElectron, numberOfElectronBounds, numberOfOrbit, errorDecisionOrbit
SHARE NUMERIC energyMem, iterationError, convergenceErrorMax, dampingFactor
SHARE NUMERIC mixing,broadening
SHARE NUMERIC cosAx,sinAx,cosAy,sinAy,cx0,cy0,cz0 !--- 3D graphics
SHARE NUMERIC sdEnergy(20)        ! electron state energy
SHARE NUMERIC sdState(20,200,200) ! electron states 0...20 0:ground state
SHARE NUMERIC occupation(20)      ! occupation of orbit
SHARE NUMERIC wrk(200,200)        ! state work space in steepestDescent
SHARE NUMERIC vv(200,200)         ! effective potential
SHARE NUMERIC vvext(200,200)      ! external potential
SHARE NUMERIC vvh(200,200)        ! Hartree potential
SHARE NUMERIC vvx(200,200)        ! exchange potential
SHARE NUMERIC vvc(200,200)        ! correlation potential
SHARE NUMERIC rho(200,200)        ! charge density
SHARE NUMERIC psicol(64,64)       ! MAT PLOT CELL matrix for psi(x,y)
SHARE NUMERIC vcol(64,64)         ! MAT PLOT CELL matrix  for potential V(x,y)
LET NNx = 64                  ! max number of sdState(,NNx,NNy)
LET NNy = 64                  !
LET dx = 0.25                 ! (au) x-division
LET dy = 0.25
LET lz = 12.0                 ! (au) v=dx*dy*lz
LET iterCount = 0             ! sd iteration count
LET numberOfElectron = 4      !
LET numberOfElectronBounds = 6! selection bounds OF numberOfElectron
LET numberOfOrbit = 10        !
LET energyMem = 0.0
LET iterationError = 1.0
LET convergenceErrorMax = 1.0e-5
LET dampingFactor = 0.01     ! steepest descent method
LET mixing = 0.5             ! charge mixing in setRho()
LET broadening = 0.01        ! (au) level broadening IN setOccupation

! ---------- set initial condition

EXTERNAL SUB setInitialCondition(stateMax,vIndex,nElectron)   !public
   DECLARE EXTERNAL SUB setInitialState,setExternalPotential,initDraw
   LET iterCount = 0
   LET numberOfElectron = nElectron
   CALL setInitialState(stateMax)
   CALL setExternalPotential(vIndex)
   CALL initDraw
   ! set window
   SET WINDOW 0,500, 0,500
END SUB

EXTERNAL SUB setNumberOfElectron(nElectron)   !public
   LET numberOfElectron = nElectron
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 NNx-2
            LET sdState(ist,i,j) = RND-0.5
         NEXT j
      NEXT i
      CALL normalizeState(ist)
   NEXT ist
END SUB

EXTERNAL SUB setExternalPotential(vIndex)
   LET x0 = 0.5*NNx*dx
   LET y0 = 0.5*NNy*dy
   FOR i=0 TO NNx-1
      FOR j=0 TO NNy-1
         LET x = i*dx
         LET y = j*dy
         IF vIndex=0 THEN
            LET vvext(i,j) = 0.5*((x-x0)*(x-x0)+(y-y0)*(y-y0))
         ELSEIF vIndex=1 THEN
            LET r = SQR((x-x0)*(x-x0)+(y-y0)*(y-y0))
            IF r>5 THEN LET vvext(i,j)=20 ELSE LET vvext(i,j)=0
         ELSE
            LET vvext(i,j)=0
         END IF
      NEXT j
   NEXT i
END SUB

EXTERNAL SUB setInitialOccupation(nOrbit, nElectron)
   LET occ = 1.0*nElectron/nOrbit
   FOR iState=0 TO nOrbit
      LET occupation(iState) = occ
   NEXT iState
END SUB

EXTERNAL SUB initDraw !--- set set color pallet and MAT PLOT matrix vcol(,)
   FOR i = 0 TO 50 !--- set color pallet
      SET COLOR MIX( 40+i) 0.02*i,0.02*i,0     ! red for |psi> +
      SET COLOR MIX(100+i) 0,0.02*i,0.02*i     ! blue for |psi> -
      SET COLOR MIX(160+i) 0,0.02*i,0     ! green for V(x)
   NEXT i
   FOR i=0 TO NNx-1 !--- set vcol(,)
      FOR j=0 TO NNy-1
         LET col = 0.02*vvext(i,j)
         IF col>1 THEN LET col = 1
         IF col<0 THEN LET col = 0
         LET vcol(i,j) = 160+INT(col*50)
      NEXT j
   NEXT i
END SUB

! ---------- iterate LDA

EXTERNAL SUB iterateLDA(stateMax, iterMax) !public
   DECLARE EXTERNAL FUNCTION steepestDescent
   DECLARE EXTERNAL SUB setElectronDensity,setEffectivePotential,solveKohnSham,sortState,setOccupation
   LET errorDecisionOrbit = int((numberOfElectron-1)/2)
   CALL setElectronDensity
   CALL setEffectivePotential
   CALL solveKohnSham(numberOfOrbit,iterMax)
   CALL sortState(numberOfOrbit)
   CALL setOccupation(numberOfOrbit,numberOfElectron)
   LET iterationError = sdEnergy(errorDecisionOrbit) - energyMem
   LET energyMem = sdEnergy(errorDecisionOrbit)
END SUB

!--- (2) set electron density rho <-- sdState(), occupation()
EXTERNAL SUB setElectronDensity
   FOR i=0 TO NNx-1
      FOR j=0 TO NNy-1
         LET rho(i,j) = rho(i,j)*(1.0-mixing)
         FOR ie=0 TO numberOfOrbit-1
            IF occupation(ie)>0.0 THEN
               LET rho(i,j) = rho(i,j) + mixing*occupation(ie)*(sdState(ie,i,j)*sdState(ie,i,j))/lz
            END IF
         NEXT ie
      NEXT j
   NEXT i
END SUB

!--- (3) set effective potential <-- electron density
EXTERNAL SUB setEffectivePotential
   DECLARE EXTERNAL SUB poisson,setVxc

   CALL poisson(20) !setVh
   CALL setVxc
   FOR i=0 TO NNx-1
      FOR j=0 TO NNy-1
         LET vv(i,j) = vvext(i,j)+vvh(i,j)+vvx(i,j)+vvc(i,j)
      NEXT j
   NEXT i
END SUB

EXTERNAL SUB poisson(iterMax)
   LET h2 = 4.0*PI*dx*dx
   LET omegav4 = 1.8/4.0
   FOR iter=0 TO iterMax-1
      FOR i=1 TO NNx-2
         FOR j=1 TO NNy-2
            LET vvh(i,j) = vvh(i,j)+omegav4*(vvh(i+1,j)+vvh(i-1,j)+vvh(i,j+1)+vvh(i,j-1)-4.0*vvh(i,j) +h2*rho(i,j))
         NEXT j
      NEXT i
   NEXT iter
END SUB

EXTERNAL SUB setVxc  !set exchage and correlation potential (Perdew and Zunger)
   LET c1 = -0.984745022
   FOR i=0 TO NNx-1
      FOR j=1 TO NNy-2
         LET rh = rho(i,j)
         LET rh3 = rh^0.33333333
         LET vvx(i,j) = c1*rh3
         LET rs = 0.6204/(rh3+1.0e-20)
         IF rs>=1.0 THEN
            LET sqrtrs = SQR(rs)
            LET ec = -0.1423/(1.0+1.0529*sqrtrs+0.3334*rs)
            LET vvc(i,j) = ec*(1.0+1.22838*sqrtrs+0.4445*rs)/(1.0+1.0529*sqrtrs+0.3334*rs)
         ELSE
            LET vvc(i,j) = -0.05837-0.0084*rs +(0.0311+0.00133*rs)*LOG(rs)
         END IF
      NEXT j
   NEXT i
END SUB

EXTERNAL FUNCTION eeCorrelation(rh)
   LET r = 0.6204/(rh^0.33333333+1.0e-20)
   IF r>=1.0 THEN
      LET ec = -0.1423/(1.0+1.0529*SQR(r)+0.3334*r)
   ELSE
      LET ec = -0.0480-0.0116*r+(0.0311+0.0020*r)*LOG(r)
   END IF
   LET eeCorrelation = ec
END FUNCTION

!--- (4) solve Kohn-Sham equation
EXTERNAL SUB solveKohnSham(stateMax, iterMax)
   DECLARE EXTERNAL FUNCTION steepestDescent
   DECLARE EXTERNAL SUB GramSchmidt,sortState
   FOR i=0 TO iterMax-1
      FOR ist=0 TO stateMax-1
         LET sdEnergy(ist) = steepestDescent(ist,dampingFactor)
      NEXT ist
      CALL GramSchmidt(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
         LET pij = sdState(ist,i,j)
         LET wrk(i,j) = (4*pij-sdState(ist,i+1,j)-sdState(ist,i-1,j)-sdState(ist,i,j-1)-sdState(ist,i,j+1))/h2+(vv(i,j)-ei)*pij
      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
         LET sdState(ist,i,j) = sdState(ist,i,j)-damp*wrk(i,j)
      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
         LET pij = sdState(ist,i,j)
         LET s = s+pij*((4.0*pij-sdState(ist,i+1,j)-sdState(ist,i-1,j)-sdState(ist,i,j-1)-sdState(ist,i,j+1))/h2+vv(i,j)*pij)
         LET sn = sn+pij*pij
      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
               LET sdState(istate,i,j) = sdState(istate,i,j) - s*sdState(ist,i,j)
            NEXT j
         NEXT i
      NEXT ist
      CALL normalizeState(istate)
   NEXT iState
END SUB

!--- (5) sort state
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
               LET w = sdState(ist,i,j)
               LET sdState(ist,i,j) = sdState(ist+1,i,j)
               LET sdState(ist+1,i,j) = w
            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

!--- (6) set occupation
EXTERNAL SUB setOccupation(stateMax, nElectron)
   DECLARE EXTERNAL FUNCTION trialOcc,FermiDirac
   LET eUpper = sdEnergy(stateMax-1)+1.0
   LET eLower = sdEnergy(0)-1.0
   FOR i=0 TO stateMax-1
      IF sdEnergy(i)>eUpper THEN LET eUpper = sdEnergy(i)
      IF sdEnergy(i)<eLower THEN LET eLower = sdEnergy(i)
   NEXT i

   DO WHILE (eUpper-eLower>1.0e-12)
      LET eFermi = (eUpper+eLower)/2.0
      LET ntrial = trialOcc(stateMax, eFermi)
      IF ntrial<nElectron THEN
         LET eLower = eFermi
      ELSE
         LET eUpper = eFermi
      END IF
   LOOP
   LET eFermi = (eUpper+eLower)/2.0
   FOR i=0 TO stateMax-1
      LET occupation(i) = 2.0*FermiDirac(sdEnergy(i), eFermi)
      IF (occupation(i)<0.0001) THEN LET occupation(i) = 0.0
      IF (2.0-occupation(i)<0.0001) THEN LET occupation(i) = 2.0
   NEXT i
END SUB

EXTERNAL FUNCTION trialOcc(stateMax, eFermi)
   DECLARE EXTERNAL FUNCTION FermiDirac
   LET s = 0.0
   FOR i=0 TO stateMax-1
      LET s = s + 2.0*FermiDirac(sdEnergy(i), eFermi)
   NEXT i
   LET trialOcc = s
END FUNCTION

EXTERNAL FUNCTION FermiDirac(ee, ef)
   LET et = broadening !level width
   LET a = (ee-ef)/et
   IF a>100 THEN LET ret = 0.0 ELSE LET ret = 1.0/(EXP(a)+1.0)
   LET FermiDirac = ret
END FUNCTION

! ---------- utility

EXTERNAL FUNCTION innerProduct(ist,jst)  !--- <ist|jst>
   LET s = 0
   FOR i=1 TO NNx-2
      FOR j=1 TO NNy-2
         LET s = s + sdState(ist,i,j)*sdState(jst,i,j)
      NEXT j
   NEXT i
   LET innerProduct = s*dx*dy
END FUNCTION

EXTERNAL SUB normalizeState(ist)
   LET s = 0
   FOR i=1 TO NNx-2
      FOR j=1 TO NNy-2
         LET s = s + sdState(ist,i,j)*sdState(ist,i,j)*dx*dy
      NEXT j
   NEXT i
   LET a = SQR(1/s)
   FOR i=1 TO NNx-2
      FOR j=1 TO NNy-2
         LET sdState(ist,i,j) = a*sdState(ist,i,j)
      NEXT j
   NEXT i
END SUB

EXTERNAL FUNCTION totalEnergy
   DECLARE EXTERNAL FUNCTION eeCorrelation
   LET sei = 0.0
   FOR i=0 TO numberOfOrbit-1
      LET sei = sei + occupation(i)*sdEnergy(i)
   NEXT i
   LET s = 0.0
   FOR i=1 TO NNx-2
      FOR j=1 TO NNy-2
         LET s = s + (-0.5*vvh(i,j)-0.25*vvx(i,j)+eeCorrelation(rho(i,j))-vvc(i,j))*rho(i,j)
      NEXT j
   NEXT i
   LET s = s*dx*dy
   LET totalEnergy = sei + s
END FUNCTION

! ---------- 3D graphics
 

Re: 関数実行するたびカウントアップ

 投稿者:Bsitumonn  投稿日:2017年 8月26日(土)20時34分42秒
返信・引用
  > No.4397[元記事へ]

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

丁寧な説明ありがとうございます。
 

Re: 関数実行するたびカウントアップ

 投稿者:SHIRAISHI Kazuo  投稿日:2017年 8月26日(土)09時31分2秒
返信・引用
  関数segmentを呼ぶ度にjを1つずつ増やす方法を教えてください。

jをモジュールのPUBLIC変数にすれば可能です。

主プログラムに
DECLARE EXTERNAL SUB m.segment
DECLARE EXTERNAL NUMERIC m.j
を書くと,以後,jと書くだけでモジュール内の変数jの参照が可能です。

そして,副プログラムsegmentをモジュールm内に書いてください。

MODULE m
DECLARE PUBLIC SUB segment
DECLARE PUBLIC NUMERIC j
LET j=0       ! 主プログラム実行前に実行される
EXTERNAL SUB segment(z1,z2,col)
・・・・・・
LET j=j+1
END SUB
END MODULE


なお,Full BASICの関数は引数はすべて値渡しですが,
副プログラムは実引数に変数を書くと参照渡しになるで,
副プログラムなら
EXTERNAL SUB segment(z1,z2,col,j)
・・・・・・
LET j=j+1
END SUB
のようにjを引数に追加しておく手も使えます。


2個目の質問
「cinputは、DECLARE EXTERNALの宣言がなくて、なぜいきなり呼び出せるでしょうか。」
は,オプションメニューの「互換性」の設定によります。JIS互換を選択すれば,DECLARE EXTERNAL SUB宣言の省略が許されなくなります。
 

関数実行するたびカウントアップ

 投稿者:Bsitumonn  投稿日:2017年 8月25日(金)21時10分46秒
返信・引用
  初めてこの掲示板を利用します。
失礼なこともあるかもしれませんが、よろしくお願いします。

関数segmentを呼ぶ度にjを1つずつ増やす方法を教えてください。
下記がおおまかなプログラムです。

OPTION ARITHMETIC complex
LET  j=1
DECLARE EXTERNAL SUB segment
LET w=5

SET WINDOW -w,w,-w,w
DRAW grid(0.5,0.5)
PRINT "複素平面上の三角形を関数 w=z^3 で写す"
PRINT "3点の座標を入力してください。"

CALL cinput(z1)
CALL cinput(z2)
CALL cinput(z3)
CALL segment(z1,z2,2)
CALL segment(z2,z3,3)
CALL segment(z3,z1,4)

END

EXTERNAL SUB segment(z1,z2,col)
OPTION ARITHMETIC COMPLEX
DEF f(z)=z^3
LET  i=SQR(-1)

SET LINE COLOR col
REM 線分を描く
SET LINE width 3
PLOT LINES

REM 線分上の点の写像による像を描く
SET LINE width 2

FOR t=0 TO 1 STEP 0.2
   LET  z=(1-t)*z1+t*z2
   PRINT z,j
   PLOT LINES : re(z),im(z);

NEXT t

LET j=j+1

END SUB

REM 複素数の入力
EXTERNAL SUB cinput(z)
OPTION ARITHMETIC COMPLEX
INPUT PROMPT "実部虚部を入力: ": x, y
LET  z=complex(x,y)
END SUB

質問ばかりで悪いのですが、cinputは、DECLARE EXTERNALの宣言がなくて、なぜいきなり呼び出せるでしょうか。教えていただけると幸いです。
お返事おねがいします。
 

[023]イオン結晶(2次元)

 投稿者:mike  投稿日:2017年 8月 4日(金)09時33分26秒
返信・引用
  イオン結晶の分子動力学法プログラム023NaClIonMD2D.basを公開します。
 NaClのような陽イオンと陰イオンからなる系において、電子殻の相互作用による引力や斥力に加え、静電的な
引力や斥力が働きます。静電気力は1/r^2に比例するため、比較的遠くまでその力が及びます。このため、分子動力学の
高速化(O(N))に必要な近くの粒子だけで力の計算をすることが難しくなります。
 しかしながら、多くの場合、陽イオンと陰イオンは近くに存在し、中性に近くなるので、互いに静電気力を打ち消す方向に
働くため、静電気力は急激に減少します。これをDebye遮蔽(1/r^2)*exp(-ar)といいます。
 本プログラムではDebye遮蔽を前提として力の計算をしています。(粒子登録法によるO(N)の高速化を行っています)

表示の説明:
赤色は陽イオン、青色は陰イオンをあらわします。
[K]を押し、温度を上げていくと結晶は不安定になりやがて融解します。(融点は2次元のため現実とは異なります)
[J]を押し、温度を下げていくと再び結晶に戻ります。
[1]-[7]はイオン結晶の種類を選択します。
表示は[D]のキーを押すことで変更できます。

試験環境:
本プログラムは十進BASIC 6.6.3.3 / macOS 10.7.5 でテストしました。

----------


!
! ========= molecular dynamics 2D - ion ==========
!
! 023NaClIonMD2D.bas
!   Copyright(C) 2017 Mitsuru Ikeuchi
!   Released under the MIT license ( https://opensource.org/licenses/MIT )
!
! ver 0.0.1   2017.08.03  created
!
OPTION ARITHMETIC NATIVE
DECLARE EXTERNAL SUB imd2d.setInitialCondition, imd2d.moveParticles, imd2d.drawParticles
DECLARE EXTERNAL FUNCTION INKEY$
LET tempMode = 0    !tempMode: 0:adiabatic  1:constant temperature
LET contTemp = 300  !contTemp: controled constant temperature(K)
LET ddTemp = 0      !contTemp = contTemp + ddTemp
LET drawMode = 1    !drawMode: 0:bond  1:circle+bond 2:velocitySpace
LET pauseFlag = 0       !if pauseFlag=0, CALL moveParticles(tempMode,contTemp)

!setInitialCondition(material,boxSizeInNM,contTemp)
CALL setInitialCondition(1,6.0,contTemp) !material = 1:NaCl 2:MgO 3:CaO 4:BaO 5:NaF 6:KF 7:KCl

DO
   LET contTemp = contTemp + ddTemp
   IF contTemp>3000 THEN LET contTemp = 3000
   IF contTemp<10 THEN LET contTemP = 10
   IF pauseFlag=0 THEN CALL moveParticles(tempMode,contTemp) ELSE WAIT DELAY 0.05
   CALL drawParticles(tempMode,contTemp,drawMode)
   LET S$=INKEY$
   IF S$="1" OR S$="2" OR S$="3" OR S$="4" OR S$="5" OR S$="6" OR S$="7" THEN
      LET material = VAL(S$)
      LET tempMode = 0
      LET contTemp = 300
      CALL setInitialCondition(material,6.0,contTemp)
   ELSEIF S$="T" OR S$="t" THEN
      LET tempMode = MOD(tempMode+1,2)
      IF tempMode=0 THEN LET ddTemp = 0
   ELSEIF S$="K" OR S$="k" THEN
      LET tempMode = 1
      IF ddtemp=0 THEN LET ddTemp = 1 ELSE LET ddTemp = 0
   ELSEIF S$="J" OR S$="j" THEN
      LET tempMode = 1
      IF ddTemp=0 THEN LET ddTemp = -1 ELSE LET ddTemp = 0
   ELSEIF S$="D" OR S$="d" THEN
      LET drawMode = MOD(drawMode+1,3)
   ELSEIF S$=" " THEN
      LET pauseFlag = MOD(pauseFlag+1,2)
   END IF
LOOP UNTIL S$=CHR$(27)
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

! ---------- molecular dynamics 2D - ion ----------
!
!  method: velocity Verlet ( F=m*d^2r/dt^2 -> r(t+dt)=r(t)+v*dt,v=v+(F/m)*dt )
!    (1) vi = vi + (Fi/mi)*(0.5dt)
!    (2) ri = ri + vi*dt
!    (3) calculation Fi <- {r1,r2,...,rn} Fi=sum(Fij,j=1 to n),Fij=F(ri-rj)
!    (4) vi = vi + (Fi/mi)*(0.5dt)
!    (5) goto (1)
!
!  force: ion f(r) = fc + fr + fa
!        fc = eForceConst*zi*zj*(EXP(-r/6.5e-10)/r)*(1.0/r+1.0/6.5e-10) !Debye-screened Coulomb force
!        fr = 6.9742e-11*EXP((a-r)/b) !repulsive force
!        fa = -6.9742e-21*(c/r^6) !attractive force
!
MODULE imd2d
MODULE OPTION ARITHMETIC NATIVE
PUBLIC SUB setInitialCondition !(molKind,boxSizeInNM,xtalSizeInNM,contTemp)
PUBLIC SUB moveParticles !(tempMode,contTemp)
PUBLIC SUB drawParticles !(drawMode:0:realSpace 1:velocitySpace)
SHARE NUMERIC ionKind1,ionKind2, sysTime, dt, nMolec, xMax, yMax, rCutoff, hh
SHARE NUMERIC xx(2000),yy(2000)             ! (xx(i),yy(i))  : position of i-th particle
SHARE NUMERIC vx(2000),vy(2000)             ! (vx(i),vy(i))  : velocity of i-th particle
SHARE NUMERIC ffx(2000),ffy(2000)           ! (ffx(i),ffy(i)): total force of i-th particle
SHARE NUMERIC kind(2000),mas(2000)          ! kind(i),mas(i) : molec kind, mass of i-th particle
SHARE NUMERIC reg(2000,0 TO 100)            ! register near molec reg(i,0):number of near i-th molec
SHARE NUMERIC ion_m(0 TO 17),ion_z(0 TO 17) ! ion mass,ion charge
SHARE NUMERIC ion_a(0 TO 17),ion_b(0 TO 17) ! ion force param a,ion force param b
SHARE NUMERIC ion_c(0 TO 17),ion_r(0 TO 17) ! ion force param c, ion radius
SHARE NUMERIC ion_col(0 TO 17)              ! ion draw color
SHARE STRING ion_str$(0 TO 17)              ! ion string
SHARE NUMERIC forceTable(0 TO 17, 0 TO 17,0 TO 1001) ! force table

LET ionKind1 = 3           ! ion kind1 - 3:Na+
LET ionKind2 = 7           ! ion kind2 - 7:Cl-
LET sysTime = 0.0          ! system time (s) in the module
LET dt = 2.0*1.0e-15       ! time step (s)
LET nMolec = 100           ! number of particles
LET xMax = 6.0E-9          ! x-Box size (m)
LET yMax = 6.0E-9          ! y-Box size (m)
LET rCutoff = 1.0e-9       ! force cutoff radius (m)
LET hh = 1e-12             ! force table step (m)
LET ion_z(0)=0.0           ! if ion_z(0)=0.0 then set ion data

! ---------- set initial condition

EXTERNAL SUB setInitialCondition(material,boxSizeInNM,contTemp)
   DECLARE EXTERNAL SUB setIonData,ajustVelocity,setForceTable
   DECLARE EXTERNAL FUNCTION setNaClTypeBlock
   RANDOMIZE
   LET sysTime = 0.0
   CALL setIonData
   CALL setForceTable
   LET xMax = boxSizeInNM*1.0e-9
   LET yMax = boxSizeInNM*1.0e-9
   IF material=1 THEN !NaCl
      LET ionKind1 = 3 !Na+
      LET ionKind2 = 7 !Cl-
      LET lattice = 5.6407e-10
   ELSEIF material=2 THEN !MgO
      LET ionKind1 = 4 !Mg++
      LET ionKind2 = 1 !O--
      LET lattice = 4.212e-10*0.94 !0.94: correction factor
   ELSEIF material=3 THEN !CaO
      LET ionKind1 = 9 !Ca++
      LET ionKind2 = 1 !O--
      LET lattice = 4.80e-10*0.94 !0.94: correction factor
   ELSEIF material=4 THEN !BaO
      LET ionKind1 = 12 !Ba++
      LET ionKind2 = 1 !O--
      LET lattice = 5.536e-10
   ELSEIF material=5 THEN !NaF
      LET ionKind1 = 3 !Na+
      LET ionKind2 = 2 !F-
      LET lattice = 4.62e-10
   ELSEIF material=6 THEN !KF
      LET ionKind1 = 8 !K+
      LET ionKind2 = 2 !F-
      LET lattice = 5.34e-10
   ELSEIF material=7 THEN !KCl
      LET ionKind1 = 8 !K+
      LET ionKind2 = 7 !Cl-
      LET lattice = 6.29e-10
   END IF
   LET s = 0.5*(xMax - lattice*6)
   !            setNaClTypeBlock(ii, knd1, knd2, nx, ny, lattice, xPos, yPos)
   LET nMolec = setNaClTypeBlock(1, ionKind1, ionKind2, 6, 6, lattice, s, s)
   CALL ajustVelocity(contTemp)
   ! set window
   SET WINDOW 0,500,0,500
END SUB

EXTERNAL SUB setIonData
! ion potential data
!      0 mass,1 charge,  2  a ,     3  b , 4 c , 5  r-ion, 6 color 7 str$
   DATA 10.81,  3.0, 0.720e-10, 0.080e-10,  0.0, 0.23e-10, 13 , "B+++"   !0
   DATA 16.00, -2.0, 1.626e-10, 0.085e-10, 20.0, 1.40e-10,  2 , "O--"    !1
   DATA 19.00, -1.0, 1.565e-10, 0.085e-10, 20.0, 1.33e-10,  2 , "F-"     !2
   DATA 22.99,  1.0, 1.260e-10, 0.080e-10, 20.0, 1.02e-10,  4 , "Na+"    !3
   DATA 24.31,  2.0, 1.161e-10, 0.080e-10, 10.0, 0.72e-10,  4 , "Mg++"   !4
   DATA 26.98,  3.0, 1.064e-10, 0.080e-10,  2.0, 0.53e-10, 13 , "Al+++"  !5
   DATA 28.09,  4.0, 1.012e-10, 0.080e-10,  0.0, 0.40e-10,  7 , "Si++++" !6
   DATA 35.45, -1.0, 1.950e-10, 0.090e-10, 30.0, 1.81e-10,  2 , "Cl-"    !7
   DATA 39.10,  1.0, 1.595e-10, 0.080e-10, 15.0, 1.38e-10,  4 , "K+"     !8
   DATA 40.08,  2.0, 1.414e-10, 0.080e-10, 10.0, 1.00e-10,  4 , "Ca++"   !9
   DATA 47.88,  4.0, 1.235e-10, 0.080e-10,  0.0, 0.61e-10,  7 , "Ti++++" !10
   DATA 87.62,  2.0, 1.632e-10, 0.080e-10, 15.0, 1.16e-10,  4 , "Sr++"   !11
   DATA 137.3,  2.0, 1.820e-10, 0.080e-10, 20.0, 1.36e-10,  4 , "Ba++"   !12
   DATA 4.003,  0.0, 1.200e-10, 0.110e-10, 4.76, 1.28e-10,  3 , "He"     !13
   DATA 20.18,  0.0, 1.415e-10, 0.112e-10,11.03, 1.37e-10,  3 , "Ne"     !14
   DATA 39.95,  0.0, 1.878e-10, 0.117e-10,38.53, 1.70e-10,  3 , "Ar"     !15
   DATA 83.80,  0.0, 2.041e-10, 0.130e-10,55.33, 1.83e-10,  3 , "Kr"     !16
   DATA 131.3,  0.0, 2.258e-10, 0.145e-10,85.55, 1.99e-10,  3 , "Xe"     !17

   IF ion_z(0)=0.0 THEN
      FOR i=0 TO 17
         READ ion_m(i),ion_z(i),ion_a(i),ion_b(i),ion_c(i),ion_r(i),ion_col(i),ion_str$(i)
         LET ion_m(i) = ion_m(i)*1.661e-27
         LET ion_z(i) = ion_z(i)*1.602e-19
      NEXT i
   END IF
END SUB

EXTERNAL FUNCTION setNaClTypeBlock(ii, knd1, knd2, nx, ny, lattice, xPos, yPos)
   DECLARE EXTERNAL SUB setParticle
   LET ipp = ii
   LET a = lattice/2.0
   FOR i=0 TO 2*nx-1
      FOR j=0 TO 2*ny-1
         LET x = xPos + a*i
         LET y = yPos + a*j
         IF MOD(i+j,2)=0 THEN
            LET knd = knd1
         ELSE
            LET knd = knd2
         END IF
         CALL setParticle(ipp, knd, x, y)
         LET ipp = ipp + 1
      NEXT j
   NEXT i
   LET setNaClTypeBlock = ipp - 1
END FUNCTION

EXTERNAL SUB setParticle(i, knd, x, y)
   LET xx(i) = x
   LET yy(i) = y
   LET vx(i) = 200.0*(RND+RND+RND+RND+RND+RND-3)
   LET vy(i) = 200.0*(RND+RND+RND+RND+RND+RND-3)
   LET ffx(i) = 0.0
   LET ffy(i) = 0.0
   LET mas(i) = ion_m(knd)
   LET kind(i) = knd
END SUB

EXTERNAL SUB setForceTable
   DECLARE EXTERNAL  FUNCTION cutoff
   LET eForceConst = 1.0/(4.0*PI*8.8542e-12) !epsilon0=8.8542e-12
   FOR ki=0 TO 17
      FOR kj=0 TO 17
         LET a = ion_a(ki)+ion_a(kj)
         LET b = ion_b(ki)+ion_b(kj)
         LET c = ion_c(ki)*ion_c(kj)*1.0e-60
         LET zi = ion_z(ki)
         LET zj = ion_z(kj)
         FOR ir=10 TO 1001
            LET r = ir*hh
            LET fc = eForceConst*zi*zj*(EXP(-r/6.5e-10)/r)*(1.0/r+1.0/6.5e-10) !Debye-screened Coulomb force
            LET fr = 6.9742e-11*EXP((a-r)/b) !repulsive force
            LET fa = -6.9742e-21*(c/(r*r*r*r*r*r)) !attractive force
            LET forceTable(ki,kj,ir) = cutoff(r)*(fc + fr + fa)
         NEXT ir
         FOR ir=0 TO 9
            LET forceTable(ki,kj,ir) = forceTable(ki,kj,10)
         NEXT ir
      NEXT kj
   NEXT ki
END SUB

EXTERNAL FUNCTION cutoff(r)
   IF r>0 AND r<0.8*rCutoff THEN
      LET ret = 1
   ELSEIF r>=0.8*rCutoff AND r<rCutoff THEN
      LET ret = 0.5+0.5*COS(PI*(r-0.8*rCutoff)/(0.2*rCutoff))
   else
      LET ret = 0
   END IF
   LET cutoff = ret
END FUNCTION

! ---------- move particles

EXTERNAL SUB moveParticles(tempMode,contTemp) !tempMode 0:adiabatic  1:constant-temp
   DECLARE EXTERNAL SUB moveParticlesDT,ajustVelocity,registerNearMolec
   IF (tempMode=1) THEN CALL ajustVelocity(contTemp)
   CALL registerNearMolec
   FOR i=1 TO 20
      CALL moveParticlesDT
   NEXT i
END SUB

EXTERNAL SUB moveParticlesDT ! velocity Verlet method
   DECLARE EXTERNAL SUB calcForce
   FOR i=1 TO nMolec
      LET a = 0.5*dt/mas(i)
      LET vx(i) = vx(i)+a*ffx(i)
      LET vy(i) = vy(i)+a*ffy(i)
      LET xx(i) = xx(i)+vx(i)*dt
      LET yy(i) = yy(i)+vy(i)*dt
   NEXT i
   CALL calcForce
   FOR i=1 TO nMolec
      LET a = 0.5*dt/mas(i)
      LET vx(i) = vx(i)+a*ffx(i)
      LET vy(i) = vy(i)+a*ffy(i)
   NEXT i
   LET sysTime=sysTime+dt
END SUB

EXTERNAL SUB calcForce
   DECLARE EXTERNAL FUNCTION force,boundaryForce
   LET s = 0.5*3.418e-10
   FOR i=1 TO nMolec
      LET ffx(i) = 0.0
      LET ffy(i) = 0.0
   NEXT i
   FOR i=1 TO nMolec-1
      FOR k=1 TO reg(i,0)-1
         LET j = reg(i,k)
         LET xij = xx(i)-xx(j)
         LET yij = yy(i)-yy(j)
         LET rij = SQR(xij*xij+yij*yij)
         IF rij<rCutoff THEN
            LET f = force(rij,kind(i),kind(j))
            LET fxij = f*xij/rij
            LET fyij = f*yij/rij
            LET ffx(i) = ffx(i)+fxij
            LET ffy(i) = ffy(i)+fyij
            LET ffx(j) = ffx(j)-fxij
            LET ffy(j) = ffy(j)-fyij
         END IF
      NEXT k
   NEXT i
   FOR i=1 TO nMolec  ! boundary force
      LET ffx(i) = ffx(i)+boundaryForce(xx(i)+s)+boundaryForce(xx(i)-xMax-s)
      LET ffy(i) = ffy(i)+boundaryForce(yy(i)+s)+boundaryForce(yy(i)-yMax-s)
   NEXT i
END SUB

EXTERNAL FUNCTION force(r,ki,kj) !force(r) <-- forceTable - linear interporation
   LET ir = INT(r/hh)
   LET a = r - ir*hh
   LET force = ((hh-a)*forceTable(ki,kj,ir) + a*forceTable(ki,kj,ir+1))/hh
END FUNCTION

EXTERNAL FUNCTION boundaryForce(r)
   LET r6 = (3.418e-10/r)^6
   LET boundaryForce = (24.0*(0.5*1.711e-21)*r6*(2.0*r6-1.0)/r)
END FUNCTION

EXTERNAL SUB registerNearMolec
   LET rCut = rCutoff+20*2000*dt
   LET rcut2 = rCut*rCut
   FOR i=1 TO nMolec-1
      LET k = 1
      FOR j=i+1 TO nMolec
         LET r2 = (xx(i)-xx(j))*(xx(i)-xx(j))+(yy(i)-yy(j))*(yy(i)-yy(j))
         IF (r2<rcut2) THEN
            LET reg(i,k) = j
            LET k = k + 1
         END IF
      NEXT j
      LET reg(i,0) = k
   NEXT i
END SUB

EXTERNAL FUNCTION maxNearMolec
   LET mx = 0
   FOR i=1 TO nMolec-1
      IF mx<reg(i,0) THEN LET mx = reg(i,0)
   NEXT i
   LET maxNearMolec = mx-1
END FUNCTION

! ---------- utility

EXTERNAL FUNCTION systemTemprature
   LET kB = 1.38e-23 ! Boltzman's constant (J/K)
   LET ek= 0.0       !kinetic energy (J)
   FOR i=1 TO nMolec
      LET ek = ek + 0.5*mas(i)*(vx(i)^2+vy(i)^2)
   NEXT i
   LET systemTemprature = ek/(nMolec*kB) !2D: E/N=kT, 3D: E/N=(3/2)kT
END FUNCTION

EXTERNAL SUB ajustVelocity(temp)
   DECLARE EXTERNAL FUNCTION systemTemprature
   LET r = sqr(temp/systemTemprature)
   FOR i=1 TO nMolec
      LET vx(i) = r*vx(i)
      LET vy(i) = r*vy(i)
   NEXT i
END SUB

! ---------- draw particles

EXTERNAL SUB drawParticles(tempMode,contTemp,drawMode)
   DECLARE EXTERNAL FUNCTION systemTemprature,maxNearMolec
   DECLARE EXTERNAL sub realSpace,velocitySpace,plotBond
   SET DRAW MODE HIDDEN
   CLEAR
   IF drawMode=0 OR drawMode=1 THEN !--- 0:circle  1:circle+bond
      call plotBond(drawMode)
   ELSEIF drawMode=2 THEN
      call velocitySpace
   END IF

   !--- draw caption
   SET TEXT HEIGHT 10
   SET TEXT COLOR 1 ! black
   LET tmp$ = "adiabatic   constantTemp  "
   PLOT TEXT, AT  50, 70 ,USING "time =#####.## (ps)   temp =####.## (K)":sysTime*1E12,systemTemprature
   PLOT TEXT, AT  50, 55 :ion_str$(ionKind1)&","&ion_str$(ionKind2)
   PLOT TEXT, AT 150, 55 ,USING "N =####":nMolec
   PLOT TEXT, AT  50, 40 ,USING "tempMode = ## ":tempMode
   PLOT TEXT, AT 200, 40 :tmp$(tempMode*12+1:tempMode*12+12)
   PLOT TEXT, AT  50, 25 ,USING "controled Temperature =####.# (K)":contTemp
   PLOT TEXT, AT  50, 10 :"2-dimensional ion - molecular dynamics"
   PLOT TEXT, AT  50,480 :"[esc] exit  [space]pause/go  [D]change draw mode"
   PLOT TEXT, AT  50,460 :"[1]NaCl [2]MgO [3]CaO [4]BaO [5]NaF [6]KF [7]KCl"
   PLOT TEXT, AT  50,440 :"[T] toggle 0:adiabatic/1:constant temperature"
   PLOT TEXT, AT  50,420 :"temp control -> [J]coolDown/stop  [k]heatUp/stop"
   SET DRAW MODE EXPLICIT
END SUB

EXTERNAL sub plotBond(drawMode)
   LET boxSize = 300
   LET mag = boxSize/xMax
   LET xp = 100
   LET yp = 100
   SET LINE COLOR 1 ! black : !--- box
   PLOT LINES: xp,yp; boxSize+xp,yp; boxSize+xp,boxSize+yp; xp,boxSize+yp; xp,yp
   SET TEXT HEIGHT 6
   SET TEXT COLOR 1 ! black
   PLOT TEXT, AT xp,boxSize+2+yp ,USING  "box size =##.# x ##.# (nm)":xMax*1e9,yMax*1e9
   FOR i=1 TO nMolec
      SET LINE COLOR ion_col(kind(i))
      DRAW circle WITH SCALE(ion_r(kind(i))*mag)*SHIFT(xx(i)*mag+xp,yy(i)*mag+yp)
   NEXT i
   IF drawMode=1 THEN
      FOR i=1 TO nMolec-1
         FOR k=1 TO reg(i,0)-1
            LET j = reg(i,k)
            LET r = SQR((xx(i)-xx(j))*(xx(i)-xx(j))+(yy(i)-yy(j))*(yy(i)-yy(j)))
            LET r0 = (ion_r(kind(i))+ion_r(kind(j)))
            IF r<1.1*r0 AND kind(i)<>kind(j) THEN
               SET LINE COLOR 8 !gray
               PLOT LINES: xx(i)*mag+xp,yy(i)*mag+yp;xx(j)*mag+xp,yy(j)*mag+yp
            END IF
         NEXT k
      NEXT i
   END IF
END sub

EXTERNAL sub velocitySpace
   LET boxSize = 300
   LET xp = 100
   LET yp = 100
   SET LINE COLOR 1 !black : axis
   PLOT LINES: xp,boxSize/2+yp; boxSize+xp,boxSize/2+yp !vx-axis
   PLOT LINES: boxSize/2+xp,yp; boxSize/2+xp,boxSize+yp !vy-axis
   SET TEXT HEIGHT 6
   SET TEXT COLOR 1 ! black
   PLOT TEXT, AT boxSize+xp,boxSize/2+yp: "vx"
   PLOT TEXT, AT boxSize+xp,boxSize/2-12+yp: "2000m/s"
   PLOT TEXT, AT boxSize/2-12+xp,boxSize+yp: "vy  2000m/s"
   PLOT TEXT, AT boxSize/2-8+xp,boxSize/2-10+yp: "0"
   PLOT TEXT, AT xp,boxSize+8+yp: "velocity space (vx,vy)"
   LET mag = boxSize/4000
   FOR i=1 TO nMolec
      IF kind(i)=ionKind1 THEN SET LINE COLOR 4 ELSE SET LINE COLOR 2  !4:red, 2:blue
      DRAW circle WITH SCALE(5)*SHIFT(vx(i)*mag+boxSize/2+xp,vy(i)*mag+boxSize/2+yp)
   NEXT i
END sub

END MODULE
 

[022]多電子系の電子密度とエネルギー

 投稿者:mike  投稿日:2017年 7月28日(金)09時45分33秒
返信・引用
  定常状態の多電子系の電子密度とエネルギーを求めるプログラム(022sampleLDA1D.bas)を公開します。
 1電子の定常状態の電子のエネルギーと電子状態を求める問題は、[008]最急降下法で解くことができました。
多電子系の場合、厳密に電子状態を求めるのは、電子数が多くなるとほとんど不可能になります。
ここでは、1つの電子に注目し他の電子は平均場として扱う近似により多電子系の問題を解く、密度汎関数法により
電子密度とエネルギーを求めます。
 密度汎関数法(Density Functional Theory)では、電子密度はそれぞれの電子軌道から決まる電子の存在確率の
総計から求められ、1電子はSchroedingerの方程式に似たKohn-Sham方程式(ポテンシャルの中に多電子の寄与がある)
に従うと仮定します。実効ポテンシャルVeffへの多電子の寄与は電子密度の関数になります。
        Veff = Vext + VH + Vx + Vc
Vextは外部ポテンシャル、VHは静電ポテンシャル、Vxは交換ポテンシャル、Vcは相関ポテンシャルです。
 本プログラムでは電子密度を仮定し、Kohn-Sham方程式を最急降下法で解き、各電子軌道の総計から電子密度を求め、
これを繰り返すことで電子密度と電子軌道、エネルギーを求めます。

表示の説明:
黒の線は外部ポテンシャル、灰色は電子密度を表します。[D]を押すことで表示を変更できます。

試験環境:
本プログラムは十進BASIC 6.6.3.3 / macOS 10.7.5 でテストしました。

----------------

!
! ========= RSDFT - local density approximation 1D ==========
!
! 022sampleLDA1D.bas
!   Copyright(C) 2017 Mitsuru Ikeuchi
!   Released under the MIT license ( https://opensource.org/licenses/MIT )
!
! ver 0.0.1   2017.07.28  created
!
OPTION ARITHMETIC NATIVE
DECLARE EXTERNAL SUB lda1d.setInitialCondition,lda1d.setNumberOfElectron,lda1d.iterateLDA,lda1d.drawState
DECLARE EXTERNAL FUNCTION INKEY$

LET stateMax = 10  !state 0,1,...,stateMax-1
LET vIndex = 0     !0:harmonic potential,  1:quantum well
LET nElectron = 4
LET iterMax = 10   ! 10 = iteration in iterateLDA
LET dispMode = 1   !0:Vext+rho, 1:Vext+rho+orbit, 2:rho+Vext+Veff+Vh+Vx+Vc
CALL setInitialCondition(stateMax,vIndex,nElectron)

DO
   CALL iterateLDA(stateMax, iterMax)
   CALL drawState(dispMode)
   LET S$=INKEY$
   IF S$="D" OR S$="d" THEN
      LET dispMode = mod(dispMode+1,3)
   ELSEIF S$="1" OR S$="2" OR S$="3" OR S$="4" OR S$="5" OR S$="6" THEN
      LET nElectron = VAL(s$)
      CALL setNumberOfElectron(nElectron)
   ELSEIF S$="7" THEN
      LET vIndex = 0
      CALL setInitialCondition(stateMax,vIndex,nElectron)
   ELSEIF S$="8" THEN
      LET vIndex = 1
      CALL setInitialCondition(stateMax,vIndex,nElectron)
   END IF
LOOP UNTIL S$=chr$(27)
END

EXTERNAL FUNCTION INKEY$ !from decimal BASIC library
OPTION ARITHMETIC NATIVE
SET ECHO "OFF"
LET S$=""
CHARACTER INPUT NOWAIT: S$
LET INKEY$=S$
END FUNCTION

! ---------- RSDFT - local density approximation 1D ----------
!
! - real space density functional theory - local density approximation
! - solve Kohn-Sham equation - successive approximation
! - Vxc : LDA(local density approximation)
!         J. P. Perdew and A. Zunger; Phys. Rev., B23, 5048 (1981)
!
! procedure
! (1) given: trial |i>, occupation(i)
!
! (2) set electron density rho
!     rho <-- |i>, occupation[(i), mixing rho(iter-1)
!
! (3) set effective potential
!      Veff = Vext + Vh + Vx + Vc
!      Vh <-- rho (Poisson eq. ,SOR iteration)
!      Vx,Vc <-- rho (LDA:Perdew-Zunger)
!
! (4) solve Kohn-Sham equation (successive approximation)
!      |i> steepest descent method: |i(next)> = |i> - dump{H-E}|i>
!      E(i) <-- <i|H|i>
!      {|0>,..,|i>,..,|N>} orthogonallization : Gram-Schmidt
!
! (5) sort state
!      sort orbit by E(i)
!
! (6) set occupation
!      occupation(i) <-- E(i)
!
! (7) goto (2)
!
MODULE lda1d
MODULE OPTION ARITHMETIC NATIVE
OPTION BASE 0
PUBLIC SUB setInitialCondition, setNumberOfElectron, iterateLDA, drawState
SHARE NUMERIC NNx, dx, lylz, iterCount
SHARE NUMERIC numberOfElectron, numberOfElectronBounds, numberOfOrbit, errorDecisionOrbit
SHARE NUMERIC energyMem, iterationError, convergenceErrorMax, dampingFactor
SHARE NUMERIC mixing,broadening
SHARE NUMERIC sdEnergy(20)    ! electron state energy
SHARE NUMERIC sdState(20,400) ! electron states 0...20 0:ground state
SHARE NUMERIC occupation(20)  ! occupation of orbit
SHARE NUMERIC wrk(400)        ! state work space in steepestDescent
SHARE NUMERIC vv(400)         ! effective potential
SHARE NUMERIC vvext(400)      ! external potential
SHARE NUMERIC vvh(400)        ! Hartree potential
SHARE NUMERIC vvx(400)        ! exchange potential
SHARE NUMERIC vvc(400)        ! correlation potential
SHARE NUMERIC rho(400)        ! charge density
LET NNx = 64                  ! max number of sdState(,NNx,NNy)
LET dx = 0.25                 ! (au) x-division
LET lylz = 16.0*16.0          ! (au) v=dx*ly*lz
LET iterCount = 0             ! sd iteration count
LET numberOfElectron = 4      !
LET numberOfElectronBounds = 6! selection bounds OF numberOfElectron
LET numberOfOrbit = 10        !
LET energyMem = 0.0
LET iterationError = 1.0
LET convergenceErrorMax = 1.0e-5
LET dampingFactor = 0.03     ! steepest descent method
LET mixing = 0.5             ! charge mixing in setRho()
LET broadening = 0.001       ! (au) level broadening IN setOccupation

! ---------- set initial condition

EXTERNAL SUB setInitialCondition(stateMax,vIndex,nElectron)   !public
   DECLARE EXTERNAL SUB setInitialState,setExternalPotential
   LET iterCount = 0
   LET numberOfElectron = nElectron
   CALL setInitialState(stateMax)
   CALL setExternalPotential(vIndex)
   ! set window
   SET WINDOW 0,500, 0,500
END SUB

EXTERNAL SUB setNumberOfElectron(nElectron)   !public
   LET iterCount = 0
   LET numberOfElectron = nElectron
END SUB

EXTERNAL SUB setInitialState(stateMax)
   DECLARE EXTERNAL SUB normalizeState
   RANDOMIZE
   FOR ist=0 TO stateMax-1
      FOR i=1 TO NNx-2
         LET sdState(ist,i) = RND-0.5
      NEXT i
      LET sdState(ist,0) = 0
      LET sdState(ist,NNx-1) = 0
      CALL normalizeState(ist)
   NEXT ist
END SUB

EXTERNAL SUB setExternalPotential(vIndex)
   LET x0 = 0.5*NNx*dx
   IF vIndex=0 THEN !--- hermonic
      FOR i=0 TO NNx-1
         LET x = i*dx
         LET vvext(i) = MIN(0.5*(x-x0)^2,24.5)
      NEXT i
   ELSEIF vIndex=1 THEN !--- well
      FOR i=0 TO NNx-1
         LET x = i*dx
         IF ABS(x-x0)<4 THEN LET vvext(i) = 0 ELSE LET vvext(i) = 18
      NEXT i
   END IF
END SUB

EXTERNAL SUB setInitialOccupation(nOrbit, nElectron)
   LET occ = 1.0*nElectron/nOrbit
   FOR iState=0 TO nOrbit
      LET occupation(iState) = occ
   NEXT iState
END SUB

! ---------- iterate LDA

EXTERNAL SUB iterateLDA(stateMax, iterMax) !public
   DECLARE EXTERNAL FUNCTION steepestDescent
   DECLARE EXTERNAL SUB setElectronDensity,setEffectivePotential,solveKohnSham,sortState,setOccupation
   LET errorDecisionOrbit = (numberOfElectron-1)/2
   CALL setElectronDensity
   CALL setEffectivePotential
   CALL solveKohnSham(numberOfOrbit,iterMax)
   CALL sortState(numberOfOrbit)
   CALL setOccupation(numberOfOrbit,numberOfElectron)
   LET iterationError = sdEnergy(errorDecisionOrbit) - energyMem
   LET energyMem = sdEnergy(errorDecisionOrbit)
END SUB

!--- (2) set electron density rho <-- sdState(), occupation()
EXTERNAL SUB setElectronDensity
   FOR i=0 TO NNx-1
      LET rho(i) = rho(i)*(1.0-mixing)
      FOR ie=0 TO numberOfOrbit-1
         IF occupation(ie)>0.0 THEN
            LET rho(i) = rho(i) + mixing*occupation(ie)*(sdState(ie,i)*sdState(ie,i))/lylz
         END IF
      NEXT ie
   NEXT i
END SUB

!--- (3) set effective potential <-- electron density
EXTERNAL SUB setEffectivePotential
   DECLARE EXTERNAL SUB poisson,setVxc

   CALL poisson(20) !setVh
   CALL setVxc
   FOR i=0 TO NNx-1
      LET vv(i) = vvext(i)+vvh(i)+vvx(i)+vvc(i)
   NEXT i
END SUB

EXTERNAL SUB poisson(iterMax)
   LET h2 = 4.0*PI*dx*dx
   LET omegav2 = 0.5*1.8
   FOR iter=0 TO iterMax-1
      FOR i=1 TO NNx-2
         LET vvh(i) = vvh(i)+omegav2*(vvh(i+1)+vvh(i-1)-2.0*vvh(i) +h2*rho(i))
      NEXT i
   NEXT iter
END SUB

EXTERNAL SUB setVxc  !set exchage and correlation potential (Perdew and Zunger)
   LET c1 = -0.984745022
   FOR i=0 TO NNx-1
      LET rh = rho(i)
      LET rh3 = rh^0.33333333
      LET vvx(i) = c1*rh3
      LET rs = 0.6204/(rh3+1.0e-20)
      IF rs>=1.0 THEN
         LET sqrtrs = SQR(rs)
         LET ec = -0.1423/(1.0+1.0529*sqrtrs+0.3334*rs)
         LET vvc(i) = ec*(1.0+1.22838*sqrtrs+0.4445*rs)/(1.0+1.0529*sqrtrs+0.3334*rs)
      ELSE
         LET vvc(i) = -0.05837-0.0084*rs +(0.0311+0.00133*rs)*LOG(rs)
      END IF
   NEXT i
END SUB

EXTERNAL FUNCTION eeCorrelation(rh)
   LET r = 0.6204/(rh^0.33333333+1.0e-20)
   IF r>=1.0 THEN
      LET ec = -0.1423/(1.0+1.0529*SQR(r)+0.3334*r)
   ELSE
      LET ec = -0.0480-0.0116*r+(0.0311+0.0020*r)*LOG(r)
   END IF
   LET eeCorrelation = ec
END FUNCTION

!--- (4) solve Kohn-Sham equation
EXTERNAL SUB solveKohnSham(stateMax, iterMax)
   DECLARE EXTERNAL FUNCTION steepestDescent
   DECLARE EXTERNAL SUB GramSchmidt,sortState
   FOR i=0 TO iterMax-1
      FOR ist=0 TO stateMax-1
         LET sdEnergy(ist) = steepestDescent(ist,dampingFactor)
      NEXT ist
      CALL GramSchmidt(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>
      LET wrk(i) = (2*sdState(ist,i)-sdState(ist,i+1)-sdState(ist,i-1))/h2+(vv(i)-ei)*sdState(ist,i)
   NEXT i
   FOR i=1 TO NNx-2                          !---  |ist(next)> = |ist> - damp*|wrk>  ( norm(|ist(next)>) <>1 )
      LET sdState(ist,i) = sdState(ist,i)-damp*wrk(i)
   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
      LET s = s+sdState(ist,i)*((2*sdState(ist,i)-sdState(ist,i+1)-sdState(ist,i-1))/h2+vv(i)*sdState(ist,i))
      LET sn = sn + sdState(ist,i)*sdState(ist,i)
   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
            LET sdState(istate,i) = sdState(istate,i) - s*sdState(ist,i)
         NEXT i
      NEXT ist
      CALL normalizeState(istate)
   NEXT iState
END SUB

!--- (5) sort state
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
            LET w = sdState(ist,i)
            LET sdState(ist,i) = sdState(ist+1,i)
            LET sdState(ist+1,i) = w
         NEXT i
         LET w = sdEnergy(ist)
         LET sdEnergy(ist) = sdEnergy(ist+1)
         LET sdEnergy(ist+1) = w
      END IF
   NEXT ist
END SUB

!--- (6) set occupation
EXTERNAL SUB setOccupation(stateMax, nElectron)
   DECLARE EXTERNAL FUNCTION trialOcc,FermiDirac
   LET eUpper = sdEnergy(stateMax-1)+1.0
   LET eLower = sdEnergy(0)-1.0
   FOR i=0 TO stateMax-1
      IF sdEnergy(i)>eUpper THEN LET eUpper = sdEnergy(i)
      IF sdEnergy(i)<eLower THEN LET eLower = sdEnergy(i)
   NEXT i

   DO WHILE (eUpper-eLower>1.0e-12)
      LET eFermi = (eUpper+eLower)/2.0
      LET ntrial = trialOcc(stateMax, eFermi)
      IF ntrial<nElectron THEN
         LET eLower = eFermi
      ELSE
         LET eUpper = eFermi
      END IF
   LOOP
   LET eFermi = (eUpper+eLower)/2.0
   FOR i=0 TO stateMax-1
      LET occupation(i) = 2.0*FermiDirac(sdEnergy(i), eFermi)
      IF (occupation(i)<0.0001) THEN LET occupation(i) = 0.0
      IF (2.0-occupation(i)<0.0001) THEN LET occupation(i) = 2.0
   NEXT i
END SUB

EXTERNAL FUNCTION trialOcc(stateMax, eFermi)
   DECLARE EXTERNAL FUNCTION FermiDirac
   LET s = 0.0
   FOR i=0 TO stateMax-1
      LET s = s + 2.0*FermiDirac(sdEnergy(i), eFermi)
   NEXT i
   LET trialOcc = s
END FUNCTION

EXTERNAL FUNCTION FermiDirac(ee, ef)
   LET et = broadening !level width
   LET a = (ee-ef)/et
   IF a>100 THEN LET ret = 0.0 ELSE LET ret = 1.0/(EXP(a)+1.0)
   LET FermiDirac = ret
END FUNCTION

! ---------- utility

EXTERNAL FUNCTION innerProduct(ist,jst)  !--- <ist|jst>
   LET s = 0
   FOR i=1 TO NNx-2
      LET s = s + sdState(ist,i)*sdState(jst,i)
   NEXT i
   LET innerProduct = s*dx
END FUNCTION

EXTERNAL SUB normalizeState(ist)
   LET s = 0
   FOR i=1 TO NNx-2
      LET s = s + sdState(ist,i)*sdState(ist,i)*dx
   NEXT i
   LET a = SQR(1/s)
   FOR i=1 TO NNx-2
      LET sdState(ist,i) = a*sdState(ist,i)
   NEXT i
END SUB

EXTERNAL FUNCTION totalEnergy
   DECLARE EXTERNAL FUNCTION eeCorrelation
   LET sei = 0.0
   FOR i=0 TO numberOfOrbit-1
      LET sei = sei + occupation(i)*sdEnergy(i)
   NEXT i
   LET s = 0.0
   FOR i=1 TO NNx-1
      LET s = s + (-0.5*vvh(i)-0.25*vvx(i)+eeCorrelation(rho(i))-vvc(i))*rho(i)
   NEXT i
   LET s = s*dx
   LET totalEnergy = sei + s
END FUNCTION

! ---------- drawState

EXTERNAL SUB drawState(dispMode)  !public
   DECLARE EXTERNAL SUB dispInnerProduct,plotfn
   DECLARE EXTERNAL FUNCTION totalEnergy
   LET sc = 20
   LET xp = 50
   LET yp = 180
   LET vMag = 10
   LET stMag = 100

   SET DRAW MODE HIDDEN
   CLEAR
   SET LINE COLOR 1 ! black : PLOT x-axis
   PLOT LINES: xp,yp;dx*NNx*sc+xp,yp
   !---plot rho
   SET AREA COLOR 8 ! gray : PLOT rho;
   FOR i=0 TO NNx-2
      PLOT AREA: dx*i*sc+xp,yp;dx*i*sc+xp,rho(i)*20000+yp;dx*(i+1)*sc+xp,rho(i+1)*20000+yp;dx*(i+1)*sc+xp,yp
   NEXT i
   CALL plotFN(vvext,sc,vMag,xp,yp,0,1) !black, plot external potential vvext(x)
   !
   IF dispmode=2 THEN !---plot Vext,Veff(),vvh(),vvxc()x10
      SET TEXT HEIGHT 6
      CALL plotFN(vv,sc,vMag,xp,yp,0,10)    !dark green, plot effective potential vv(x)
      CALL plotFN(vvh,sc,vMag,xp,yp,0,2)    !blue, plot Hartree potential vvh(x)
      CALL plotFN(vvx,sc,vMag*10,xp,yp,0,4) !red, plot exchange potential vvx(x)
      CALL plotFN(vvc,sc,vMag*10,xp,yp,0,7) !magenta, plot correlation potential vvc(x)
      SET TEXT COLOR 1
      PLOT TEXT, AT 50, yp+65 :"Vext()"
      SET TEXT COLOR 10
      PLOT TEXT, AT 50, yp+50 :"Veff()"
      SET TEXT COLOR 2
      PLOT TEXT, AT 50, yp+35 :"Vh()"
      SET TEXT COLOR 4
      PLOT TEXT, AT 50, yp+20 :"Vx() x 10"
      SET TEXT COLOR 7
      PLOT TEXT, AT 50, yp+5 :"Vc() x 10"
   END IF
   IF dispMode=1 THEN !--- plot Vext(),|i>
      SET TEXT HEIGHT 5
      FOR ist=0 TO 4
         IF sdEnergy(ist)<12 THEN
            SET LINE COLOR 1+ist
            SET TEXT COLOR 1+ist
            FOR i=0 TO NNx-1 !plot wave function |psi(x,t)>
               PLOT LINES: i*dx*sc+xp, sdState(ist,i)*stMag+sdEnergy(ist)*20+yp;
            NEXT i
            PLOT LINES
            PLOT TEXT, AT xp-20,sdEnergy(ist)*20+yp :"|"&STR$(ist)&">"
         END IF
      NEXT ist
   END IF
   CALL dispInnerProduct(0,dx*NNx*sc+xp+20,yp)
   !--- caption
   SET TEXT HEIGHT 10
   SET TEXT COLOR 1 ! black
   PLOT TEXT, AT 50,115 ,USING "iteration count =###### error =#.##########":iterCount,iterationError
   PLOT TEXT, AT 50,100 ,USING "total energy =###.##########":totalEnergy
   PLOT TEXT, AT 50, 85 ,USING "E0 =###.########## Occ =#.#####":sdEnergy(0),occupation(0)
   PLOT TEXT, AT 50, 70 ,USING "E1 =###.########## Occ =#.#####":sdEnergy(1),occupation(1)
   PLOT TEXT, AT 50, 55 ,USING "E2 =###.########## Occ =#.#####":sdEnergy(2),occupation(2)
   PLOT TEXT, AT 50, 40 ,USING "E3 =###.########## Occ =#.#####":sdEnergy(3),occupation(3)
   PLOT TEXT, AT 50, 25 ,USING "E4 =###.########## Occ =#.#####":sdEnergy(4),occupation(4)
   PLOT TEXT, AT 50, 10 :"RS-DFT - Local Density Approximation 1D"
   PLOT TEXT, AT 50,470 :"[esc] exit            [D] chage dispMode"
   PLOT TEXT, AT 50,455 :"[1] ... [6] number of electron"
   PLOT TEXT, AT 50,440 :"[7] hermonics k*x^2   [8] quantum well"
   SET DRAW MODE EXPLICIT
END SUB

EXTERNAL SUB dispInnerProduct(ist,xp,yp)
   DECLARE EXTERNAL FUNCTION innerProduct
   SET TEXT HEIGHT 5
   SET TEXT COLOR 1 ! black
   FOR jst=0 TO numberOfOrbit-1
      PLOT TEXT, AT xp,yp+15*jst ,USING "("&STR$(ist)&"|"&STR$(jst)&") = -%.###^^^^":innerProduct(ist,jst)
   NEXT jst
   PLOT TEXT, AT xp,yp+15*10 :"(i|j) inner product"
END SUB

EXTERNAL SUB plotFN(a(),sc,mag,xp,yp,offset,col)
   SET LINE COLOR col
   FOR i=0 TO NNx-1
      PLOT LINES: dx*i*sc+xp,a(i)*mag+offset+yp;
   NEXT i
   PLOT LINES
END SUB

END MODULE
 

新1000桁モード

 投稿者:しばっち  投稿日:2017年 7月16日(日)16時31分22秒
返信・引用 編集済
  十進BASICの1000桁モードでは根性足りなかったので、新たな1000桁モードを設定して、FACT(1億)やSUPER-FACT(1000万)等
を計算してみました。(※ご注意 それなりに計算時間がかかります)
ttmathライブラリーを使用しています。

LET MODE=0
SELECT CASE MODE
CASE 0
   LET T=TIME
   LET X=FACT(12)
   CALL POW(X,X,RESULT$) !'X^N
   PRINT RESULT$
   PRINT TIME-T
   !2.16136854137681821022184125634148248819277439276856614943668779759374015927471583106592012734245952019394970811589706056678555509602666892262853286963399272624718192094412355351855126976834619525882212314189739932773415168630756169900513658510182323332373473499827943602913845866193736546382782945259951879659321847521838639531976483388098572301935336416235494360371625249430543171712585820018204557600757659419294850500417373875788626819825199950337907403249688891357217862083761274130969418499307206106124775263607305819648633066733354877683318544283412749937741320358307160031045916238726656414998645969668486063799136387061569535403249857358744409290597217065826574844652429571085685396226665454889018425548129950363668417167709598931948788021370222023373478957578453869374042172871420547662647531104912996278899652524570542599912981658810875761529791724632424261687434245565068590484867082473812801455874989871262294862658856398056257782919304106334636474909267649153798992665460574519126362677282548109387566698144e+4157895294
   ! 6.00000000122236E-2
CASE 1
   LET T=TIME
   LET N=12
   CALL POWFACT(N,RESULT$) !' (N!)^(N!)
   PRINT RESULT$
   PRINT TIME-T
   !2.16136854137681821022184125634148248819277439276856614943668779759374015927471583106592012734245952019394970811589706056678555509602666892262853286963399272624718192094412355351855126976834619525882212314189739932773415168630756169900513658510182323332373473499827943602913845866193736546382782945259951879659321847521838639531976483388098572301935336416235494360371625249430543171712585820018204557600757659419294850500417373875788626819825199950337907403249688891357217862083761274130969418499307206106124775263607305819648633066733354877683318544283412749937741320358307160031045916238726656414998645969668486063799136387061569535403249857358744409290597217065826574844652429571085685396226665454889018425548129950363668417167709598931948788021370222023373478957578453869374042172871420547662647531104912996278899652524570542599912981658810875761529791724632424261687434245565068590484867082473812801455874989871262294862658856398056257782919304106334636474909267649153798992665460574519126362677282548109387566698144e+4157895294
   ! 5.99999999976717E-2
CASE 2
   LET T=TIME
   LET N=100000000
   CALL LFACT(N,RESULT$) !' 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 *...* (N-1) * N
   PRINT RESULT$
   PRINT TIME-T
   !1.617203794921462386338773185612804043292374530648735079789813462666737154407504386580500876929190785730659306997709526379072095316295179654167179073513009006634129337747061591098685297097001582503472866301385800430032150945662267422861318875599783688513582747708092357961161833291289737062618637215019469899671053443152756646656487690303084660412647294491564580376647772129380120207037940015370921432087635413424927860657596349753458123900398747077870811819146068196799068285805276780342699291909785817326004701598706171485200668534925511351721008280251885258366413790982774377048289170918554487709474512123782332682445769981037301850337091283801728542366418127850072350755894754427590475945724598891919272598883437250343540668568318110852233598880090166164148355169172329626842127353660159845036832550429220759703016224818725548782786942271616990913661405300333388635357873739104632007105374412990250693891528936041921842151622873513865541985938245791021807671861353403161524245569590956406240525270162218772164532514879e+756570556
   !  682.120000000003
CASE 3
   LET T=TIME
   LET N=10
   CALL DOUBLEFACT(N,RESULT$) !' (N!)!
   PRINT RESULT$
   PRINT TIME-T
   !9.05199383547993459066223334541125889552198008664676868116314684779771059925518196453479041306538539527348189982329395778742235835651919605262124430846948778576764869559997043450588358061928072056124053527199058676605381825380218545252893543490674379874682902832717881991350611759327438286087082893716928362722869514941184566224960239661223087820620654988693997209045226577413230551325238368145865183939219774507733773526942850610284588755921757033509597723710534792930741070551439957695535284309143461400603472683366857009617354558034731296338086943875630076158484363480326013561877788402404321547742395651713816120381647441147623336006797943428300341263884257264797238025712931763921792106198924511782304221981368602255187251376620790571651695312481368830266300012667706495970524334606110206318956844192728556930489510830894015380706494724588299721349272612934424245289039282358777312680822791908369313735461594023948162097772384089698577582027209192522851246823823690027225358916556410761710770943207249755384549728357e+22228103
   ! 24.3600000000006
CASE 4
   LET T=TIME
   LET N=10000000
   CALL SUPERFACT(N,RESULT$) !' 2! * 3! * 4! * 5! * 6! * 7! *...* (N-1)! * N!
   PRINT RESULT$
   PRINT TIME-T
   !7.57878063697632117934715725034689887850671089443406030271039160509233584481021403298076663336844755155915921583136881963888696144823115311974475153652320052328567687659313122147525926979022065775096814149543780943831056005639617612529134106891342492205025312880072030110524448270901118158556779555922457929864716128158926304472044534907981834912866729470185803624587740859802958368016614391696499517054601914699582368295366201916833259354724597499035431744478719547535112048420085182135883183139640075048794059050861979118486547704710303906662116743881780796260260607839867444471378390091965608185269472213075704781769727546420860395320469358545090857923155193441929471103031866923296551225088390710078737054395326253048465352243419256199723124194271758472079381545351720056285172965339286261523760459249713701343864777816002666323040324461827576855665076282487404193204553774230151726249949958799269190261940238395540732979558437696555627226806719086098708609655790391738629544475963536580065298185018462037113821528992e+317427983505213
   ! 519.520000000004
CASE 5
   LET T=TIME
   LET N=1000000
   CALL HYPERFACT(N,RESULT$) !' 2^2 * 3^3 * 4^4 * 5^5 * 6^6 * 7^7 *...* N^N
   PRINT RESULT$
   PRINT TIME-T
   !6.23843686069241708059018283397497948396261032852709731523416295394277029282373398540713998575786733916012178463248223335823564645677499951494357410270839598310702765729112864006872878210039192370858067631070477998029933192482095612942476222728193754789852046900542637350103218384857392408531746283262373870131749000590618579432431865685268475514639336933206256949748969056305754856297026364320874076956431488052821375046045014638177390507897654960469730334797389190814704272105505552145212856118896574432939226143059057550973687604971378550419066289193097660009196503168010347848334618250277928124959562240318573539372896774337790426664723720446835857225634766789091996454216139135180235792568662383244534617370354298146218164699181831236845991726175424266619170479398316619116762592250453383101685685868627057862815641617259768873342240171953069022117764243910963990193443081142526456850279800179901219903311799202766853622229530523474152006391596501481823711359357236984255316292254520259892322621873560620869753255813e+2891429379524
   ! 1174.98
CASE 6
   LET T=TIME
   LET N=100
   CALL ULTRAFACT(N,RESULT$) !' 2!^2! * 3!^3! * 4!^4! * 5!^5! *...* (N!)^(N!) (※勝手に定義しました)
   PRINT RESULT$
   PRINT TIME-T
   !3.307080960022281488149685181352176531550030000680108534032717033505076494123338098391176644134271120923715544294326669107069219167043537761730096895592066843656679310779891205726619124126162841982779761492884571376859784533052299617920561168804808615757107870620134900273753181513620811433142712532260956074341430989241833002971570815290974543714913248254203286851480130951709685133299150786700397824107119153942159169775818240783061500294959931370412314159465834265700234915886098417613309965099149322128100572876169863206201780236823398411259753726291735981632356338054952252766969178844309737447458012208298750966505513188884915009461058844365777426551026604473623309791600908689872794148411741116573428111446864855182302393146522821144797551410563113777231534035768396111470722190923370955356992550575488022965531058889991735748065665507261010356836086850738831636682033967132801946334926843099424816607830075479676752197329292022559261874282671751989524980904001955577795260859806643024209758925103865456149835089471e+14889769765873800652572367200821121336307811992644550853047524406583971799601810165018835448400401070603264427270304269629350367510485185954426320887983405520084
   ! 1.95000000000437
CASE 7
   CALL LMAX(RESULT$)
   PRINT RESULT$
END SELECT
END

EXTERNAL  SUB POW(X,N,RESULT$)
LET B$=REPEAT$(CHR$(0),3000)
CALL POW1000(ABS(X),ABS(N),B$)
FOR I=LEN(B$) TO 1 STEP -1
   IF B$(I:I)<="9" AND B$(I:I)>="0" THEN
      EXIT FOR
   END IF
NEXT I
LET RESULT$=B$(1:I)

SUB POW1000(X,N,Y$)
   ASSIGN ".\DLL\ipow1000_128.dll","ipow1000" !'指数部 128bit 仮数部 3392bit(1000桁)
   !' ASSIGN ".\DLL\ipow1000_256.dll","ipow1000" !'指数部 256bit 仮数部 3392bit(1000桁)
END SUB
END SUB

EXTERNAL  SUB LFACT(N,RESULT$)
LET B$=REPEAT$(CHR$(0),3000)
CALL FACT1000(ABS(N),B$)
FOR I=LEN(B$) TO 1 STEP -1
   IF B$(I:I)<="9" AND B$(I:I)>="0" THEN
      EXIT FOR
   END IF
NEXT I
LET RESULT$=B$(1:I)

SUB FACT1000(N,Y$)
   ASSIGN ".\DLL\ipow1000_256.dll","fact1000"
END SUB
END SUB

EXTERNAL  SUB SUPERFACT(N,RESULT$)
LET B$=REPEAT$(CHR$(0),3000)
CALL SUPERFACT1000(ABS(N),B$)
FOR I=LEN(B$) TO 1 STEP -1
   IF B$(I:I)<="9" AND B$(I:I)>="0" THEN
      EXIT FOR
   END IF
NEXT I
LET RESULT$=B$(1:I)

SUB SUPERFACT1000(N,Y$)
   ASSIGN ".\DLL\ipow1000_512.dll","superfact1000"
END SUB
END SUB

EXTERNAL  SUB HYPERFACT(N,RESULT$)
LET B$=REPEAT$(CHR$(0),3000)
CALL HYPERFACT1000(ABS(N),B$)
FOR I=LEN(B$) TO 1 STEP -1
   IF B$(I:I)<="9" AND B$(I:I)>="0" THEN
      EXIT FOR
   END IF
NEXT I
LET RESULT$=B$(1:I)

SUB HYPERFACT1000(N,Y$)
   ASSIGN ".\DLL\ipow1000_512.dll","hyperfact1000"
END SUB
END SUB

EXTERNAL  SUB ULTRAFACT(N,RESULT$)
LET B$=REPEAT$(CHR$(0),3000)
CALL ULTRAFACT1000(ABS(N),B$)
FOR I=LEN(B$) TO 1 STEP -1
   IF B$(I:I)<="9" AND B$(I:I)>="0" THEN
      EXIT FOR
   END IF
NEXT I
LET RESULT$=B$(1:I)

SUB ULTRAFACT1000(N,Y$)
   ASSIGN ".\DLL\ipow1000_1024.dll","ultrafact1000"
END SUB
END SUB

EXTERNAL  SUB DOUBLEFACT(N,RESULT$)
LET B$=REPEAT$(CHR$(0),3000)
CALL DOUBLEFACT1000(N,B$)
FOR I=LEN(B$) TO 1 STEP -1
   IF B$(I:I)<="9" AND B$(I:I)>="0" THEN
      EXIT FOR
   END IF
NEXT I
LET RESULT$=B$(1:I)

SUB DOUBLEFACT1000(N,Y$)
   ASSIGN ".\DLL\ipow1000_128.dll","doublefact1000"
END SUB
END SUB

EXTERNAL  SUB POWFACT(N,RESULT$)
LET B$=REPEAT$(CHR$(0),3000)
CALL POWFACT1000(ABS(N),B$)
FOR I=LEN(B$) TO 1 STEP -1
   IF B$(I:I)<="9" AND B$(I:I)>="0" THEN
      EXIT FOR
   END IF
NEXT I
LET RESULT$=B$(1:I)

SUB POWFACT1000(N,Y$)
   ASSIGN ".\DLL\ipow1000_128.dll","powfact1000"
END SUB
END SUB

EXTERNAL  SUB LMAX(RESULT$)
LET B$=REPEAT$(CHR$(0),3000)
CALL MAX1000(B$)
FOR I=LEN(B$) TO 1 STEP -1
   IF B$(I:I)<="9" AND B$(I:I)>="0" THEN
      EXIT FOR
   END IF
NEXT I
LET RESULT$=B$(1:I)

SUB MAX1000(Y$)
   ASSIGN ".\DLL\ipow1000_3072.dll","max1000"
END SUB
END SUB
------------------------------------------------------------------------------------------
                                              ipow1000.cpp

#include <ttmath/ttmath.h>
#include <string>
#include <boost/lexical_cast.hpp>
using namespace std;
using boost::lexical_cast;

typedef ttmath::Big<TTMATH_BITS(EXP), TTMATH_BITS(3392)> bigfloat;
typedef ttmath::UInt<1024> bigint;

bigfloat ipow(unsigned int x,unsigned int n)
{
    bigfloat a,v;
    a=x;
    v=1;
    while(n!=0) {
        if((n % 2)==1) v=v*a;
        a=a*a;
        n=n>>1;
    }
    return v;
}

bigfloat ipow(bigint x,bigint n)
{
    bigfloat a,v;
    a=x;
    v=1;
    while(n!=0) {
        if((n % 2)==1) v=v*a;
        a=a*a;
        n=n>>1;
    }
    return v;
}

bigint fact(unsigned int n)
{
    bigint a;
    unsigned int i;
    a=1;
    for(i=2;i<=n;i++) a*=i;
    return a;
}

bigfloat fact(bigint n)
{
    bigfloat a;
    bigint i;
    a=1;
    for(i=2;i<=n;i++) a*=i;
    return a;
}

bigfloat superfact(unsigned int n)
{
    bigfloat a,b;
    a=1;
    b=1;
    for(unsigned int i=2; i<=n; i++)
    {
        a*=i;
        b*=a;
    }
    return b;
}

bigfloat hyperfact(unsigned int n)
{
    bigfloat a;
    a=1;
    for(unsigned int i=2; i<=n; i++) a*=ipow(i,i);
    return a;
}

bigfloat ultrafact(unsigned int n)
{
    bigfloat a;
    bigint x;
    a=1;
    for(unsigned int i=2; i<=n; i++)
    {
        x=fact(i);
        a*=ipow(x,x);
    }
    return a;
}

bigfloat powfact(unsigned int n)
{
    bigfloat a;
    bigint x;
    x=fact(n);
    a=ipow(x,x);
    return a;
}

bigfloat doublefact(unsigned int n)
{
    bigfloat a;
    a=fact(fact(n));
    return a;
}

bigfloat multifact(unsigned int n,unsigned int m)
{
    bigfloat v;
    v=1;
    for(unsigned int i=n; i>=1; i-=m) v*=i;
    return v;
}

extern "C"  __declspec(dllexport)  void ipow1000(unsigned int x,unsigned int n,char *b)
{
    bigfloat y;
    string s;
    y=ipow(x,n);
    s=boost::lexical_cast<string>(y);
    strcpy(b,s.c_str());
}

extern "C"  __declspec(dllexport)  void fact1000(unsigned int n,char *b)
{
    bigfloat y;
    string s;
    y=fact((bigint)n);
    s=boost::lexical_cast<string>(y);
    strcpy(b,s.c_str());
}

extern "C"  __declspec(dllexport)  void multifact1000(unsigned int n,unsigned int m,char *b)
{
    bigfloat y;
    string s;
    y=multifact(n,m);
    s=boost::lexical_cast<string>(y);
    strcpy(b,s.c_str());
}

extern "C"  __declspec(dllexport)  void superfact1000(unsigned int n,char *b)
{
    bigfloat y;
    string s;
    y=superfact(n);
    s=boost::lexical_cast<string>(y);
    strcpy(b,s.c_str());
}

extern "C"  __declspec(dllexport)  void hyperfact1000(unsigned int n,char *b)
{
    bigfloat y;
    string s;
    y=hyperfact(n);
    s=boost::lexical_cast<string>(y);
    strcpy(b,s.c_str());
}

extern "C"  __declspec(dllexport)  void ultrafact1000(unsigned int n,char *b)
{
    bigfloat y;
    string s;
    y=ultrafact(n);
    s=boost::lexical_cast<string>(y);
    strcpy(b,s.c_str());
}

extern "C"  __declspec(dllexport)  void powfact1000(unsigned int n,char *b)
{
    bigfloat y;
    string s;
    y=powfact(n);
    s=boost::lexical_cast<string>(y);
    strcpy(b,s.c_str());
}

extern "C"  __declspec(dllexport)  void doublefact1000(unsigned int n,char *b)
{
    bigfloat y;
    string s;
    y=doublefact(n);
    s=boost::lexical_cast<string>(y);
    strcpy(b,s.c_str());
}

extern "C"  __declspec(dllexport)  void max1000(char *b)
{
    bigfloat y;
    string s;
    y.SetMax();
    s=boost::lexical_cast<string>(y);
    strcpy(b,s.c_str());
}

VC++2012にてコンパイルしました。下記よりダウンロードしてください。(math.zip)

http://fast-uploader.com/file/7055745358295/

ダウンロードパス:shibacchi

有効期限は本日より1ヶ月です。


 

[021] 分子動力学プログラムの高速化

 投稿者:mike  投稿日:2017年 7月16日(日)10時21分52秒
返信・引用
   Morseポテンシャルによる分子動力学法(2次元)の粒子登録法による高速化プログラム[015]のさらなる高速化
プログラム(021fasterMMD2D.bas)を公開します。
 単純な粒子登録法はO(N^2)なので、他の部分がO(N)でも粒子数Nが大きくなるにつれて粒子登録法の部分が効いてきて、
やがて、実行時間はO(N^2)に近づきます。
 粒子登録の前に、全体を格子状の領域に区分し、粒子に近い領域から近い粒子を登録することにより、粒子登録をO(N)に
することができます。

実行時間の実測:
decimal BASIC 6.6.3.3 /MacOS 10.7.5/ corei7(2.7GHz)
fast   :従来の単純な粒子登録法(EXTERNAL SUB registerNearMolec)
faster :前処理+粒子登録法(EXTERNAL SUB preRegistration, EXTERNAL SUB registration)
      Molec   fast    (time/N)   faster  (time/N)
Fe   N= 535    10.02s (0.0187)    10.00s (0.0187)
Fe   N= 952    18.03s (0.0189)    17.01s (0.0179)
Fe   N=2139    48.00s (0.0224)    36.04s (0.0168)
Fe   N=3802   104.61s (0.0275)    63.99s (0.0168)
Fe   N=7187   272.17s (0.0379)   120.93s (0.0168)

前処理付きの粒子登録法(faster)では、ほぼ粒子数Nに比例した実行時間になっています。

試験環境:
本プログラムは十進BASIC 6.6.3.3 / macOS 10.7.5 でテストしました。


-------------

!
! ========= molecular dynamics 2D - Morse potential ==========
!
! 021fasterMMD2D.bas
!   Copyright(C) 2017 Mitsuru Ikeuchi
!   Released under the MIT license ( https://opensource.org/licenses/MIT )
!
! ver 0.0.1   2017.07.16  created
!
OPTION ARITHMETIC NATIVE
DECLARE EXTERNAL SUB mmd2d.setInitialCondition, mmd2d.moveParticles, mmd2d.drawParticles
DECLARE EXTERNAL FUNCTION INKEY$
LET tempMode = 0    !tempMode: 0:adiabatic  1:constant temperature
LET contTemp = 300  !contTemp: controled constant temperature(K)
LET drawMode = 1    !drawMode: 0:bond  1:circle+bond 2:velocitySpace


!setInitialCondition(molKind,boxSizeInNM,xtalSizeInNM,contTemp)
! decimal BASIC 6.6.2.2 /MacOS 10.7.5/ corei7(2.7GHz)  nMolec   fast     faster
!CALL setInitialCondition(3, 8.0, 6.0,contTemp) !3:Fe  N= 535    10.02s   10.00s
!CALL setInitialCondition(3,12.0, 8.0,contTemp) !3:Fe  N= 952    18.03s   17.01s
CALL setInitialCondition(3,16.0,12.0,contTemp)  !3:Fe  N=2139    48.00s   36.04s
!CALL setInitialCondition(3,20.0,16.0,contTemp) !3:Fe  N=3802   104.61s   63.99s
!CALL setInitialCondition(3,24.0,22.0,contTemp) !3:Fe  N=7187   272.17s  120.93s

INPUT  PROMPT "choice 0:fast 1:faster ":a$
IF VAL(a$)=1 THEN
   LET fasterSW=1
   PRINT "faster mode"
ELSE
   LET fasterSW=0
   PRINT "fast mode"
END IF

LET t0 = TIME
FOR it=1 TO 20
   CALL moveParticles(tempMode,contTemp,fasterSW)
   CALL drawParticles(tempMode,contTemp,drawMode)
NEXT it
PRINT TIME-t0
stop

END

! ---------- molecular dynamics 2D - Morse potential ----------
!
!  method: velocity Verlet ( F=m*d^2r/dt^2 -> r(t+dt)=r(t)+v*dt,v=v+(F/m)*dt )
!    (1) vi = vi + (Fi/mi)*(0.5dt)
!    (2) ri = ri + vi*dt
!    (3) calculation Fi <- {r1,r2,...,rn} Fi=sum(Fij,j=1 to n),Fij=F(ri-rj)
!    (4) vi = vi + (Fi/mi)*(0.5dt)
!    (6) goto (1)
!  potential: Morse V(r) = D*((1-EXP(-A*(r-r0)))^2-1)
!                        = D*(EXP(-2*A*(r-r0))-2*EXP(-A*(r-r0)))
!    (D:dissociation energy, r0:bond length, A:width parameter { A=SQR(k/(2*D)) }
!             force F(r) = -dV(r)/dr
!                        = 2*D*A*y*(y-1), y=EXP(-A*(r-r0))
MODULE mmd2d
MODULE OPTION ARITHMETIC NATIVE
PUBLIC SUB setInitialCondition !(molKind,boxSizeInNM,xtalSizeInNM,contTemp)
PUBLIC SUB moveParticles !(tempMode,contTemp)
PUBLIC SUB drawParticles !(drawMode:0:realSpace 1:velocitySpace)
SHARE NUMERIC molecKind,sysTime,dt,nMolec,xMax,yMax, Nsx,Nsy, mass,Dmrs,Amrs,r0mrs, rCutoff, hh
SHARE NUMERIC xx(8000),yy(8000)           ! (xx(i),yy(i))  : position of i-th particle
SHARE NUMERIC vx(8000),vy(8000)           ! (vx(i),vy(i))  : velocity of i-th particle
SHARE NUMERIC ffx(8000),ffy(8000)         ! (ffx(i),ffy(i)): total force of i-th particle
SHARE NUMERIC reg(8000,0 TO 100)         ! register near molec reg(i,0):number of near i-th molec
SHARE NUMERIC section(0 TO 101,0 TO 101,0 TO 20) !use pre-registration
SHARE NUMERIC molecData(0 TO 18,0 TO 3) ! molecule 0:mass, 1:epsilon, 2:sigma, 3:dt
SHARE NUMERIC forceTable(0 TO 1200)     ! force table
SHARE STRING molecSTR$(0 TO 18)         ! molec string
LET molecKind = 3          ! 3:Fe
LET sysTime = 0.0          ! system time (s) in the module
LET dt = 5.0*1.0e-15       ! time step (s)
LET nMolec = 100           ! number of particles
LET xMax = 6.0E-9          ! x-Box size (m)
LET yMax = 6.0E-9          ! y-Box size (m)
LET Nsx = 80
LET Nsy = 80
LET mass = 55.847*1.661e-27! mass of Fe (kg)
LET Dmrs = 0.4174*1.602e-19! D of Morse potential (J)   : energy of dissociation
LET Amrs = 1.3885e10       ! A of Morse potential (1/m) : width parameter
LET r0mrs = 2.845e-10      ! r0 of Morse potential (m)  : bond length
LET rCutoff = 1.0e-9       ! force cutoff radius (m)
LET hh = 1e-12             ! force table step (m)
LET molecData(0,0)=0       ! if molecData(0,0)=0 then set molecData(,)

! ---------- set initial condition

EXTERNAL SUB setInitialCondition(molKind,boxSizeInNM,xtalSizeInNM,contTemp)
   DECLARE EXTERNAL SUB setMolecData,setMolecules,ajustVelocity,setForceTable
   DECLARE EXTERNAL FUNCTION setCrystalBlock
   RANDOMIZE
   CALL setMolecData
   LET molecKind = molKind
   LET mass = molecData(molecKind,0)
   LET Dmrs = molecData(molecKind,1)
   LET Amrs = molecData(molecKind,2)
   LET r0mrs = molecData(molecKind,3)
   LET sysTime = 0.0
   LET xMax = boxSizeInNM*1.0e-9
   LET yMax = boxSizeInNM*1.0e-9
   ! set particles
   LET s = 0.5*(boxSizeInNM-xtalSizeInNM)*1.0e-9
   LET nMolec = setCrystalBlock(1, s, s, xtalSizeInNM*1.0e-9, xtalSizeInNM*1.0e-9, PI/4)
   CALL ajustVelocity(contTemp)
   LET rCutoff = MIN(1.0e-9, 3.0*r0mrs)
   CALL setForceTable
   ! set window
   SET WINDOW 0,500,0,500
END SUB

EXTERNAL SUB setMolecData
! Morse potential data
!       0:mass(AU) 1:D(eV)  2:A(m^-1)  3:r0(m)
   DATA   183.85 , 0.9906, 1.4116e10, 3.032e-10 !  0 W
   DATA    95.94 , 0.8032, 1.5079e10, 2.976e-10 !  1 Mo
   DATA    51.996, 0.4414, 1.5721e10, 2.754e-10 !  2 Cr
   DATA    55.847, 0.4174, 1.3885e10, 2.845e-10 !  3 Fe
   DATA    58.71 , 0.4205, 1.4199e10, 2.780e-10 !  4 Ni
   DATA    26.98 , 0.2703, 1.1646e10, 3.253e-10 !  5 Al
   DATA   207.19 , 0.2348, 1.1836e10, 3.733e-10 !  6 Pb
   DATA    63.54 , 0.3429, 1.3588e10, 2.866e-10 !  7 Cu
   DATA   107.87 , 0.3323, 1.3690e10, 3.115e-10 !  8 Ag
   DATA    39.948, 0.0104, 1.3400e10, 3.816e-10 !  9 Ar
   DATA   200.59 , 0.0734, 1.4900e10, 3.255e-10 ! 10 Hg
   DATA    40.08 , 0.1623, 0.8054e10, 4.569e-10 ! 11 Ca
   DATA    87.62 , 0.1513, 0.7878e10, 4.988e-10 ! 12 Sr
   DATA   137.34 , 0.1416, 0.6570e10, 5.373e-10 ! 13 Ba
   DATA    22.99 , 0.0633, 0.5900e10, 5.336e-10 ! 14 Na
   DATA    39.102, 0.0542, 0.4977e10, 6.369e-10 ! 15 K
   DATA    20.183, 0.0031, 1.6500e10, 3.076e-10 ! 16 Ne
   DATA    83.80 , 0.0141, 1.2500e10, 4.097e-10 ! 17 Kr
   DATA   131.30 , 0.0200, 1.2400e10, 4.467e-10 ! 18 Xe

   DATA "W" ,"Mo","Cr","Fe","Ni","Al","Pb","Cu","Ag","Ar","Hg"
   DATA "Ca","Sr","Ba","Na","K" ,"Ne","Kr","Xe"

   IF molecData(0,0)=0 THEN
      MAT READ molecData
      FOR i=0 TO 18
         LET molecData(i,0) = molecData(i,0)*1.661e-27 !mass(AU)--> (kg)
         LET molecData(i,1) = molecData(i,1)*1.602e-19 !eps(eV) --> (J)
      NEXT i
      MAT READ molecSTR$
   END IF
END SUB

EXTERNAL function setCrystalBlock(ii, x0, y0, xLen, yLen, theta)
   DECLARE EXTERNAL SUB setParticle
   LET iip = ii
   LET a = 0.98*r0mrs
   LET b = 0.866025*a
   LET leng = xLen
   IF (leng<yLen) THEN LET leng = yLen
   LET leng = 1.5*leng
   LET nx = INT(leng/b) + 1
   LET ny = INT(leng/a) + 1
   LET sth = SIN(theta)
   LET cth = COS(theta)
   FOR i=0 TO nx-1
      LET x = b*i - leng/2.0
      FOR j=0 TO ny-1
         LET y = a*j - leng/2.0
         IF MOD(i,2)=1 THEN LET y = y + 0.5*a
         LET xp = x0 + xLen/2.0 + cth*x - sth*y
         LET yp = y0 + yLen/2.0 + sth*x + cth*y
         IF (xp>=x0 AND xp<=x0+xLen AND yp>=y0 AND yp<=y0+yLen) THEN
            CALL setParticle(iip, xp, yp)
            LET iip = iip + 1
         END IF
      NEXT j
   NEXT i
   LET setCrystalBlock = iip - 1
end function

EXTERNAL SUB setParticle(i, x, y)
   LET xx(i) = x
   LET yy(i) = y
   LET vx(i) = 200.0*(RND+RND+RND+RND+RND+RND-3)
   LET vy(i) = 200.0*(RND+RND+RND+RND+RND+RND-3)
   LET ffx(i) = 0.0
   LET ffy(i) = 0.0
END SUB

EXTERNAL SUB setForceTable
   DECLARE EXTERNAL  FUNCTION cutoff
   FOR ir=10 TO 1200
      LET r = ir*hh
      LET y = EXP(-Amrs*(r-r0mrs))
      LET forceTable(ir) = cutoff(r)*2.0*Dmrs*Amrs*y*(y-1)
   NEXT ir
   FOR ir=0 TO 9
      LET forceTable(ir) = forceTable(10)
   NEXT ir
END SUB

EXTERNAL FUNCTION cutoff(r)
   IF r>0 AND r<0.8*rCutoff THEN
      LET ret = 1
   ELSEIF r>=0.8*rCutoff AND r<rCutoff THEN
      LET ret = 0.5+0.5*COS(PI*(r-0.8*rCutoff)/(0.2*rCutoff))
   else
      LET ret = 0
   END IF
   LET cutoff = ret
END FUNCTION

! ---------- move particles

EXTERNAL SUB moveParticles(tempMode,contTemp,fasterSW) !tempMode 0:adiabatic  1:constant-temp
   DECLARE EXTERNAL SUB moveParticlesDT,ajustVelocity,registerNearMolec,registration
   IF (tempMode=1) THEN CALL ajustVelocity(contTemp)
   IF fasterSW=1 THEN
      CALL registration !faster registration
   ELSE
      CALL registerNearMolec !fast registration
   END IF
   FOR i=1 TO 20
      CALL moveParticlesDT
   NEXT i
END SUB

EXTERNAL SUB moveParticlesDT ! velocity Verlet method
   DECLARE EXTERNAL SUB calcForce
   LET a = 0.5*dt/mass
   FOR i=1 TO nMolec
   !LET a = 0.5*dt/mass(i)
      LET vx(i) = vx(i)+a*ffx(i)
      LET vy(i) = vy(i)+a*ffy(i)
      LET xx(i) = xx(i)+vx(i)*dt
      LET yy(i) = yy(i)+vy(i)*dt
   NEXT i
   CALL calcForce
   FOR i=1 TO nMolec
   !LET a = 0.5*dt/mass(i)
      LET vx(i) = vx(i)+a*ffx(i)
      LET vy(i) = vy(i)+a*ffy(i)
   NEXT i
   LET sysTime=sysTime+dt
END SUB

EXTERNAL SUB calcForce
   DECLARE EXTERNAL FUNCTION force,boundaryForce
   LET s = 0.5*3.418e-10
   FOR i=1 TO nMolec
      LET ffx(i) = 0.0
      LET ffy(i) = 0.0
   NEXT i
   FOR i=1 TO nMolec-1
      FOR k=1 TO reg(i,0)-1
         LET j = reg(i,k)
         LET xij = xx(i)-xx(j)
         LET yij = yy(i)-yy(j)
         LET rij = SQR(xij*xij+yij*yij)
         IF rij<rCutoff THEN
            LET f = force(rij)
            LET fxij = f*xij/rij
            LET fyij = f*yij/rij
            LET ffx(i) = ffx(i)+fxij
            LET ffy(i) = ffy(i)+fyij
            LET ffx(j) = ffx(j)-fxij
            LET ffy(j) = ffy(j)-fyij
         END IF
      NEXT k
   NEXT i
   FOR i=1 TO nMolec  ! boundary force
      LET ffx(i) = ffx(i)+boundaryForce(xx(i)+s)+boundaryForce(xx(i)-xMax-s)
      LET ffy(i) = ffy(i)+boundaryForce(yy(i)+s)+boundaryForce(yy(i)-yMax-s)
   NEXT i
END SUB

!EXTERNAL FUNCTION force(r) ! force(r) = -dV(r)/dr
!   LET y = EXP(-Amrs*(r-r0mrs))
!   LET force = 2*Dmrs*Amrs*y*(y-1)
!END FUNCTION

EXTERNAL FUNCTION force(r) !force(r) <-- forceTable - linear interporation
   LET ir = INT(r/hh)
   LET a = r - ir*hh
   LET force = ((hh-a)*forceTable(ir) + a*forceTable(ir+1))/hh
END FUNCTION

EXTERNAL FUNCTION boundaryForce(r)
   LET r6 = (3.418e-10/r)^6
   LET boundaryForce = (24.0*(0.5*1.711e-21)*r6*(2.0*r6-1.0)/r)
END FUNCTION

EXTERNAL SUB registerNearMolec !fast registeration
   LET rCut = rCutoff+20*2000*dt
   LET rcut2 = rCut*rCut
   FOR i=1 TO nMolec-1
      LET k = 1
      FOR j=i+1 TO nMolec
         LET r2 = (xx(i)-xx(j))*(xx(i)-xx(j))+(yy(i)-yy(j))*(yy(i)-yy(j))
         IF (r2<rcut2) THEN
            LET reg(i,k) = j
            LET k = k + 1
         END IF
      NEXT j
      LET reg(i,0) = k
   NEXT i
END SUB

EXTERNAL FUNCTION maxNearMolec
   LET mx = 0
   FOR i=1 TO nMolec-1
      IF mx<reg(i,0) THEN LET mx = reg(i,0)
   NEXT i
   LET maxNearMolec = mx-1
END FUNCTION


EXTERNAL SUB registration !faster registration
   DECLARE EXTERNAL SUB preRegistration
   CALL preRegistration
   LET rreg = rCutoff+20*2000*dt
   LET rreg2 = rreg*rreg
   FOR ipp=1 TO nMolec-1
      LET kp = 1
      LET i0 = INT(Nsx*(xx(ipp)-rreg)/xMax)
      IF (i0<0) THEN LET i0 = 0
      LET i1 = INT(Nsx*(xx(ipp)+rreg)/xMax )
      IF (i1>=Nsx) THEN LET i1 = Nsx-1
      LET j0 = INT(Nsy*(yy(ipp)-rreg)/yMax )
      IF (j0<0) THEN LET j0 = 0
      LET j1 = INT(Nsy*(yy(ipp)+rreg)/yMax )
      IF (j1>=Nsy) THEN LET j1 = Nsy-1
      FOR i=i0 TO i1
         FOR j=j0 TO j1
            FOR iq=1 TO section(i,j,0)
               LET jp = section(i,j,iq)
               IF (jp>ipp) THEN
                  LET r2=(xx(ipp)-xx(jp))*(xx(ipp)-xx(jp))+(yy(ipp)-yy(jp))*(yy(ipp)-yy(jp))
                  IF (r2<rreg2) THEN
                     LET reg(ipp,kp) = jp
                     LET kp = kp + 1
                  END IF
               END IF
            NEXT iq
         NEXT j
      NEXT i
      LET reg(ipp,0) = kp
   NEXT ipp
END SUB

EXTERNAL SUB preRegistration
   FOR i=0 TO Nsx-1
      FOR j=0 TO Nsy-1
         LET section(i,j,0) = 0
      NEXT j
   NEXT i
   FOR ipp=1 TO nMolec
      LET i = INT(Nsx*xx(ipp)/xMax)
      IF i>=Nsx THEN LET i = Nsx-1
      LET j = INT(Nsy*yy(ipp)/yMax)
      IF j>=Nsy THEN LET j = Nsy-1
      LET iq = section(i,j,0) + 1
      LET section(i,j,0) = iq
      LET section(i,j,iq) = ipp
   NEXT ipp
END SUB

! ---------- utility

EXTERNAL FUNCTION systemTemprature
   LET kB = 1.38e-23 ! Boltzman's constant (J/K)
   LET ek= 0.0       !kinetic energy (J)
   FOR i=1 TO nMolec
      LET ek = ek + 0.5*mass*(vx(i)^2+vy(i)^2)
   NEXT i
   LET systemTemprature = ek/(nMolec*kB) !2D: E/N=kT, 3D: E/N=(3/2)kT
END FUNCTION

EXTERNAL SUB ajustVelocity(temp)
   DECLARE EXTERNAL FUNCTION systemTemprature
   LET r = sqr(temp/systemTemprature)
   FOR i=1 TO nMolec
      LET vx(i) = r*vx(i)
      LET vy(i) = r*vy(i)
   NEXT i
END SUB

! ---------- draw particles

EXTERNAL SUB drawParticles(tempMode,contTemp,drawMode)
   DECLARE EXTERNAL FUNCTION systemTemprature,maxNearMolec
   DECLARE EXTERNAL sub realSpace,velocitySpace,plotBond
   SET DRAW MODE HIDDEN
   CLEAR
   IF drawMode=0 OR drawMode=1 THEN !--- 0:circle  1:circle+bond
      call plotBond(drawMode)
   ELSEIF drawMode=2 THEN
      call velocitySpace
   END IF

   !--- draw caption
   SET TEXT HEIGHT 10
   SET TEXT COLOR 1 ! black
   LET tmp$ = "adiabatic   constantTemp  "
   PLOT TEXT, AT  50, 70 ,USING "time =#####.## (ps)   temp =####.## (K)":sysTime*1E12,systemTemprature
   PLOT TEXT, AT  50, 55 ,USING molecSTR$(molecKind)&"  N =####":nMolec
   PLOT TEXT, AT  50, 40 ,USING "tempMode = ## ":tempMode
   PLOT TEXT, AT 200, 40 :tmp$(tempMode*12+1:tempMode*12+12)
   PLOT TEXT, AT  50, 25 ,USING "controled Temperature =####.# (K)":contTemp
   PLOT TEXT, AT  50, 10 :"2-dimensional molecular dynamics"
   SET DRAW MODE EXPLICIT
END SUB

EXTERNAL sub plotBond(drawMode)
   LET boxSize = 300
   LET mag = boxSize/xMax
   LET xp = 100
   LET yp = 100
   SET LINE COLOR 1 ! black : !--- box
   PLOT LINES: xp,yp; boxSize+xp,yp; boxSize+xp,boxSize+yp; xp,boxSize+yp; xp,yp
   SET TEXT HEIGHT 6
   SET TEXT COLOR 1 ! black
   PLOT TEXT, AT xp,boxSize+2+yp ,USING  "box size =##.# x ##.# (nm)":xMax*1e9,yMax*1e9
   FOR i=1 TO nMolec
      SET LINE COLOR 8 ! gray
      DRAW circle WITH SCALE(r0mrs/2.0*mag)*SHIFT(xx(i)*mag+xp,yy(i)*mag+yp)
   NEXT i
   IF drawMode=1 THEN
      FOR i=1 TO nMolec-1
         FOR k=1 TO reg(i,0)-1
            LET j = reg(i,k)
            LET r = SQR((xx(i)-xx(j))*(xx(i)-xx(j))+(yy(i)-yy(j))*(yy(i)-yy(j)))
            IF r<1.2*r0mrs THEN
               IF r<0.93*r0mrs THEN
                  SET LINE COLOR 4 !red 13 !'oleave
               ELSEIF r<1.03*r0mrs THEN
                  SET LINE COLOR 3 !green
               ELSEIF r<1.08*r0mrs THEN
                  SET LINE COLOR 2 !blue
               ELSE
                  SET LINE COLOR 8 !gray
               END IF
               PLOT LINES: xx(i)*mag+xp,yy(i)*mag+yp;xx(j)*mag+xp,yy(j)*mag+yp
            END IF
         NEXT k
      NEXT i
   END IF
END sub

EXTERNAL sub velocitySpace
   LET boxSize = 300
   LET xp = 100
   LET yp = 100
   SET LINE COLOR 1 !black : axis
   PLOT LINES: xp,boxSize/2+yp; boxSize+xp,boxSize/2+yp !vx-axis
   PLOT LINES: boxSize/2+xp,yp; boxSize/2+xp,boxSize+yp !vy-axis
   SET TEXT HEIGHT 6
   SET TEXT COLOR 1 ! black
   PLOT TEXT, AT boxSize+xp,boxSize/2+yp: "vx"
   PLOT TEXT, AT boxSize+xp,boxSize/2-12+yp: "2000m/s"
   PLOT TEXT, AT boxSize/2-12+xp,boxSize+yp: "vy  2000m/s"
   PLOT TEXT, AT boxSize/2-8+xp,boxSize/2-10+yp: "0"
   PLOT TEXT, AT xp,boxSize+8+yp: "velocity space (vx,vy)"
   LET mag = boxSize/4000
   FOR i=1 TO nMolec
      SET LINE COLOR 2 ! blue
      DRAW circle WITH SCALE(5)*SHIFT(vx(i)*mag+boxSize/2+xp,vy(i)*mag+boxSize/2+yp)
   NEXT i
END sub

END MODULE
 

[020]周期的境界条件下の電子状態(3次元)

 投稿者:mike  投稿日:2017年 7月10日(月)09時40分29秒
返信・引用
   3次元-周期的境界条件下の下で定常状態の電子のエネルギーと電子状態を求めるプログラム
(020harmonicsPSD2D.bas)を公開します。これは、最急降下法による3次元の定常状態の電子状態を求めるプログラム
[013]に周期的境界条件を付けたものです。

表示の説明:
 3次元-周期的境界条件下の自由空間(V(x,y,z)=0)では基底状態は|0>=constant, E(0)=0 になります。
x、y、z方向の周期が同じ場合、E1=E2=E3=E4=E5=E6で縮退しているため、|1>..|6>はこの固有空間の中の
任意の直交する関数の組でよいため、形が一意には決まりません。

試験環境:
本プログラムは十進BASIC 6.6.3.3 / macOS 10.7.5, 十進BASIC Ver 7.8.0 / windows 10 でテストしました。


-----------

!
! ========= periodic steepest descent method 3D ==========
!
! 020periodicPSD3D.bas
!   Copyright(C) 2017 Mitsuru Ikeuchi
!   Released under the MIT license ( https://opensource.org/licenses/MIT )
!
! ver 0.0.1   2017.07.10  created
!
OPTION ARITHMETIC NATIVE
DECLARE EXTERNAL SUB psd3d.setInitialCondition,psd3d.SDiteration,psd3d.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 = 3     !vIndex=0: V=1/r, 1: V=(1/2)r^2, 2:V=0(r<=4) or V=16(r>4) 3:V=0
LET stateMax = 6   !state 0,1,...,stateMax-1
LET drawMode = 0   !drawMode=0:density3D,  1:grid |i> in (x,y,0)
LET ist = 0        !display state |i>
CALL setInitialCondition(stateMax,vIndex)

DO
   LET Ay = Ay + ddAy
   CALL SDiteration(stateMax, 2) !2 - iteration in SDiteration()
   CALL drawState(ist,Ax,Ay,drawMode)
   LET S$=INKEY$
   IF S$="0" OR S$="1" OR S$="2" OR S$="3"  OR S$="4" OR S$="5" THEN
      LET ist = VAL(S$)
   ELSEIF S$="W" OR S$="w" THEN
      if ddAy=0.0 then LET ddAy=-pi/180 else LET ddAy = 0.0
   ELSEIF S$="E" OR S$="e" THEN
      IF ddAy=0.0 THEN LET ddAy= PI/180 ELSE LET ddAy = 0.0
   ELSEIF S$="D" OR S$="d" THEN
      LET drawMode = MOD(drawMode+1,2)
   END IF
LOOP UNTIL S$=CHR$(27)
END

EXTERNAL FUNCTION INKEY$ !from decimal BASIC library
OPTION ARITHMETIC NATIVE
SET ECHO "OFF"
LET S$=""
CHARACTER INPUT NOWAIT: S$
LET INKEY$=S$
END FUNCTION

! ---------- periodic steepest descent method 3D ----------
!
! 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 psd3d
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 = 16                         ! x-max number of sdState(,NNx,NNy,NNz)
LET NNy = 16                         ! y-max number of sdState(,NNx,NNy,NNz)
LET NNz = 16                         ! 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
   CALL setInitialState(stateMax)
   CALL setPotential(vIndex)
   CALL setBox
   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=0 TO NNx-1
         FOR j=0 TO NNy-1
            FOR k=0 TO NNz-1
               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
            IF vIndex=0 THEN
               LET v0 = 10.0
               IF SQR((x-8)*(x-8)+(z-8)*(z-8))<2 THEN LET V0 = 0.0
               LET vv(i,j,k) = v0
            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=0 TO NNx-1                          !---  |wrk> = (H-E_ist)|ist>
      LET ip1 = MOD(i+1,NNx)
      LET im1 = MOD(i-1+NNx,NNx)
      FOR j=0 TO NNy-1
         LET jp1 = MOD(j+1,NNy)
         LET jm1 = MOD(j-1+NNy,NNy)
         FOR k=0 TO NNz-1
            LET kp1 = MOD(k+1,NNz)
            LET km1 = MOD(k-1+NNz,NNz)
            LET kesdState = (6*sdState(ist,i,j,k)-sdState(ist,ip1,j,k)-sdState(ist,im1,j,k)&
&-sdState(ist,i,jp1,k)-sdState(ist,i,jm1,k)-sdState(ist,i,j,kp1)-sdState(ist,i,j,km1))/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=0 TO NNx-1                          !---  |ist(next)> = |ist> - damp*|wrk>  ( norm(|ist(next)>) <>1 )
      FOR j=0 TO NNy-1
         FOR k=0 TO NNz-1
            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=0 TO NNx-1                          !---  |wrk> = (H-E_ist)|ist>
      LET ip1 = MOD(i+1,NNx)
      LET im1 = MOD(i-1+NNx,NNx)
      FOR j=0 TO NNy-1
         LET jp1 = MOD(j+1,NNy)
         LET jm1 = MOD(j-1+NNy,NNy)
         FOR k=0 TO NNz-1
            LET kp1 = MOD(k+1,NNz)
            LET km1 = MOD(k-1+NNz,NNz)
            LET kesdState = (6*sdState(ist,i,j,k)-sdState(ist,ip1,j,k)-sdState(ist,im1,j,k)&
&-sdState(ist,i,jp1,k)-sdState(ist,i,jm1,k)-sdState(ist,i,j,kp1)-sdState(ist,i,j,km1))/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=0 TO NNx-1
            FOR j=0 TO NNy-1
               FOR k=0 TO NNz-1
                  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=0 TO NNx-1
      FOR j=0 TO NNy-1
         FOR k=0 TO NNz-1
            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=0 TO NNx-1
      FOR j=0 TO NNy-1
         FOR k=0 TO NNz-1
            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=0 TO NNx-1
      FOR j=0 TO NNy-1
         FOR k=0 TO NNz-1
            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)
   DECLARE EXTERNAL SUB drawDensity3D,setRotateXY,drawDensity3D,drawStateGrid
   SET DRAW MODE HIDDEN
   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 :"periodic steepest descent method 3D"
   PLOT TEXT, AT 50,470 :"[0],[1],...,[5] view state   [esc] exit "
   PLOT TEXT, AT 50,455 :"[D] change drawMode"
   PLOT TEXT, AT 50,440 :"[W] rotateLeft/stop [E] rotateRight/stop"
   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 ii=0 TO NNx-1
      LET i=ii
      IF (pzApex(1)-pzApex(0)<0) THEN LET i=NNx-ii-1
      FOR jj=0 TO NNy-1
         LET j=jj
         IF (pzApex(2)-pzApex(0)<0) THEN LET j=NNy-jj-1
         FOR kk=0 TO NNz-1
            LET k=kk
            IF (pzApex(4)-pzApex(0)<0) THEN LET k=NNz-kk-1
            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 kk
      NEXT jj
   NEXT ii
   !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)":MOD(xRotateAngle*180/PI,360),MOD(yRotateAngle*180/PI,360)
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
 

[019]周期的境界条件下の分子動力学

 投稿者:mike  投稿日:2017年 7月 2日(日)05時11分40秒
返信・引用
  周期的境界条件の付いた分子動力学法プログラム(019periodicMMD2D.bas)を公開します。
周期的境界条件は、結晶をあつかう場合、表面の影響がないようにするため有効な方法です。
ここでは[017]のプログラムに周期的境界条件を付けました。
(物理的には、あまり意味がありません。プログラムの雛形です。)

表示の説明:
鉄原子は暗い水色、アルミニウム原子は茶色の円であらわします。
表示は[1]のキーを押すことで変更できます。

試験環境:
本プログラムは十進BASIC 6.6.3.3 / macOS 10.7.5, 十進BASIC Ver 7.8.0 / windows 10 でテストしました。

----------------

!
! ========= molecular dynamics 2D - Morse potential ==========
!
! 019periodicMMD2D.bas
!   Copyright(C) 2017 Mitsuru Ikeuchi
!   Released under the MIT license ( https://opensource.org/licenses/MIT )
!
! ver 0.0.1   2017.07.01  created
!
OPTION ARITHMETIC NATIVE
DECLARE EXTERNAL SUB mmd2d.setInitialCondition, mmd2d.moveParticles, mmd2d.drawParticles
DECLARE EXTERNAL FUNCTION INKEY$
LET tempMode = 0    !tempMode: 0:adiabatic  1:constant temperature
LET contTemp = 300  !contTemp: controled constant temperature(K)
LET drawMode = 1    !drawMode: 0:bond  1:circle+bond 2:velocitySpace
DIM menu$(8)
LET menu$(1) = "Fe Al"
LET menu$(2) = "Fe Mo"
LET menu$(3) = "Fe Hg"
LET menu$(4) = "W Pb"
LET menu$(5) = "Ni Cr"
LET menu$(6) = "Ag Cu"
LET menu$(7) = "Al Ar"
LET menu$(8) = "continue"

!setInitialCondition(kind1,kind2,boxSizeInNM,contTemp)
!molKind: 0:W 1:Mo 2:Cr 3:Fe 4:Ni 5:Al 6:Pb 7:Cu 8:Ag 9:Ar 10:Hg
CALL setInitialCondition(3,5,6.0,contTemp)

DO
   CALL moveParticles(tempMode,contTemp)
   CALL drawParticles(tempMode,contTemp,drawMode)
   LET S$=INKEY$
   IF S$="0" THEN
      LOCATE VALUE (1),RANGE 0 TO 4000 : temp
      IF temp>1 THEN
         LET tempMode = 1
         LET contTemp = temp
      ELSE
         LET tempMode = 0
      END IF
   ELSEIF S$="1" THEN
      LET drawMode = MOD(drawMode+1,3)
   ELSEIF S$="9" THEN
      LOCATE CHOICE (menu$) :nmenu
      IF nmenu=1 THEN      !--- Fe Al box=6x6nm
         CALL setInitialCondition(3,5,6.0,contTemp)
      ELSEIF nmenu=2 THEN  !--- Fe Mo
         CALL setInitialCondition(3,2,6.0,contTemp)
      ELSEIF nmenu=3 THEN  !--- Fe Hg
         CALL setInitialCondition(3,20,6.0,contTemp)
      ELSEIF nmenu=4 THEN  !--- W Pb
         CALL setInitialCondition(0,6,6.0,contTemp)
      ELSEIF nmenu=5 THEN  !--- Ni Cr
         CALL setInitialCondition(4,2,6.0,contTemp)
      ELSEIF nmenu=6 THEN  !--- Ag Cu
         CALL setInitialCondition(8,7,6.0,contTemp)
      ELSEIF nmenu=7 THEN  !--- Al Ar
         CALL setInitialCondition(5,17,6.0,contTemp)
      ELSEIF nmenu=8 THEN  !contine
      !
      END IF

   END IF
LOOP UNTIL S$=CHR$(27)

END

EXTERNAL FUNCTION INKEY$  !from decimal BASIC library
OPTION ARITHMETIC NATIVE
SET ECHO "OFF"
LET S$=""
CHARACTER INPUT NOWAIT: S$
LET INKEY$=S$
END FUNCTION

! ---------- molecular dynamics 2D - Morse potential ----------
!
!  method: velocity Verlet ( F=m*d^2r/dt^2 -> r(t+dt)=r(t)+v*dt,v=v+(F/m)*dt )
!    (1) vi = vi + (Fi/mi)*(0.5dt)
!    (2) ri = ri + vi*dt
!    (3) calculation Fi <- {r1,r2,...,rn} Fi=sum(Fij,j=1 to n),Fij=F(ri-rj)
!    (4) vi = vi + (Fi/mi)*(0.5dt)
!    (6) goto (1)
!  potential: Morse V(r) = D*((1-EXP(-A*(r-r0)))^2-1)
!                        = D*(EXP(-2*A*(r-r0))-2*EXP(-A*(r-r0)))
!    (D:dissociation energy, r0:bond length, A:width parameter { A=SQR(k/(2*D)) }
!             force F(r) = -dV(r)/dr
!                        = 2*D*A*y*(y-1), y=EXP(-A*(r-r0))
MODULE mmd2d
MODULE OPTION ARITHMETIC NATIVE
PUBLIC SUB setInitialCondition !(molKind,boxSizeInNM,xtalSizeInNM,contTemp)
PUBLIC SUB moveParticles !(tempMode,contTemp)
PUBLIC SUB drawParticles !(drawMode:0:realSpace 1:velocitySpace)
SHARE NUMERIC molecKind1,molecKind2, sysTime, dt, nMolec, xMax, yMax, rCutoff, hh
SHARE NUMERIC xx(1000),yy(1000)           ! (xx(i),yy(i))  : position of i-th particle
SHARE NUMERIC vx(1000),vy(1000)           ! (vx(i),vy(i))  : velocity of i-th particle
SHARE NUMERIC ffx(1000),ffy(1000)         ! (ffx(i),ffy(i)): total force of i-th particle
SHARE NUMERIC kind(1000),mas(1000)        ! kind(i),mas(i) : molec kind, mass of i-th particle
SHARE NUMERIC reg(1000,0 TO 100)         ! register near molec reg(i,0):number of near i-th molec
SHARE NUMERIC molecData(0 TO 20,0 TO 3) ! molecule 0:mass, 1:epsilon, 2:sigma, 3:dt
SHARE STRING molecSTR$(0 TO 20)         ! molec string
SHARE NUMERIC forceTable(0 TO 20, 0 TO 20,0 TO 1001) ! force table

LET molecKind1 = 3         ! molecule kind1 - 3:Fe
LET molecKind2 = 1         ! molecule kind2 - 1:Mo
LET sysTime = 0.0          ! system time (s) in the module
LET dt = 5.0*1.0e-15       ! time step (s)
LET nMolec = 100           ! number of particles
LET xMax = 6.0E-9          ! x-Box size (m)
LET yMax = 6.0E-9          ! y-Box size (m)
LET rCutoff = 1.0e-9       ! force cutoff radius (m)
LET hh = 1e-12             ! force table step (m)
LET molecData(0,0)=0       ! if molecData(0,0)=0 then set molecData(,)

! ---------- set initial condition

EXTERNAL SUB setInitialCondition(kind1,kind2,boxSizeInNM,contTemp)
   DECLARE EXTERNAL SUB setMolecData,ajustVelocity,setForceTable
   DECLARE EXTERNAL FUNCTION setCrystalBlock
   RANDOMIZE
   CALL setMolecData
   LET molecKind1 = kind1
   LET molecKind2 = kind2
   LET sysTime = 0.0
   LET xMax = boxSizeInNM*1.0e-9
   LET yMax = boxSizeInNM*1.0e-9
   ! set particles
   LET s = 0.07*xMax
   LET xtalSize = 0.4*xMax
   LET ss = 0.53*xMax
   LET ii = setCrystalBlock(1, kind1, s, s, xtalSize, xtalSize, PI/4)
   LET ii = setCrystalBlock(ii+1, kind2, s, ss, xtalSize, xtalSize, 0)
   LET ii = setCrystalBlock(ii+1, kind2, ss, s, xtalSize, xtalSize, 0)
   LET nMolec = setCrystalBlock(ii+1, kind1, ss, ss, xtalSize, xtalSize, 0)
   CALL ajustVelocity(contTemp)
   LET rCutoff = MIN(1.0e-9, 3.0*MAX(molecData(kind1,3),molecData(kind2,3)))
   CALL setForceTable
   ! set window
   SET WINDOW 0,500,0,500
END SUB

EXTERNAL SUB setMolecData
! Morse potential data
!       0:mass(AU) 1:D(eV)  2:A(m^-1)  3:r0(m)
   DATA   183.85 , 0.9906, 1.4116e10, 3.032e-10 !  0 W
   DATA    95.94 , 0.8032, 1.5079e10, 2.976e-10 !  1 Mo
   DATA    51.996, 0.4414, 1.5721e10, 2.754e-10 !  2 Cr
   DATA    55.847, 0.4174, 1.3885e10, 2.845e-10 !  3 Fe
   DATA    58.71 , 0.4205, 1.4199e10, 2.780e-10 !  4 Ni
   DATA    26.98 , 0.2703, 1.1646e10, 3.253e-10 !  5 Al
   DATA   207.19 , 0.2348, 1.1836e10, 3.733e-10 !  6 Pb
   DATA    63.54 , 0.3429, 1.3588e10, 2.866e-10 !  7 Cu
   DATA   107.87 , 0.3323, 1.3690e10, 3.115e-10 !  8 Ag
   DATA    40.08 , 0.1623, 0.8054e10, 4.569e-10 !  9 Ca
   DATA    87.62 , 0.1513, 0.7878e10, 4.988e-10 ! 10 Sr
   DATA   137.34 , 0.1416, 0.6570e10, 5.373e-10 ! 11 Ba
   DATA    22.99 , 0.0633, 0.5900e10, 5.336e-10 ! 12 Na
   DATA    39.102, 0.0542, 0.4977e10, 6.369e-10 ! 13 K
   DATA    85.47 , 0.0464, 0.4298e10, 7.207e-10 ! 14 Rb
   DATA   132.905, 0.0449, 0.4157e10, 7.557e-10 ! 15 Cs
   DATA    20.183, 0.0031, 1.6500e10, 3.076e-10 ! 16 Ne
   DATA    39.948, 0.0104, 1.3400e10, 3.816e-10 ! 17 Ar
   DATA    83.80 , 0.0141, 1.2500e10, 4.097e-10 ! 18 Kr
   DATA   131.30 , 0.0200, 1.2400e10, 4.467e-10 ! 19 Xe
   DATA   200.59 , 0.0734, 1.4900e10, 3.255e-10 ! 20 Hg

   DATA "W" ,"Mo","Cr","Fe","Ni","Al","Pb","Cu","Ag","Ca","Sr"
   DATA "Ba","Na","K" ,"Rb","Cs","Ne","Ar","Kr","Xe","Hg"

   IF molecData(0,0)=0 THEN
      MAT READ molecData
      FOR i=0 TO 20
         LET molecData(i,0) = molecData(i,0)*1.661e-27 !mass(AU)--> (kg)
         LET molecData(i,1) = molecData(i,1)*1.602e-19 !eps(eV) --> (J)
      NEXT i
      MAT READ molecSTR$
   END IF
END SUB

EXTERNAL FUNCTION setCrystalBlock(ii, knd, x0, y0, xLen, yLen, theta)
   DECLARE EXTERNAL SUB setParticle
   LET iip = ii
   LET a = 0.98*molecData(knd,3) !r0 of knd
   LET b = 0.866025*a
   LET leng = xLen
   IF (leng<yLen) THEN LET leng = yLen
   LET leng = 1.5*leng
   LET nx = INT(leng/b) + 1
   LET ny = INT(leng/a) + 1
   LET sth = SIN(theta)
   LET cth = COS(theta)
   FOR i=0 TO nx-1
      LET x = b*i - leng/2.0
      FOR j=0 TO ny-1
         LET y = a*j - leng/2.0
         IF MOD(i,2)=1 THEN LET y = y + 0.5*a
         LET xp = x0 + xLen/2.0 + cth*x - sth*y
         LET yp = y0 + yLen/2.0 + sth*x + cth*y
         IF (xp>=x0 AND xp<=x0+xLen AND yp>=y0 AND yp<=y0+yLen) THEN
            CALL setParticle(iip, knd, xp, yp)
            LET iip = iip + 1
         END IF
      NEXT j
   NEXT i
   LET setCrystalBlock = iip - 1
end function

EXTERNAL SUB setParticle(i, knd, x, y)
   LET xx(i) = x
   LET yy(i) = y
   LET vx(i) = 200.0*(RND+RND+RND+RND+RND+RND-3)
   LET vy(i) = 200.0*(RND+RND+RND+RND+RND+RND-3)
   LET ffx(i) = 0.0
   LET ffy(i) = 0.0
   LET mas(i) = molecData(knd,0)
   LET kind(i) = knd
END SUB

EXTERNAL SUB setForceTable
   DECLARE EXTERNAL  FUNCTION cutoff
   FOR ki=0 TO 20
      FOR kj=0 TO 20
         LET dd = SQR(molecData(ki,1)*molecData(kj,1))
         LET aa = 0.5*(molecData(ki,2)+molecData(kj,2))
         LET r0 = 0.5*(molecData(ki,3)+molecData(kj,3))
         FOR ir=10 TO 1001
            LET r = ir*hh
            LET y = EXP(-aa*(r-r0))
            LET forceTable(ki,kj,ir) = cutoff(r)*2.0*dd*aa*y*(y-1)
         NEXT ir
         FOR ir=0 TO 9
            LET forceTable(ki,kj,ir) = forceTable(ki,kj,10)
         NEXT ir
      NEXT kj
   NEXT ki
END SUB

EXTERNAL FUNCTION cutoff(r)
   IF r>0 AND r<0.8*rCutoff THEN
      LET ret = 1
   ELSEIF r>=0.8*rCutoff AND r<rCutoff THEN
      LET ret = 0.5+0.5*COS(PI*(r-0.8*rCutoff)/(0.2*rCutoff))
   else
      LET ret = 0
   END IF
   LET cutoff = ret
END FUNCTION

! ---------- move particles

EXTERNAL SUB moveParticles(tempMode,contTemp) !tempMode 0:adiabatic  1:constant-temp
   DECLARE EXTERNAL SUB moveParticlesDT,ajustVelocity,registerNearMolec
   IF (tempMode=1) THEN CALL ajustVelocity(contTemp)
   CALL registerNearMolec
   FOR i=1 TO 20
      CALL moveParticlesDT
   NEXT i
END SUB

EXTERNAL SUB moveParticlesDT ! velocity Verlet method
   DECLARE EXTERNAL SUB calcForce
   FOR i=1 TO nMolec
      LET a = 0.5*dt/mas(i)
      LET vx(i) = vx(i)+a*ffx(i)
      LET vy(i) = vy(i)+a*ffy(i)
      LET xx(i) = xx(i)+vx(i)*dt
      LET yy(i) = yy(i)+vy(i)*dt
   NEXT i
   CALL calcForce
   FOR i=1 TO nMolec
      LET a = 0.5*dt/mas(i)
      LET vx(i) = vx(i)+a*ffx(i)
      LET vy(i) = vy(i)+a*ffy(i)
   NEXT i
   FOR i=1 TO nMolec !periodic condition
      IF xx(i)<0.0 THEN LET xx(i) = xx(i) + xMax
      IF xx(i)>xMax THEN LET xx(i) = xx(i) - xMax
      IF yy(i)<0.0 THEN LET yy(i) = yy(i) + yMax
      IF yy(i)>yMax THEN LET yy(i) = yy(i) - yMax
   NEXT i
   LET sysTime=sysTime+dt
END SUB

EXTERNAL SUB calcForce
   DECLARE EXTERNAL FUNCTION force,boundaryForce
   LET s = 0.5*3.418e-10
   FOR i=1 TO nMolec
      LET ffx(i) = 0.0
      LET ffy(i) = 0.0
   NEXT i
   FOR i=1 TO nMolec-1
      FOR k=1 TO reg(i,0)-1
         LET j = reg(i,k)
         LET xij = xx(i)-xx(j)
         IF xij>0.5*xMax THEN LET xij = xij - xMax  !x-periodic
         IF xij<-0.5*xMax THEN LET xij = xij + xMax !x-periodic
         LET yij = yy(i)-yy(j)
         IF yij>0.5*yMax THEN LET yij = yij - yMax  !y-periodic
         IF yij<-0.5*yMax THEN LET yij = yij + yMax !y-periodic
         LET rij = SQR(xij*xij+yij*yij)
         IF rij<rCutoff THEN
            LET f = force(rij,kind(i),kind(j))
            LET fxij = f*xij/rij
            LET fyij = f*yij/rij
            LET ffx(i) = ffx(i)+fxij
            LET ffy(i) = ffy(i)+fyij
            LET ffx(j) = ffx(j)-fxij
            LET ffy(j) = ffy(j)-fyij
         END IF
      NEXT k
   NEXT i
END SUB

EXTERNAL FUNCTION force(r,ki,kj) !force(r) <-- forceTable - linear interporation
   LET ir = INT(r/hh)
   LET a = r - ir*hh
   LET force = ((hh-a)*forceTable(ki,kj,ir) + a*forceTable(ki,kj,ir+1))/hh
END FUNCTION

EXTERNAL SUB registerNearMolec
   LET rCut = rCutoff+20*2000*dt
   LET rcut2 = rCut*rCut
   FOR i=1 TO nMolec-1
      LET k = 1
      FOR j=i+1 TO nMolec
         LET xij = xx(i)-xx(j)
         IF xij>0.5*xMax THEN LET xij = xij - xMax  !x-periodic
         IF xij<-0.5*xMax THEN LET xij = xij + xMax !x-periodic
         LET yij = yy(i)-yy(j)
         IF yij>0.5*yMax THEN LET yij = yij - yMax  !y-periodic
         IF yij<-0.5*yMax THEN LET yij = yij + yMax !y-periodic
         LET r2 = xij*xij+yij*yij
         IF (r2<rcut2) THEN
            LET reg(i,k) = j
            LET k = k + 1
         END IF
      NEXT j
      LET reg(i,0) = k
   NEXT i
END SUB

EXTERNAL FUNCTION maxNearMolec
   LET mx = 0
   FOR i=1 TO nMolec-1
      IF mx<reg(i,0) THEN LET mx = reg(i,0)
   NEXT i
   LET maxNearMolec = mx-1
END FUNCTION

! ---------- utility

EXTERNAL FUNCTION systemTemprature
   LET kB = 1.38e-23 ! Boltzman's constant (J/K)
   LET ek= 0.0       !kinetic energy (J)
   FOR i=1 TO nMolec
      LET ek = ek + 0.5*mas(i)*(vx(i)^2+vy(i)^2)
   NEXT i
   LET systemTemprature = ek/(nMolec*kB) !2D: E/N=kT, 3D: E/N=(3/2)kT
END FUNCTION

EXTERNAL SUB ajustVelocity(temp)
   DECLARE EXTERNAL FUNCTION systemTemprature
   LET r = sqr(temp/systemTemprature)
   FOR i=1 TO nMolec
      LET vx(i) = r*vx(i)
      LET vy(i) = r*vy(i)
   NEXT i
END SUB

! ---------- draw particles

EXTERNAL SUB drawParticles(tempMode,contTemp,drawMode)
   DECLARE EXTERNAL FUNCTION systemTemprature,maxNearMolec
   DECLARE EXTERNAL sub realSpace,velocitySpace,plotBond
   SET DRAW MODE HIDDEN
   CLEAR
   IF drawMode=0 OR drawMode=1 THEN !--- 0:circle  1:circle+bond
      call plotBond(drawMode)
   ELSEIF drawMode=2 THEN
      call velocitySpace
   END IF

   !--- draw caption
   SET TEXT HEIGHT 10
   SET TEXT COLOR 1 ! black
   LET tmp$ = "adiabatic   constantTemp  "
   PLOT TEXT, AT  50, 70 ,USING "time =#####.## (ps)   temp =####.## (K)":sysTime*1E12,systemTemprature
   PLOT TEXT, AT  50, 55 ,USING molecSTR$(molecKind1)&","&molecSTR$(molecKind2)&"  N =####":nMolec
   PLOT TEXT, AT  50, 40 ,USING "tempMode = ## ":tempMode
   PLOT TEXT, AT 200, 40 :tmp$(tempMode*12+1:tempMode*12+12)
   PLOT TEXT, AT  50, 25 ,USING "controled Temperature =####.# (K)":contTemp
   PLOT TEXT, AT  50, 10 :"periodic molecular dynamics 2D"
   PLOT TEXT, AT  50,480 :"[esc] exit    [0] temp > value "
   PLOT TEXT, AT  50,460 :"[1] dispMode  [9] menu > select"
   SET DRAW MODE EXPLICIT
END SUB

EXTERNAL sub plotBond(drawMode)
   LET boxSize = 300
   LET mag = boxSize/xMax
   LET xp = 100
   LET yp = 100
   SET LINE COLOR 8 !gray !--- box
   PLOT LINES: xp,yp; boxSize+xp,yp; boxSize+xp,boxSize+yp; xp,boxSize+yp; xp,yp
   SET TEXT HEIGHT 6
   SET TEXT COLOR 1 ! black
   PLOT TEXT, AT xp,boxSize+2+yp ,USING  "box size =##.# x ##.# (nm)":xMax*1e9,yMax*1e9
   FOR i=1 TO nMolec
      IF kind(i)=molecKind1 THEN SET LINE COLOR 11 ELSE SET LINE COLOR 12
      DRAW circle WITH SCALE(molecData(kind(i),3)/2.0*mag)*SHIFT(xx(i)*mag+xp,yy(i)*mag+yp)
   NEXT i
   IF drawMode=1 THEN
      FOR i=1 TO nMolec-1
         FOR k=1 TO reg(i,0)-1
            LET j = reg(i,k)
            LET r = SQR((xx(i)-xx(j))*(xx(i)-xx(j))+(yy(i)-yy(j))*(yy(i)-yy(j)))
            LET r0 = 0.5*(molecData(kind(i),3)+molecData(kind(j),3))
            IF r<1.2*r0 THEN
               IF r<0.93*r0 THEN
                  SET LINE COLOR 4 !red
               ELSEIF r<1.03*r0 THEN
                  SET LINE COLOR 3 !green
               ELSEIF r<1.08*r0 THEN
                  SET LINE COLOR 2 !blue
               ELSE
                  SET LINE COLOR 8 !gray
               END IF
               PLOT LINES: xx(i)*mag+xp,yy(i)*mag+yp;xx(j)*mag+xp,yy(j)*mag+yp
            END IF
         NEXT k
      NEXT i
   END IF
END sub

EXTERNAL sub velocitySpace
   LET boxSize = 300
   LET xp = 100
   LET yp = 100
   SET LINE COLOR 1 !black : axis
   PLOT LINES: xp,boxSize/2+yp; boxSize+xp,boxSize/2+yp !vx-axis
   PLOT LINES: boxSize/2+xp,yp; boxSize/2+xp,boxSize+yp !vy-axis
   SET TEXT HEIGHT 6
   SET TEXT COLOR 1 ! black
   PLOT TEXT, AT boxSize+xp,boxSize/2+yp: "vx"
   PLOT TEXT, AT boxSize+xp,boxSize/2-12+yp: "2000m/s"
   PLOT TEXT, AT boxSize/2-12+xp,boxSize+yp: "vy  2000m/s"
   PLOT TEXT, AT boxSize/2-8+xp,boxSize/2-10+yp: "0"
   PLOT TEXT, AT xp,boxSize+8+yp: "velocity space (vx,vy)"
   LET mag = boxSize/4000
   FOR i=1 TO nMolec
      IF kind(i)=molecKind1 THEN SET LINE COLOR 11 ELSE SET LINE COLOR 12
      DRAW circle WITH SCALE(5)*SHIFT(vx(i)*mag+boxSize/2+xp,vy(i)*mag+boxSize/2+yp)
   NEXT i
END sub

END MODULE
 

[018]周期的境界条件下の電子状態(2次元)

 投稿者:mike  投稿日:2017年 6月26日(月)09時33分32秒
返信・引用
   2次元-周期的境界条件下の下で定常状態の電子のエネルギーと電子状態を求めるプログラム
(018harmonicsPSD2D.bas)を公開します。これは、最急降下法による2次元の定常状態の電子状態を求めるプログラム
[011]に周期的境界条件を付けたものです。

表示の説明:
 2次元-周期的境界条件下の自由空間(V(x,y)=0)では基底状態は|0>=constant, E(0)=0 になります。
x、y方向の周期が同じ場合、E1=E2=E3=E4で縮退しているため、|1>..|4>はこの固有空間の中の
任意の直交する関数の組でよいため、形が一意には決まりません。
E5=E6=E7=E8も縮退しています。

試験環境:
本プログラムは十進BASIC 6.6.3.1 / macOS 10.7.5, 十進BASIC Ver 7.8.0 / windows 10 でテストしました。


------------

!
! ========= periodic steepest descent method 2D ==========
!
! 018periodicPSD2D.bas
!   Copyright(C) 2017 Mitsuru Ikeuchi
!   Released under the MIT license ( https://opensource.org/licenses/MIT )
!
! ver 0.0.1   2017.06.26  created
!
OPTION ARITHMETIC NATIVE
DECLARE EXTERNAL SUB psd2d.setInitialCondition,psd2d.SDiteration,psd2d.drawState
DECLARE EXTERNAL FUNCTION INKEY$

LET stateMax = 10  !state 0,1,...,stateMax-1
LET vIndex = 0     !0:free space,  1:2d-crystal
LET iterMax = 5    !5 -iteration in SDiteration()
LET drawMode = 1   !0:draw3DPsi, 1:drawPsi
CALL setInitialCondition(stateMax,vIndex)

DO
   CALL SDiteration(stateMax, iterMax)
   CALL drawState(drawMode,vIndex)
   LET S$=INKEY$
   IF S$="6" THEN CALL setInitialCondition(stateMax,0) !Hermonics
   IF S$="7" THEN CALL setInitialCondition(stateMax,1) !Well
   IF S$="8" THEN LET drawMode = 0 !0:draw3DPsi
   IF S$="9" THEN LET drawMode = 1 !1:drawPsi
LOOP UNTIL S$=CHR$(27)
END

EXTERNAL FUNCTION INKEY$  !from decimal BASIC library
OPTION ARITHMETIC NATIVE
SET ECHO "OFF"
LET S$=""
CHARACTER INPUT NOWAIT: S$
LET INKEY$=S$
END FUNCTION

! ---------- periodic steepest descent method 2D ----------
!
! 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> - dampingFactor*(H-ei)|i>
! (5) orthogonalization { |0>,|1>,..,|i>,.. }  (Gram-Schmidt)
! (6) goto (2)
!
MODULE psd2d
MODULE OPTION ARITHMETIC NATIVE
OPTION BASE 0
PUBLIC SUB setInitialCondition, SDiteration, drawState
SHARE NUMERIC NNx, NNy, dx, dy, iterCount
SHARE NUMERIC cosAx,sinAx,cosAy,sinAy,cx0,cy0,cz0 !--- 3D graphics
SHARE NUMERIC sdEnergy(20)        ! electron state energy
SHARE NUMERIC sdState(20,200,200) ! electron states 0...20 0:ground state
SHARE NUMERIC wrk(200,200)        ! state work space in steepestDescent
SHARE NUMERIC vv(200,200)         ! external potential
SHARE NUMERIC psicol(64,64)       ! MAT PLOT CELL matrix for psi(x,y)
SHARE NUMERIC vcol(64,64)         ! MAT PLOT CELL matrix  for potential V(x,y)
LET NNx = 64                      ! max number of sdState(,NNx,NNy)
LET NNy = 64                      ! max number of sdState(,NNx,NNy)
LET dx = 0.25                     ! (au) x division
LET dy = 0.25                     ! (au) y division
LET iterCount = 0                 ! sd iteration count

! ---------- set initial condition

EXTERNAL SUB setInitialCondition(stateMax,vIndex)   !public
   DECLARE EXTERNAL SUB setInitialState,setPotential,initDraw
   LET iterCount = 0
   CALL setInitialState(stateMax)
   CALL setPotential(vIndex)
   CALL initDraw
   ! set window
   LET xMargin = 60
   LET yMargin = 120
   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=0 TO NNx-1
         FOR j=0 TO NNx-1
            LET sdState(ist,i,j) = RND-0.5
         NEXT j
      NEXT i
      CALL normalizeState(ist)
   NEXT ist
END SUB

EXTERNAL SUB setPotential(vIndex)
   LET x0 = 0.5*NNx*dx
   LET y0 = 0.5*NNy*dy
   FOR i=0 TO NNx-1
      FOR j=0 TO NNy-1
         LET x = i*dx
         LET y = j*dy
         IF vIndex=0 THEN
            LET vv(i,j) = 0
         ELSEIF vIndex=1 THEN
            IF (x-x0)*(y-y0)>0 THEN LET vv(i,j)=10 ELSE LET vv(i,j)=0
         ELSE
            LET vv(i,j)=0
         END IF
      NEXT j
   NEXT i
END SUB

EXTERNAL SUB initDraw !--- set set color pallet and MAT PLOT matrix vcol(,)
   FOR i = 0 TO 50 !--- set color pallet
      SET COLOR MIX( 40+i) 0.02*i,0.02*i,0     ! red for |psi> +
      SET COLOR MIX(100+i) 0,0.02*i,0.02*i     ! blue for |psi> -
      SET COLOR MIX(160+i) 0,0.02*i,0     ! green for V(x)
   NEXT i
   FOR i=0 TO NNx-1 !--- set vcol(,)
      FOR j=0 TO NNy-1
         LET col = 0.02*vv(i,j)
         IF col>1 THEN LET col = 1
         IF col<0 THEN LET col = 0
         LET vcol(i,j) = 160+INT(col*50)
      NEXT j
   NEXT i
END SUB

! ---------- steepest descent iteration

EXTERNAL SUB SDiteration(stateMax, iterMax) !public
   DECLARE EXTERNAL FUNCTION steepestDescent
   DECLARE EXTERNAL SUB GramSchmidt,sortState
   LET damp = 0.01 !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=0 TO NNx-1                          !---  |wrk> = (H-E_ist)|ist>
      LET ip1 = MOD(i+1,NNx)
      LET im1 = MOD(i-1+NNx,NNx)
      FOR j=0 TO NNy-1
         LET jp1 = MOD(j+1,NNy)
         LET jm1 = MOD(j-1+NNy,NNy)
         LET pij = sdState(ist,i,j)
         LET wrk(i,j) = (4*pij-sdState(ist,ip1,j)-sdState(ist,im1,j)-sdState(ist,i,jm1)-sdState(ist,i,jp1))/h2+(vv(i,j)-ei)*pij
      NEXT j
   NEXT i
   FOR i=0 TO NNx-1                          !---  |ist(next)> = |ist> - damp*|wrk>  ( norm(|ist(next)>) <>1 )
      FOR j=0 TO NNy-1
         LET sdState(ist,i,j) = sdState(ist,i,j)-damp*wrk(i,j)
      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=0 TO NNx-1
      LET ip1 = MOD(i+1,NNx)
      LET im1 = MOD(i-1+NNx,NNx)
      FOR j=0 TO NNy-1
         LET jp1 = MOD(j+1,NNy)
         LET jm1 = MOD(j-1+NNy,NNy)
         LET pij = sdState(ist,i,j)
         LET pij = sdState(ist,i,j)
         LET s = s+pij*((4*pij-sdState(ist,ip1,j)-sdState(ist,im1,j)-sdState(ist,i,jm1)-sdState(ist,i,jp1))/h2+vv(i,j)*pij)
         LET sn = sn+pij*pij
      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=0 TO NNx-1
            FOR j=0 TO NNy-1
               LET sdState(istate,i,j) = sdState(istate,i,j) - s*sdState(ist,i,j)
            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
               LET w = sdState(ist,i,j)
               LET sdState(ist,i,j) = sdState(ist+1,i,j)
               LET sdState(ist+1,i,j) = w
            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=0 TO NNx-1
      FOR j=0 TO NNy-1
         LET s = s + sdState(ist,i,j)*sdState(jst,i,j)
      NEXT j
   NEXT i
   LET innerProduct = s*dx*dy
END FUNCTION

EXTERNAL SUB normalizeState(ist)
   LET s = 0
   FOR i=0 TO NNx-1
      FOR j=0 TO NNy-1
         LET s = s + sdState(ist,i,j)*sdState(ist,i,j)*dx*dy
      NEXT j
   NEXT i
   LET a = SQR(1/s)
   FOR i=0 TO NNx-1
      FOR j=0 TO NNy-1
         LET sdState(ist,i,j) = a*sdState(ist,i,j)
      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) !shift*xRotateAx*yRotateAy*(shift^-1)
   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

! ---------- drawState

EXTERNAL SUB drawState(drawMode,vIndex)
   DECLARE EXTERNAL SUB setRotateXYParameters,drawState3D,draw3DPsix6,drawPsix6
   IF vIndex=0 THEN LET v$="free space" ELSE LET v$="2D-crystal"
   SET DRAW MODE HIDDEN
   CLEAR
   IF drawMode=0 THEN
      CALL setRotateXYParameters(-PI/6,-PI/12,NNx/2,NNy/2,0)
      CALL draw3DPsix6
   ELSEIF drawMode=1 THEN
      CALL drawPsix6
   END IF
   !--- caption
   SET TEXT HEIGHT 10
   SET TEXT COLOR 1 ! black
   PLOT TEXT, AT 50,100 ,USING "iteration count =##### ":iterCount
   PLOT TEXT, AT 50, 85 ,USING "E0 =###.###########  E5 =###.###########(au)":sdEnergy(0),sdEnergy(5)
   PLOT TEXT, AT 50, 70 ,USING "E1 =###.###########  E6 =###.###########(au)":sdEnergy(1),sdEnergy(6)
   PLOT TEXT, AT 50, 55 ,USING "E2 =###.###########  E7 =###.###########(au)":sdEnergy(2),sdEnergy(7)
   PLOT TEXT, AT 50, 40 ,USING "E3 =###.###########  E8 =###.###########(au)":sdEnergy(3),sdEnergy(8)
   PLOT TEXT, AT 50, 25 ,USING "E4 =###.###########  E9 =###.###########(au)":sdEnergy(4),sdEnergy(9)
   PLOT TEXT, AT 50, 10 :"periodic steepest descent method 2D"
   PLOT TEXT, AT 50,480 :"potential -> [6] free space [7] 2d-crystal "
   PLOT TEXT, AT 50,465 :"display   -> [8] psi3D      [9] psi  "
   PLOT TEXT, AT 50,450 :"[esc] exit "
   PLOT TEXT, AT 50,420 :"potential : "&v$
   SET DRAW MODE EXPLICIT
END SUB

EXTERNAL SUB draw3DPsix6
   DECLARE EXTERNAL SUB drawState3D
   CALL drawState3D(0,2,0.5, 30,120) !(ist,sc,zMag,xShift,yShift)
   CALL drawState3D(1,2,0.5,180,120)
   CALL drawState3D(2,2,0.5,330,120)
   CALL drawState3D(3,2,0.5, 30,270)
   CALL drawState3D(4,2,0.5,180,270)
   CALL drawState3D(5,2,0.5,330,270)
   SET TEXT HEIGHT 10
   PLOT TEXT, AT  30,130:"|0>"
   PLOT TEXT, AT 180,130:"|1>"
   PLOT TEXT, AT 330,130:"|2>"
   PLOT TEXT, AT  30,280:"|3>"
   PLOT TEXT, AT 180,280:"|4>"
   PLOT TEXT, AT 330,280:"|5>"
END SUB

EXTERNAL SUB drawState3D(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)*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)),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)*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)),xShift,yShift) !(x,y,z)
      NEXT j
      PLOT LINES
   NEXT i
END SUB

EXTERNAL SUB drawPsix6
   DECLARE EXTERNAL SUB drawPsi
   CALL drawPsi(0, 30,135) !(ist,xPos,yPos)
   CALL drawPsi(1,180,135)
   CALL drawPsi(2,330,135)
   CALL drawPsi(3, 30,285)
   CALL drawPsi(4,180,285)
   CALL drawPsi(5,330,285)
   SET TEXT HEIGHT 10
   PLOT TEXT, AT  30,120:"|0>"
   PLOT TEXT, AT 180,120:"|1>"
   PLOT TEXT, AT 330,120:"|2>"
   PLOT TEXT, AT  30,270:"|3>"
   PLOT TEXT, AT 180,270:"|4>"
   PLOT TEXT, AT 330,270:"|5>"
END SUB

EXTERNAL SUB drawPsi(ist,xPos,yPos)  !drawMode=0
   MAT psicol = vcol !--- psicol(,) <- Vcol(,) setted SUB initDraw
   FOR i=0 TO NNx-1 !--- set psicol(,) for MAT PLOT CELLS
      FOR j=0 TO NNy-1
         LET psi = sdState(ist,i,j)
         IF abs(psi)>0.002 THEN
            IF psi>=0 THEN
               LET col = psi*10
               IF col>1 THEN LET col = 1
               IF col<0 THEN LET col = 0
               LET psicol(i,j) = 40+INT(col*50)
            ELSE
               LET col = -psi*10
               IF col>1 THEN LET col = 1
               IF col<0 THEN LET col = 0
               LET psicol(i,j) = 100+INT(col*50)
            END IF
         END IF
      NEXT j
   NEXT i
   MAT PLOT CELLS,IN xPos,yPos; xPos+128,yPos+128 :psicol
END SUB

END MODULE
 

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

 投稿者:SHIRAISHI Kazuo  投稿日:2017年 6月24日(土)08時49分45秒
返信・引用
  十進BASIC32ビット版は,手続き呼び出し時にCPUのSPレジスタを操作して受け渡し場所を確保しているのですが,それがMAC-OSのお気に召さなかったようです。
 

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

 投稿者:SHIRAISHI Kazuo  投稿日:2017年 6月23日(金)23時01分21秒
返信・引用
  ご報告ありがとうございました。

Windows版,Linux版では問題なく,Macだけの問題のようです。

 

実行時内部エラーの報告

 投稿者:mike  投稿日:2017年 6月23日(金)10時41分53秒
返信・引用
  実行時内部エラーが発生しました。

access violation エラーです。

decimalBASIC 6.6.3.1 / MacOS 10.7.5 / MACmini (intel core-i7 2.7GHz 4GB)

再起動後に実行してもエラーが再現します。

下記のプログラムです。
コメントアウトして発生場所を絞り込むと、
SUB setInitialCondition で発生しているものと推定できます。

   !set density
   FOR i=0 TO NNp-1
      LET s = 0.0
      FOR j=0 TO NNp-1
         LET s = s + mass(j)*kernel(distance(i,j), (hh(i)+hh(j))/2.0)
      NEXT j
      LET density(i) = s
   NEXT i

の部分を下記のように変更するとエラーは回避できそうですが

   !set density
   FOR i=0 TO NNp-1
      LET s = 0.0
      FOR j=0 TO NNp-1
         LET r = distance(i,j)
         LET s = s + mass(j)*kernel(r, (hh(i)+hh(j))/2.0)
      NEXT j
      LET density(i) = s
   NEXT i


まだ未完ですが、全体のプログラムを載せます。

-------------

!
! ========= smoothed particle hydrodynamics 2D ==========
!
! 0xxgasSPH2D.bas
!
!    Copyright(C) 2017 Mitsuru Ikeuchi
!    Released under the MIT license ( https://opensource.org/licenses/MIT )
!
!    ver 0.0.0  2017.06.22 created
!
OPTION ARITHMETIC NATIVE
DECLARE EXTERNAL SUB sph2d.setInitialCondition, sph2d.timeEvolution, sph2d.drawParticles

CALL setInitialCondition

LET t0 = TIME
FOR it=1 TO 1000
!CALL timeEvolution(10)
!CALL drawParticles
NEXT it
PRINT TIME - t0
END

! ---------- smoothed particle hydrodynamics 2D ----------
!
! - particle base Lagrangian method
!    Monaghan; "Smoothed Particle Hydrodynamics"
!              Annu. Rev. Astron. Astrophys.1992. 30:543-74
!
!    W(x-xi,h) : kernel weight function (q=|x-xi|/h)
!      = aKernel*(1-1.5*q^2+0.75*q^3)  (q<1)
!      = aKernel*0.25*(2-q)^3  (1<=q<2)
!      = 0  (q>=2)
!
!        aKernel = 2/3/h (1D)
!                = 10/(7Pi)/h^2 (2D)
!                = 1/Pi/h^3 (3D)
!
!    f(x) --> sum[mj/rhoj*f(xj)*W(x-xj,h), j]
!    grad f(x) --> -sum[mj/rhoj*f(xj)*grad W(x-xj,h), j]
!
! - time step
!   (1) registration - set section  (not implemented in this code)
!   (2) set density
!   (3) set pressure
!   (4) Verlet step1 (t = t + dt/2)
!   (5) set acceleration ax ay and power <-- x(t+dt/2),y(t+dt/2)
!   (6) Verlet step2 (t = (t + dt/2) + dt/2)
!   (7) set boundary
!   goto (1)
!
MODULE sph2d
MODULE OPTION ARITHMETIC NATIVE
MODULE OPTION BASE 0
PUBLIC SUB setInitialCondition,timeEvolution,drawParticles
SHARE NUMERIC sysTime, dt, xMax, yMax, hh0, NNp
SHARE NUMERIC xx(5000),yy(5000)  ! (xx(i),yy(i))  : position of i-th particle
SHARE NUMERIC vx(5000),vy(5000)  ! (vx(i),vy(i))  : velocity of i-th particle
SHARE NUMERIC ax(5000),ay(5000)  ! (ax(i),ay(i))  : total force/mass of i-th particle
SHARE NUMERIC hh(5000),mass(5000),cp(5000)
SHARE NUMERIC density(5000),pressure(5000),energy(5000),power(5000)
LET sysTime = 0.0  ! (s) system time
LET dt = 0.0001    ! (s) time step
LET xMax = 5.0     ! (m) x-Box size
LET yMax = 5.0     ! (m) y-Box size
LET hh0 = 0.25     ! (m) h of weight function
LET NNp = 400      ! number of smoothed particle

! ---------- set initial condition

EXTERNAL SUB setInitialCondition
!  i,j, s
   DECLARE EXTERNAL FUNCTION kernel,distance
   RANDOMIZE
   LET sysTime = 0.0
   !set smoothed particles
   FOR i=0 TO NNp-1
      LET mass(i) = (0.029/0.0224)/4.0 !air
      LET cp(i) = 1000.0 !air
      LET xx(i) = 0.1*MOD(i,20)+0.5
      LET yy(i) = 0.1*INT(i/20)+1.0
      LET vx(i) = 0.0
      LET vy(i) = 0.0
      LET ax(i) = 0.0
      LET ay(i) = 0.0
      LET hh(i) = hh0
   NEXT i
   ! set energy and power
   FOR i=0 TO NNp-1
      LET energy(i) = mass(i)*cp(i)*300 ! 300=temp(K)
      LET power(i) = 0.0
   NEXT i
   !set density
   FOR i=0 TO NNp-1
      LET s = 0.0
      FOR j=0 TO NNp-1
         LET s = s + mass(j)*kernel(distance(i,j), (hh(i)+hh(j))/2.0)
      NEXT j
      LET density(i) = s
   NEXT i
   ! set window
   SET WINDOW 0,500, 0,500
END SUB

! ---------- time evolution

EXTERNAL SUB timeEvolution(nCalc)
!  i
   DECLARE EXTERNAL SUB timeStep
   FOR i=0 TO NNp-1
      CALL timeStep
      LET sysTime = sysTime + dt
   NEXT i
END SUB

EXTERNAL SUB timeStep
!  i,j, dtv2,s,rr
   DECLARE EXTERNAL SUB accCalc
   DECLARE EXTERNAL FUNCTION kernel,distance
   LET dtv2 = dt/2.0

   !(2) set density
   FOR i=0 TO NNp-1
      LET s = 0.0
      FOR j=0 TO NNp-1
         LET s = s + mass(j)*kernel(distance(i,j), (hh(i)+hh(j))/2.0)
      Next j
      LET density(i) = s
   NEXT i

   !(3) SET pressure
   FOR i=0 TO NNp-1
      LET pressure(i) = 8.31*density(i)*(energy(i)/(mass(i)*cp(i)))
   NEXT i

   !(4) Verlet step1 (t = t + dt/2)
   FOR i=0 TO NNp-1
      LET vx(i) = vx(i) + dtv2*ax(i)
      LET vy(i) = vy(i) + dtv2*ay(i)
      LET xx(i) = xx(i) + vx(i)*dt
      LET yy(i) = yy(i) + vy(i)*dt
   NEXT i
   !(5) set acceleration ax[] ay[] and power[]
   CALL accCalc
   !(6) Verlet step2 (t = t + dt/2 + dt/2)
   FOR i=0 TO NNp-1
      LET vx(i) = vx(i) + dtv2*ax(i)
      LET vy(i) = vy(i) + dtv2*ay(i)
      LET energy(i) = energy(i) + power(i)*dt
   NEXT i

   !(7) boundary
   LET rr = 1.0
   FOR i=0 TO NNp-1
      if (xx(i) < 0.0) then
         LET xx(i) = 0.0
         LET vx(i) = -rr*vx(i)
         LET vy(i) = rr*vy(i)
      END IF
      if (xx(i) > xMax) then
         LET xx(i) = xMax
         LET vx(i) = -rr*vx(i)
         LET vy(i) = rr*vy(i)
      end if
      if (yy(i) < 0.0) then
         LET yy(i) = 0.0
         LET vx(i) = rr*vx(i)
         LET vy(i) = -rr*vy(i)
      END IF
      if (yy(i) > yMax) then
         LET yy(i) = yMax
         LET vx(i) = rr*vx(i)
         LET vy(i) = -rr*vy(i)
      END IF
   NEXT i
END SUB

! --- (5) set acceleration ax[] ay[] and power[]

EXTERNAL SUB accCalc
!  i,j, ai,aj,b,rij,gradW,gradWx,gradWy
   DECLARE EXTERNAL FUNCTION dwwvdr
   FOR i=0 TO NNp-1
      LET ax(i) = 0.0
      LET ay(i) = 0.0
      LET power(i) = 0.0
      LET ai = pressure(i)/(density(i)*density(i))
      FOR j=0 TO NNp-1
         LET rij = SQR((xx(i)-xx(j))*(xx(i)-xx(j))+(yy(i)-yy(j))*(yy(i)-yy(j)))
         IF (rij>0) THEN
            LET aj = pressure(j)/(density(j)*density(j))
            LET b = mass(j)*(ai+aj)
            LET gradW = dwwvdr(rij, (hh(i)+hh(j))/2.0)
            LET gradWx = gradW*(xx(i)-xx(j))/rij
            LET gradWy = gradW*(yy(i)-yy(j))/rij
            LET ax(i) = ax(i) - b*gradWx
            LET ay(i) = ay(i) - b*gradWy
            LET power(i) = power(i) + 0.5*b*((vx(i)-vx(j))*gradWx+(vy(i)-vy(j))*gradWy)
         END IF
      NEXT j
   NEXT i
END SUB

! --- smoothed particle

EXTERNAL FUNCTION distance(i,j)
!  x,y
   LET x = xx(i)-xx(j)
   LET y = yy(i)-yy(j)
   LET distance = SQR(x*x+y*y)
END FUNCTION

EXTERNAL FUNCTION kernel(r,h)
!  a,q,ret
   LET ret = 0.0
   LET q = r/h
   LET a = 10.0/(7.0*3.141592)/(h*h)
   if (q<1.0) then
      LET ret = a*(1.0-1.5*q*q+0.75*q*q*q)
   ELSEIF (q<2.0) THEN
      LET ret = a*0.25*(2.0-q)*(2.0-q)*(2.0-q)
   END IF
   LET kernel = ret
END FUNCTION

EXTERNAL FUNCTION dwwvdr(r,h)
!  a,q,ret
   LET ret = 0.0
   LET q = r/h
   LET a = 10.0/(7.0*3.141592)/(h*h*h)
   if (q<1.0) then
      LET ret = a*(-3.0*q+2.25*q*q)
   elseif (q<=2.0) then
      LET ret = -a*0.75*(2.0-q)*(2.0-q)
   END IF
   LET dwwvdr = ret
END FUNCTION

! ---------- utility



! ---------- draw particles

EXTERNAL SUB drawParticles
   DECLARE EXTERNAL FUNCTION systemTemprature
   LET boxsize = 300
   LET sc = boxsize/xMax
   LET xp = 100
   LET yp = 80
   LET vScate = 100*dt
   LET fScale = 1000*dt*dt/mass(0)

   SET DRAW MODE HIDDEN
   CLEAR
   SET LINE COLOR 1 ! black : wall
   PLOT LINES: xp,yp; xp+boxSize,yp; boxSize+xp,boxSize+yp; xp,boxSize+yp; xp,yp
   !draw Balls, velocity and force
   FOR i=1 TO nMolec
      SET LINE COLOR 2 ! blue : molecule
      DRAW circle WITH SCALE(sigma/2.0*sc)*SHIFT(xx(i)*sc+xp,yy(i)*sc+yp)
      SET LINE COLOR 4 ! red : velocity
      PLOT LINES: xx(i)*sc+xp,yy(i)*sc+yp;
      PLOT LINES: (xx(i)+vx(i)*vScate)*sc+xp,(yy(i)+vy(i)*vScate)*sc+yp
      !SET LINE COLOR 1 ! black : force
      !PLOT LINES: xx(i)*sc+xp,yy(i)*sc+yp;
      !PLOT LINES: (xx(i)+ffx(i)*fScale)*sc+xp,(yy(i)+ffy(i)*fScale)*sc+yp
   NEXT i
   SET TEXT HEIGHT 6
   SET LINE COLOR 4 ! red : velocity
   PLOT LINES: xp+0.6*boxsize,yp+boxsize+25;xp+0.6*boxsize+50,yp+boxsize+25
   SET TEXT COLOR 4 ! red
   PLOT TEXT, AT xp+0.6*boxsize+55,yp+boxsize+22: "velosity"
   SET LINE COLOR 1 ! black : force
   PLOT LINES: xp+0.6*boxsize,yp+boxsize+10;xp+0.6*boxsize+50,yp+boxsize+10
   SET TEXT COLOR 1 ! black
   PLOT TEXT, AT xp+0.6*boxsize+55,yp+boxsize+7: "force"

   SET TEXT HEIGHT 10
   SET TEXT COLOR 1 ! black
   PLOT TEXT, AT 50, 50 ,USING "time =####.#### ":sysTime
   PLOT TEXT, AT 50, 35 ,USING  "N =####":NNp
   PLOT TEXT, AT 50, 20 :"Ar in the box (2-dimensional molecular dynamics)"
   SET TEXT HEIGHT 6
   PLOT TEXT, AT xp,boxSize+2+yp ,USING  "box size =##.# x ##.# (m)":xMax,yMax
   SET DRAW MODE EXPLICIT
END SUB

END MODULE
 

[017]異種金属結晶の接触

 投稿者:mike  投稿日:2017年 6月19日(月)10時12分40秒
返信・引用
  Morseポテンシャルによる異種金属結晶の接触の分子動力学法プログラム017mixMMD2D.basを公開します。
異なった種類のナノ結晶を接触させます。接触すると表面エネルギーを解放し、温度が高くなります。
一般に表面エネルギー(Morseポテンシャルの場合は解離エネルギーに相当)の大きい方の表面に
小さい方の原子が覆う形が安定になります。

表示の説明:
鉄原子は暗い水色、アルミニウム原子は茶色の円であらわします。
表示は[D]のキーを押すことで変更できます。

試験環境:
本プログラムは十進BASIC 6.6.3.1 / macOS 10.7.5, 十進BASIC Ver 7.8.0 / windows 10 でテストしました。


---------------------


!
! ========= molecular dynamics 2D - Morse potential ==========
!
! 017mixMMD2D.bas
!   Copyright(C) 2017 Mitsuru Ikeuchi
!   Released under the MIT license ( https://opensource.org/licenses/MIT )
!
! ver 0.0.0  2017.04.22  created
! ver 0.0.1  2017.06.11  bug fixed
!
OPTION ARITHMETIC NATIVE
DECLARE EXTERNAL SUB mmd2d.setInitialCondition, mmd2d.moveParticles, mmd2d.drawParticles
DECLARE EXTERNAL FUNCTION INKEY$
LET tempMode = 0    !tempMode: 0:adiabatic  1:constant temperature
LET contTemp = 300  !contTemp: controled constant temperature(K)
LET drawMode = 0    !drawMode: 0:bond  1:circle+bond 2:velocitySpace
DIM menu$(8)
LET menu$(1) = "Fe Al"
LET menu$(2) = "Fe Mo"
LET menu$(3) = "Fe Hg"
LET menu$(4) = "W Pb"
LET menu$(5) = "Ni Cr"
LET menu$(6) = "Ag Cu"
LET menu$(7) = "Al Ar"
LET menu$(8) = "continue"

!setInitialCondition(kind1,kind2,boxSizeInNM,contTemp)
!molKind: 0:W 1:Mo 2:Cr 3:Fe 4:Ni 5:Al 6:Pb 7:Cu 8:Ag 17:Ar 20:Hg
CALL setInitialCondition(3,5,6.0,contTemp)

DO
   CALL moveParticles(tempMode,contTemp)
   CALL drawParticles(tempMode,contTemp,drawMode)
   LET S$=INKEY$
   IF S$="0" THEN
      LOCATE VALUE (1),RANGE 0 TO 4000 : temp
      IF temp>1 THEN
         LET tempMode = 1
         LET contTemp = temp
      ELSE
         LET tempMode = 0
      END IF
   ELSEIF S$="D" OR S$="d" THEN
      LET drawMode = MOD(drawMode+1,3)
   ELSEIF S$="9" THEN
      LOCATE CHOICE (menu$) :nmenu
      IF nmenu=1 THEN      !--- Fe Al box=6x6nm
         CALL setInitialCondition(3,5,6.0,contTemp)
      ELSEIF nmenu=2 THEN  !--- Fe Mo
         CALL setInitialCondition(3,1,6.0,contTemp)
      ELSEIF nmenu=3 THEN  !--- Fe Hg
         CALL setInitialCondition(3,20,6.0,contTemp)
      ELSEIF nmenu=4 THEN  !--- W Pb
         CALL setInitialCondition(0,6,6.0,contTemp)
      ELSEIF nmenu=5 THEN  !--- Ni Cr
         CALL setInitialCondition(4,2,6.0,contTemp)
      ELSEIF nmenu=6 THEN  !--- Ag Cu
         CALL setInitialCondition(8,7,6.0,contTemp)
      ELSEIF nmenu=7 THEN  !--- Al Ar
         CALL setInitialCondition(5,17,6.0,contTemp)
      ELSEIF nmenu=8 THEN  !contine
      END IF
   END IF
LOOP UNTIL S$=CHR$(27)
END


EXTERNAL FUNCTION INKEY$  !from decimalBASIC library
OPTION ARITHMETIC NATIVE
SET ECHO "OFF"
LET S$=""
CHARACTER INPUT NOWAIT: S$
LET INKEY$=S$
END FUNCTION

! ---------- molecular dynamics 2D - Morse potential ----------
!
!  method: velocity Verlet algorithm
!    (1) vi = vi + (Fi/mi)*(0.5dt)
!    (2) ri = ri + vi*dt
!    (3) calculation Fi <- {r1,r2,...,rn} Fi=sum(Fij,j=1 to n),Fij=F(ri-rj)
!    (4) vi = vi + (Fi/mi)*(0.5dt)
!    (6) goto (1)
!  potential: Morse V(r) = D*((1-EXP(-A*(r-r0)))^2-1)
!                        = D*(EXP(-2*A*(r-r0))-2*EXP(-A*(r-r0)))
!    (D:dissociation energy, r0:bond length, A:width parameter { A=SQR(k/(2*D)) }
!             force F(r) = -dV(r)/dr
!                        = 2*D*A*y*(y-1), y=EXP(-A*(r-r0))
MODULE mmd2d
MODULE OPTION ARITHMETIC NATIVE
PUBLIC SUB setInitialCondition !(molKind,boxSizeInNM,xtalSizeInNM,contTemp)
PUBLIC SUB moveParticles !(tempMode,contTemp)
PUBLIC SUB drawParticles !(drawMode:0:realSpace 1:velocitySpace)
SHARE NUMERIC molecKind1,molecKind2, sysTime, dt, nMolec, xMax, yMax, rCutoff, hh
SHARE NUMERIC xx(500),yy(500)           ! (xx(i),yy(i))  : position of i-th particle
SHARE NUMERIC vx(500),vy(500)           ! (vx(i),vy(i))  : velocity of i-th particle
SHARE NUMERIC ffx(500),ffy(500)         ! (ffx(i),ffy(i)): total force of i-th particle
SHARE NUMERIC kind(500),mas(500)        ! kind(i),mas(i) : molec kind, mass of i-th particle
SHARE NUMERIC reg(500,0 TO 100)         ! register near molec reg(i,0):number of near i-th molec
SHARE NUMERIC molecData(0 TO 20,0 TO 3) ! molecule 0:mass, 1:epsilon, 2:sigma, 3:dt
SHARE STRING molecSTR$(0 TO 20)         ! molec string
SHARE NUMERIC forceTable(0 TO 20, 0 TO 20,0 TO 1001) ! force table

LET molecKind1 = 3         ! molecule kind1 - 3:Fe
LET molecKind2 = 1         ! molecule kind2 - 1:Mo
LET sysTime = 0.0          ! system time (s) in the module
LET dt = 5.0*1.0e-15       ! time step (s)
LET nMolec = 100           ! number of particles
LET xMax = 6.0E-9          ! x-Box size (m)
LET yMax = 6.0E-9          ! y-Box size (m)
LET rCutoff = 1.0e-9       ! force cutoff radius (m)
LET hh = 1e-12             ! force table step (m)
LET molecData(0,0)=0       ! if molecData(0,0)=0 then set molecData(,)

! ---------- set initial condition

EXTERNAL SUB setInitialCondition(kind1,kind2,boxSizeInNM,contTemp)
   DECLARE EXTERNAL SUB setMolecData,setMolecules,ajustVelocity,setForceTable
   DECLARE EXTERNAL FUNCTION setCrystalBlock
   RANDOMIZE
   CALL setMolecData
   LET molecKind1 = kind1
   LET molecKind2 = kind2
   LET sysTime = 0.0
   LET xMax = boxSizeInNM*1.0e-9
   LET yMax = boxSizeInNM*1.0e-9
   ! set particles
   LET s = 0.1*xMax
   LET xtalSize = 0.35*xMax
   LET ss = 0.55*xMax
   LET ii = setCrystalBlock(1, kind1, s, s, xtalSize, xtalSize, PI/4)
   LET ii = setCrystalBlock(ii+1, kind2, s, ss, xtalSize, xtalSize, 0)
   LET ii = setCrystalBlock(ii+1, kind2, ss, s, xtalSize, xtalSize, 0)
   LET nMolec = setCrystalBlock(ii+1, kind1, ss, ss, xtalSize, xtalSize, 0)
   CALL ajustVelocity(contTemp)
   LET rCutoff = MIN(1.0e-9, 3.0*MAX(molecData(kind1,3),molecData(kind2,3)))
   CALL setForceTable
   SET COLOR MIX(240) 1.0,0.2,0.2 !for bond direction color
   SET COLOR MIX(241) 0.6,0.6,0.0
   SET COLOR MIX(242) 0.0,1.0,0.0
   SET COLOR MIX(243) 0.0,0.6,0.6
   SET COLOR MIX(244) 0.2,0.2,1.0
   SET COLOR MIX(245) 0.8,0.0,0.8
   ! set window
   SET WINDOW 0,500,0,500
END SUB

EXTERNAL SUB setMolecData
! Morse potential data
!       0:mass(AU) 1:D(eV)  2:A(m^-1)  3:r0(m)
   DATA   183.85 , 0.9906, 1.4116e10, 3.032e-10 !  0 W
   DATA    95.94 , 0.8032, 1.5079e10, 2.976e-10 !  1 Mo
   DATA    51.996, 0.4414, 1.5721e10, 2.754e-10 !  2 Cr
   DATA    55.847, 0.4174, 1.3885e10, 2.845e-10 !  3 Fe
   DATA    58.71 , 0.4205, 1.4199e10, 2.780e-10 !  4 Ni
   DATA    26.98 , 0.2703, 1.1646e10, 3.253e-10 !  5 Al
   DATA   207.19 , 0.2348, 1.1836e10, 3.733e-10 !  6 Pb
   DATA    63.54 , 0.3429, 1.3588e10, 2.866e-10 !  7 Cu
   DATA   107.87 , 0.3323, 1.3690e10, 3.115e-10 !  8 Ag
   DATA    40.08 , 0.1623, 0.8054e10, 4.569e-10 !  9 Ca
   DATA    87.62 , 0.1513, 0.7878e10, 4.988e-10 ! 10 Sr
   DATA   137.34 , 0.1416, 0.6570e10, 5.373e-10 ! 11 Ba
   DATA    22.99 , 0.0633, 0.5900e10, 5.336e-10 ! 12 Na
   DATA    39.102, 0.0542, 0.4977e10, 6.369e-10 ! 13 K
   DATA    85.47 , 0.0464, 0.4298e10, 7.207e-10 ! 14 Rb
   DATA   132.905, 0.0449, 0.4157e10, 7.557e-10 ! 15 Cs
   DATA    20.183, 0.0031, 1.6500e10, 3.076e-10 ! 16 Ne
   DATA    39.948, 0.0104, 1.3400e10, 3.816e-10 ! 17 Ar
   DATA    83.80 , 0.0141, 1.2500e10, 4.097e-10 ! 18 Kr
   DATA   131.30 , 0.0200, 1.2400e10, 4.467e-10 ! 19 Xe
   DATA   200.59 , 0.0734, 1.4900e10, 3.255e-10 ! 20 Hg
   DATA "W" ,"Mo","Cr","Fe","Ni","Al","Pb","Cu","Ag","Ca","Sr"
   DATA "Ba","Na","K" ,"Rb","Cs","Ne","Ar","Kr","Xe","Hg"

   IF molecData(0,0)=0 THEN
      MAT READ molecData
      FOR i=0 TO 20
         LET molecData(i,0) = molecData(i,0)*1.661e-27 !mass(AU)--> (kg)
         LET molecData(i,1) = molecData(i,1)*1.602e-19 !eps(eV) --> (J)
      NEXT i
      MAT READ molecSTR$
   END IF
END SUB

EXTERNAL SUB setGas(nMolecule)
   DECLARE EXTERNAL SUB ajustVelocity
   LET r0 = molecData(knd,3)
   FOR j=1 TO nMolecule
      LET loopCount = 0
      DO
         LET xx(j) = (xMax-2*r0)*RND + r0
         LET yy(j) = (yMax-2*r0)*RND + r0
         FOR i=1 TO j-1
            IF (xx(i)-xx(j))^2+(yy(i)-yy(j))^2 < 2*r0^2 THEN EXIT FOR
         NEXT i
         LET loopCount = loopCount + 1
         IF loopCount>1000 THEN EXIT DO
      LOOP UNTIL i>=j
      IF loopCount>1000 THEN
         LET nMolec = j - 1
         EXIT FOR
      END IF
   NEXT j
   FOR i=1 TO nMolec
      LET vx(i) = 200.0*(RND+RND+RND+RND+RND+RND-3)
      LET vy(i) = 200.0*(RND+RND+RND+RND+RND+RND-3)
      LET ffx(i) = 0.0
      LET ffy(i) = 0.0
   NEXT i
END SUB

EXTERNAL FUNCTION setCrystalBlock(ii, knd, x0, y0, xLen, yLen, theta)
   DECLARE EXTERNAL SUB setParticle
   LET iip = ii
   LET a = 0.98*molecData(knd,3) !r0 of knd
   LET b = 0.866025*a
   LET leng = xLen
   IF (leng<yLen) THEN LET leng = yLen
   LET leng = 1.5*leng
   LET nx = INT(leng/b) + 1
   LET ny = INT(leng/a) + 1
   LET sth = SIN(theta)
   LET cth = COS(theta)
   FOR i=0 TO nx-1
      LET x = b*i - leng/2.0
      FOR j=0 TO ny-1
         LET y = a*j - leng/2.0
         IF MOD(i,2)=1 THEN LET y = y + 0.5*a
         LET xp = x0 + xLen/2.0 + cth*x - sth*y
         LET yp = y0 + yLen/2.0 + sth*x + cth*y
         IF (xp>=x0 AND xp<=x0+xLen AND yp>=y0 AND yp<=y0+yLen) THEN
            CALL setParticle(iip, knd, xp, yp)
            LET iip = iip + 1
         END IF
      NEXT j
   NEXT i
   LET setCrystalBlock = iip - 1
end function

EXTERNAL SUB setParticle(i, knd, x, y)
   LET xx(i) = x
   LET yy(i) = y
   LET vx(i) = 200.0*(RND+RND+RND+RND+RND+RND-3)
   LET vy(i) = 200.0*(RND+RND+RND+RND+RND+RND-3)
   LET ffx(i) = 0.0
   LET ffy(i) = 0.0
   LET mas(i) = molecData(knd,0)
   LET kind(i) = knd
END SUB

EXTERNAL SUB setForceTable
   DECLARE EXTERNAL  FUNCTION cutoff
   FOR ki=0 TO 20
      FOR kj=0 TO 20
         LET dd = SQR(molecData(ki,1)*molecData(kj,1))
         LET aa = 0.5*(molecData(ki,2)+molecData(kj,2))
         LET r0 = 0.5*(molecData(ki,3)+molecData(kj,3))
         FOR ir=10 TO 1001
            LET r = ir*hh
            LET y = EXP(-aa*(r-r0))
            LET forceTable(ki,kj,ir) = cutoff(r)*2.0*dd*aa*y*(y-1)
         NEXT ir
         FOR ir=0 TO 9
            LET forceTable(ki,kj,ir) = forceTable(ki,kj,10)
         NEXT ir
      NEXT kj
   NEXT ki
END SUB

EXTERNAL FUNCTION cutoff(r)
   IF r>0 AND r<0.8*rCutoff THEN
      LET ret = 1
   ELSEIF r>=0.8*rCutoff AND r<rCutoff THEN
      LET ret = 0.5+0.5*COS(PI*(r-0.8*rCutoff)/(0.2*rCutoff))
   else
      LET ret = 0
   END IF
   LET cutoff = ret
END FUNCTION

! ---------- move particles

EXTERNAL SUB moveParticles(tempMode,contTemp) !tempMode 0:adiabatic  1:constant-temp
   DECLARE EXTERNAL SUB moveParticlesDT,ajustVelocity,registerNearMolec
   IF (tempMode=1) THEN CALL ajustVelocity(contTemp)
   CALL registerNearMolec
   FOR i=1 TO 20
      CALL moveParticlesDT
   NEXT i
END SUB

EXTERNAL SUB moveParticlesDT ! velocity Verlet method
   DECLARE EXTERNAL SUB calcForce
   FOR i=1 TO nMolec
      LET a = 0.5*dt/mas(i)
      LET vx(i) = vx(i)+a*ffx(i)
      LET vy(i) = vy(i)+a*ffy(i)
      LET xx(i) = xx(i)+vx(i)*dt
      LET yy(i) = yy(i)+vy(i)*dt
   NEXT i
   CALL calcForce
   FOR i=1 TO nMolec
      LET a = 0.5*dt/mas(i)
      LET vx(i) = vx(i)+a*ffx(i)
      LET vy(i) = vy(i)+a*ffy(i)
   NEXT i
   LET sysTime=sysTime+dt
END SUB

EXTERNAL SUB calcForce
   DECLARE EXTERNAL FUNCTION force,boundaryForce
   LET s = 0.5*3.418e-10
   FOR i=1 TO nMolec
      LET ffx(i) = 0.0
      LET ffy(i) = 0.0
   NEXT i
   FOR i=1 TO nMolec-1
      FOR k=1 TO reg(i,0)-1
         LET j = reg(i,k)
         LET xij = xx(i)-xx(j)
         LET yij = yy(i)-yy(j)
         LET rij = SQR(xij*xij+yij*yij)
         IF rij<rCutoff THEN
            LET f = force(rij,kind(i),kind(j))
            LET fxij = f*xij/rij
            LET fyij = f*yij/rij
            LET ffx(i) = ffx(i)+fxij
            LET ffy(i) = ffy(i)+fyij
            LET ffx(j) = ffx(j)-fxij
            LET ffy(j) = ffy(j)-fyij
         END IF
      NEXT k
   NEXT i
   FOR i=1 TO nMolec  ! boundary force
      LET ffx(i) = ffx(i)+boundaryForce(xx(i)+s)+boundaryForce(xx(i)-xMax-s)
      LET ffy(i) = ffy(i)+boundaryForce(yy(i)+s)+boundaryForce(yy(i)-yMax-s)
   NEXT i
END SUB

EXTERNAL FUNCTION force(r,ki,kj) !force(r) <-- forceTable - linear interporation
   LET ir = INT(r/hh)
   LET a = r - ir*hh
   LET force = ((hh-a)*forceTable(ki,kj,ir) + a*forceTable(ki,kj,ir+1))/hh
END FUNCTION

EXTERNAL FUNCTION boundaryForce(r)
   LET r6 = (3.418e-10/r)^6
   LET boundaryForce = (24.0*(0.5*1.711e-21)*r6*(2.0*r6-1.0)/r)
END FUNCTION

EXTERNAL SUB registerNearMolec
   LET rCut = rCutoff+20*2000*dt
   LET rcut2 = rCut*rCut
   FOR i=1 TO nMolec-1
      LET k = 1
      FOR j=i+1 TO nMolec
         LET r2 = (xx(i)-xx(j))*(xx(i)-xx(j))+(yy(i)-yy(j))*(yy(i)-yy(j))
         IF (r2<rcut2) THEN
            LET reg(i,k) = j
            LET k = k + 1
         END IF
      NEXT j
      LET reg(i,0) = k
   NEXT i
END SUB

EXTERNAL FUNCTION maxNearMolec
   LET mx = 0
   FOR i=1 TO nMolec-1
      IF mx<reg(i,0) THEN LET mx = reg(i,0)
   NEXT i
   LET maxNearMolec = mx-1
END FUNCTION

! ---------- utility

EXTERNAL FUNCTION systemTemprature
   LET kB = 1.38e-23 ! Boltzman's constant (J/K)
   LET ek= 0.0       !kinetic energy (J)
   FOR i=1 TO nMolec
      LET ek = ek + 0.5*mas(i)*(vx(i)^2+vy(i)^2)
   NEXT i
   LET systemTemprature = ek/(nMolec*kB) !2D: E/N=kT, 3D: E/N=(3/2)kT
END FUNCTION

EXTERNAL SUB ajustVelocity(temp)
   DECLARE EXTERNAL FUNCTION systemTemprature
   LET r = sqr(temp/systemTemprature)
   FOR i=1 TO nMolec
      LET vx(i) = r*vx(i)
      LET vy(i) = r*vy(i)
   NEXT i
END SUB

! ---------- draw particles

EXTERNAL SUB drawParticles(tempMode,contTemp,drawMode)
   DECLARE EXTERNAL FUNCTION systemTemprature,maxNearMolec
   DECLARE EXTERNAL sub realSpace,velocitySpace,plotBond
   SET DRAW MODE HIDDEN
   CLEAR
   IF drawMode=0 OR drawMode=1 THEN !--- 0:circle  1:circle+bond
      call plotBond(drawMode)
   ELSEIF drawMode=2 THEN
      call velocitySpace
   END IF
   !--- draw caption
   SET TEXT HEIGHT 10
   SET TEXT COLOR 1 ! black
   LET tmp$ = "adiabatic   constantTemp  "
   PLOT TEXT, AT  50, 70 ,USING "time =#####.## (ps)   temp =####.## (K)":sysTime*1E12,systemTemprature
   PLOT TEXT, AT  50, 55 ,USING molecSTR$(molecKind1)&","&molecSTR$(molecKind2)&"  N =####":nMolec
   PLOT TEXT, AT  50, 40 ,USING "tempMode = ## ":tempMode
   PLOT TEXT, AT 200, 40 :tmp$(tempMode*12+1:tempMode*12+12)
   PLOT TEXT, AT  50, 25 ,USING "controled Temperature =####.# (K)":contTemp
   PLOT TEXT, AT  50, 10 :"2-dimensional molecular dynamics"
   PLOT TEXT, AT  50,480 :"[esc] exit    [0] temp > value "
   PLOT TEXT, AT  50,460 :"[D] dispMode  [9] menu > select"
   SET DRAW MODE EXPLICIT
END SUB

EXTERNAL sub plotBond(drawMode)
   LET boxSize = 300
   LET mag = boxSize/xMax
   LET xp = 100
   LET yp = 100
   SET LINE COLOR 1 ! black : !--- box
   PLOT LINES: xp,yp; boxSize+xp,yp; boxSize+xp,boxSize+yp; xp,boxSize+yp; xp,yp
   SET TEXT HEIGHT 6
   SET TEXT COLOR 1 ! black
   PLOT TEXT, AT xp,boxSize+2+yp ,USING  "box size =##.# x ##.# (nm)":xMax*1e9,yMax*1e9
   FOR i=1 TO nMolec
      IF kind(i)=molecKind1 THEN SET LINE COLOR 11 ELSE SET LINE COLOR 12
      DRAW circle WITH SCALE(molecData(kind(i),3)/2.0*mag)*SHIFT(xx(i)*mag+xp,yy(i)*mag+yp)
   NEXT i
   IF drawMode=1 THEN
      FOR i=1 TO nMolec-1
         FOR k=1 TO reg(i,0)-1
            LET j = reg(i,k)
            LET r = SQR((xx(i)-xx(j))*(xx(i)-xx(j))+(yy(i)-yy(j))*(yy(i)-yy(j)))
            LET r0 = 0.5*(molecData(kind(i),3)+molecData(kind(j),3))
            IF r<1.2*r0 THEN
               LET yij = yy(i)-yy(j)
               IF yij=0.0 THEN LET yij = 1e-20
               LET th = 3.0*(ATN((xx(i)-xx(j))/yij)+PI/2.0)/PI
               LET col = 240 + INT((th-INT(th))*6)
               SET LINE COLOR col
               PLOT LINES: xx(i)*mag+xp,yy(i)*mag+yp;xx(j)*mag+xp,yy(j)*mag+yp
            END IF
         NEXT k
      NEXT i
   END IF
END sub

EXTERNAL sub velocitySpace
   LET boxSize = 300
   LET xp = 100
   LET yp = 100
   SET LINE COLOR 1 !black : axis
   PLOT LINES: xp,boxSize/2+yp; boxSize+xp,boxSize/2+yp !vx-axis
   PLOT LINES: boxSize/2+xp,yp; boxSize/2+xp,boxSize+yp !vy-axis
   SET TEXT HEIGHT 6
   SET TEXT COLOR 1 ! black
   PLOT TEXT, AT boxSize+xp,boxSize/2+yp: "vx"
   PLOT TEXT, AT boxSize+xp,boxSize/2-12+yp: "2000m/s"
   PLOT TEXT, AT boxSize/2-12+xp,boxSize+yp: "vy  2000m/s"
   PLOT TEXT, AT boxSize/2-8+xp,boxSize/2-10+yp: "0"
   PLOT TEXT, AT xp,boxSize+8+yp: "velocity space (vx,vy)"
   LET mag = boxSize/4000
   FOR i=1 TO nMolec
      IF kind(i)=molecKind1 THEN SET LINE COLOR 11 ELSE SET LINE COLOR 12
      DRAW circle WITH SCALE(5)*SHIFT(vx(i)*mag+boxSize/2+xp,vy(i)*mag+boxSize/2+yp)
   NEXT i
END sub
END MODULE
 

C/C++ 数学関数

 投稿者:しばっち  投稿日:2017年 6月18日(日)12時42分46秒
返信・引用
  C/C++で定義されている数学関数を定義してみました。

INPUT Z
LET L=LGAMMA(Z)
PRINT L;EXP(L)
END

EXTERNAL FUNCTION GAMMA(Z)
LET X$=PACKDBL$(Z)
IF Z>0 THEN
   LET GAMMA=GAMMA_(X$)
ELSE
   IF FP(Z)=0 THEN
      PRINT "定義域エラー"
      STOP
   END IF
   LET GAMMA=GAMMA_(X$)
END IF
FUNCTION GAMMA_(X$)
   ASSIGN ".\DLL\gamma_.dll","tgamma_",FPU
END FUNCTION
END FUNCTION

EXTERNAL FUNCTION LGAMMA(Z)
LET X$=PACKDBL$(Z)
LET LGAMMA=LGAMMA_(X$)

FUNCTION LGAMMA_(X$)
   ASSIGN ".\DLL\gamma_.dll","lgamma_",FPU
END FUNCTION
END FUNCTION

EXTERNAL FUNCTION DIGAMMA(Z)
LET X$=PACKDBL$(Z)
LET DIGAMMA=DIGAMMA_(X$)

FUNCTION DIGAMMA_(X$)
   ASSIGN ".\DLL\gamma_.dll","digamma_",FPU
END FUNCTION
END FUNCTION

EXTERNAL FUNCTION TRIGAMMA(Z)
LET X$=PACKDBL$(Z)
LET TRIGAMMA=TRIGAMMA_(X$)

FUNCTION TRIGAMMA_(X$)
   ASSIGN ".\DLL\gamma_.dll","trigamma_",FPU
END FUNCTION
END FUNCTION

EXTERNAL FUNCTION POLYGAMMA(N,Z)
LET X$=PACKDBL$(Z)
LET POLYGAMMA=POLYGAMMA_(INT(N),X$)

FUNCTION POLYGAMMA_(N,X$)
   ASSIGN ".\DLL\gamma_.dll","polygamma_",FPU
END FUNCTION
END FUNCTION
-------------------------------------------------------------------------------------------
                                        gamma_.cpp

#include <boost/math/special_functions/gamma.hpp>
#include <boost/math/special_functions/digamma.hpp>
#include <boost/math/special_functions/trigamma.hpp>
#include <boost/math/special_functions/polygamma.hpp>
using namespace boost::math;
using namespace std;

extern "C"  __declspec(dllexport)  double  tgamma_(double *a)
{
    return tgamma(*a);
}

extern "C"  __declspec(dllexport)  double  lgamma_(double *a)
{
    return lgamma(*a);
}

extern "C"  __declspec(dllexport)  double  digamma_(double *a)
{
    return digamma(*a);
}

extern "C"  __declspec(dllexport)  double  trigamma_(double *a)
{
    return trigamma(*a);
}

extern "C"  __declspec(dllexport)  double  polygamma_(int n,double *a)
{
    return polygamma(n,*a);
}



※ご注意

DLL側ではエラー処理をしていません。私自身が定義域をよく理解していませんので。
BASIC側で処理するか、範囲を超えないよう注意してください。



また、1000桁モードで使用できる関数も定義してみました。

OPTION ARITHMETIC DECIMAL_HIGH
LET X=1/5
CALL LATAN(1000,STR$(X),RESULT$)
LET S=16*VAL(RESULT$)
LET X=1/239
CALL LATAN(1000,STR$(X),RESULT$)
LET S=S-4*VAL(RESULT$)
PRINT S
PRINT PI
END

EXTERNAL  SUB LATAN(KETA,X$,RESULT$)
OPTION ARITHMETIC DECIMAL_HIGH
LET B$=REPEAT$(CHR$(0),KETA+100)
CALL ATAN1000(KETA,X$,B$)
FOR I=LEN(B$) TO 1 STEP -1
   IF B$(I:I)<="9" AND B$(I:I)>="0" THEN
      EXIT FOR
   END IF
NEXT I
LET RESULT$=B$(1:I)

SUB ATAN1000(KETA,X$,Y$)
   ASSIGN ".\DLL\atan1000.dll","atan1000"
END SUB
END SUB

EXTERNAL  SUB LACOTAN(KETA,X$,RESULT$)
OPTION ARITHMETIC DECIMAL_HIGH
LET B$=REPEAT$(CHR$(0),KETA+100)
CALL ACOTAN1000(KETA,X$,B$)
FOR I=LEN(B$) TO 1 STEP -1
   IF B$(I:I)<="9" AND B$(I:I)>="0" THEN
      EXIT FOR
   END IF
NEXT I
LET RESULT$=B$(1:I)

SUB ACOTAN1000(KETA,X$,Y$)
   ASSIGN ".\DLL\atan1000.dll","acotan1000"
END SUB
END SUB
-------------------------------------------------------------------------------------------
                                         atan1000.cpp

#include <cmath>
#include <string>
#include <boost/multiprecision/mpfr.hpp>

using namespace std;
using namespace boost::multiprecision;

extern "C"  __declspec(dllexport)  void atan1000(int keta,char *a,char *b)
{
    mpfr_float::default_precision(keta);
    mpfr_float x,y;
    string s;
    x.assign(a);
    y=atan(x);
    s=y.str(keta,true);
    strcpy(b,s.c_str());
}

extern "C"  __declspec(dllexport)  void acotan1000(int keta,char *a,char *b)
{
    mpfr_float::default_precision(keta);
    mpfr_float x,y;
    string s;
    x.assign(a);
    y=atan(1/x);
    s=y.str(keta,true);
    strcpy(b,s.c_str());
}


GMP互換ライブラリーを使用しています。

https://gmplib.org/

「mpir.dll」「mpfr.dll」をBASIC.EXEと同じフォルダーに入れてください。

http://mpir.org/
http://www.mpfr.org/


下記よりダウンロードしてください。(math.zip)

ダウンロードパス:shibacchi

http://fast-uploader.com/file/7053312569584/


なお、有効期限は本日より1ヶ月間となります。

P.S
VC++2017でコンパイルした Cファイルは問題なかったのですが、
C++ファイルはBASIC側で「DLLファイルがロードできない」等のエラーが出ました。
仕方なく諦めようかと思いましたが、VC++2010でコンパイルしてみたところ
問題なく動作しました。よってVC++2017とVC++2010でコンパイルしたDLLファイル
が混在しています。
 

レンタル掲示板
/150