QA@IT

分子動力学(速度ベルレ法)について

396 PV

このときにrx_bとry_b
がrxとryに毎ステップ更新されてしまいます。この現象を解決したいので、どうかご享受お願いします。

コード
rx_b=rx

ry_b=ry

##############################メインループ########################################

for t in range(1,nsteps+1):

t_lis.append(t*h)



#####速度ベルレ法#####

for i in range(N):

    ry_b[i]=ry[i]

    rx[i]=rx[i]+vx[i]*h+(ax[i]*h*h)/2

    ry[i]=ry[i]+vy[i]*h+(ay[i]*h*h)/2

    vx[i]=vx[i]+ax[i]*h/2

    vy[i]=vy[i]+ay[i]*h/2

ax,ay,pot=force(N,rx,ry,ax,ay,lx,ly,pot)

for i in range(N):   

    vx[i]=vx[i]+ax[i]*h/2

    vy[i]=vy[i]+ay[i]*h/2

ax,ay,pot=force(N,rx,ry,ax,ay,lx,ly,pot)

######ピストンの速度ベルレ法#######

py_p=py

pv_p=pv##ピストン速度

py+=pv*h+pa*h*h/2

########x:周期境界#######

q_in=0.0

q_out=0.0

for i in range(N):

    if rx[i]>lx:

        rx[i]-=lx

    elif rx[i]<=0.0:

        rx[i]+=lx

#########熱壁(y)###########

    if ry[i]<=(0.5):                       

        vy_l=vy[i]

        vy[i]=canon_b(temp_l)

        dq=0.5*(vy[i]**2-vy_l**2)

        if dq>=0.0:   

            q_in+=dq

        elif dq<0.0:

            q_out+=dq

        hit_l+=1

        q_in_lis.append(q_in)

        q_out_lis.append(q_out)

    elif ry[i]>(ly-0.5):

        vy_l=vy[i]

        vy[i]=-canon_b(temp_h)

        dq=0.5*(vy[i]**2-vy_l**2)

        if dq<0.0:

            q_out+=dq

        elif dq>=0.0:

            q_in+=dq

        hit_h+=1

        q_in_lis.append(q_in)

        q_out_lis.append(q_out)

q_in_sum+=q_in

q_out_sum+=q_out

q_in_list.append(q_in_sum)

q_out_list.append(q_out_sum)

##########ピストン熱壁########

for i in range(N):

    if ry_b[i]+0.5>=py and ry[i]+0.5<py:

        hit+=1

        vy_p=vy[i]

        vy[i]=canon_b(temp_p)+pv

        pa_m=(vy[i]-vy_p)/(mp*h)

        pa+=pa_m

        dq=(vy[i]**2-vy_p**2)*0.5

        if dq>=0:

            q_in+=dq

        else:

            q_out+=dq



    elif ry_b[i]-0.5<py and ry[i]-0.5>=py:

        hit+=1

        vy_p=vy[i]

        vy[i]=-canon_b(temp_p)+pv

        pa_p=(vy[i]-vy_p)/(mp*h)

        pa+=pa_p

        dq=(vy[i]**2-vy_p**2)*0.5

        if dq>=0:

            q_in+=dq

        else:

            q_out+=dq



pv+=pa*h/2

py_lis.append(py)
ウォッチ

この質問への回答やコメントをメールでお知らせします。