Python:分子動力学を理解するために(1)

分子動力学を理解するために、pythonでLJポテンシャルを持つ粒子のコードを書いてみました。pythonで書いた理由は、楽だからです。設計も何も考えずに書くには、スクリプト言語が最適ですね。以下のコードは、周期境界条件やその他諸々の効果が抜けておりおります。また、粒子がとびとびになる。笑。

#! /usr/bin/env python
import numpy as np
import random
import itertools
import scipy.spatial.distance as scipy_dist

class md:
    def __init__(self):
        atom_num=10
        step=1000
        self.step=step
        self.atom_num=atom_num
        self.dt=0.01
        self.xyz=np.zeros((step,atom_num,3))
        self.vxyz=np.zeros((step,atom_num,3))
        self.fxyz=np.zeros((step,atom_num,3))
        self.dist=scipy_dist.euclidean

    def force(self,atom_xyz,index):
        out_force=0
        for i,tmp_xyz in enumerate(atom_xyz):
            if index !=i:
                length=self.dist(atom_xyz[index],atom_xyz[i])
                out_force+=(1/length)**14-10*(1/length)**8*(atom_xyz[i]-atom_xyz[index])
        return out_force
        
    def velet(self,istep):
        for i in range(self.atom_num):
            self.xyz[istep,i,:]=2*self.xyz[istep-1,i,:]-self.xyz[istep-2,i,:]+self.dt**2/3*self.force(self.xyz[istep-1,:,:],i)
            self.vxyz[istep-2,i,:]=(self.xyz[istep-1,i,:]-self.xyz[istep-3,i,:])/self.dt
        
    def start(self):
        for i in range(3,self.step):
            self.velet(i)
            
    def init_set(self):
        for i in range(self.atom_num):
            for j in range(3):
                self.xyz[0,i,j]=random.random()+random.randint(0,10)
                self.vxyz[0,i,j]=random.random()+random.randint(0,10)
        for i in range(self.atom_num):
            self.xyz[1,i,:]=self.xyz[0,i,:]+self.dt/2*self.vxyz[0,i,:]+self.dt**2/6*self.force(self.xyz[0,:,:],i)
            self.xyz[2,i,:]=2*self.xyz[1,i,:]-self.xyz[0,i,:]+self.dt**2/3*self.force(self.xyz[1,:,:],i)
    def show(self):
        for ixyz in self.xyz:
            for i in range(self.atom_num):
                tmp_list=tuple(ixyz[i].tolist())
                print "ATOM      1  CA  MET A   1     %6.3f %6.3f %6.3f  1.00  1.00" % tmp_list
            print "END"
                
if __name__=="__main__":
    test=md()
    test.init_set()
    test.start()
    test.show()